subroutine r9sifg (x, f, g) c december 1980 edition. w. fullerton, bell labs dimension f1cs(20), f2cs(29), g1cs(21), g2cs(34) external alog, csevl, exp, inits, r1mach, sqrt c c series for f1 on the interval 2.00000e-02 to 6.25000e-02 c with weighted error 2.82e-17 c log weighted error 16.55 c significant figures required 15.36 c decimal places required 17.20 c data f1 cs( 1) / -0.1191081969 051363610 e0/ data f1 cs( 2) / -0.0247823144 996236248 e0/ data f1 cs( 3) / 0.0011910281 453357821 e0/ data f1 cs( 4) / -0.0000927027 714388562 e0/ data f1 cs( 5) / 0.0000093373 141568271 e0/ data f1 cs( 6) / -0.0000011058 287820557 e0/ data f1 cs( 7) / 0.0000001464 772071460 e0/ data f1 cs( 8) / -0.0000000210 694496288 e0/ data f1 cs( 9) / 0.0000000032 293492367 e0/ data f1 cs( 10) / -0.0000000005 206529618 e0/ data f1 cs( 11) / 0.0000000000 874878885 e0/ data f1 cs( 12) / -0.0000000000 152176187 e0/ data f1 cs( 13) / 0.0000000000 027257192 e0/ data f1 cs( 14) / -0.0000000000 005007053 e0/ data f1 cs( 15) / 0.0000000000 000940241 e0/ data f1 cs( 16) / -0.0000000000 000180014 e0/ data f1 cs( 17) / 0.0000000000 000035063 e0/ data f1 cs( 18) / -0.0000000000 000006935 e0/ data f1 cs( 19) / 0.0000000000 000001391 e0/ data f1 cs( 20) / -0.0000000000 000000282 e0/ c c series for f2 on the interval 0.00000e+00 to 2.00000e-02 c with weighted error 4.32e-17 c log weighted error 16.36 c significant figures required 14.75 c decimal places required 17.10 c data f2 cs( 1) / -0.0348409253 897013234 e0/ data f2 cs( 2) / -0.0166842205 677959686 e0/ data f2 cs( 3) / 0.0006752901 241237738 e0/ data f2 cs( 4) / -0.0000535066 622544701 e0/ data f2 cs( 5) / 0.0000062693 421779007 e0/ data f2 cs( 6) / -0.0000009526 638801991 e0/ data f2 cs( 7) / 0.0000001745 629224251 e0/ data f2 cs( 8) / -0.0000000368 795403065 e0/ data f2 cs( 9) / 0.0000000087 202677705 e0/ data f2 cs( 10) / -0.0000000022 601970392 e0/ data f2 cs( 11) / 0.0000000006 324624977 e0/ data f2 cs( 12) / -0.0000000001 888911889 e0/ data f2 cs( 13) / 0.0000000000 596774674 e0/ data f2 cs( 14) / -0.0000000000 198044313 e0/ data f2 cs( 15) / 0.0000000000 068641396 e0/ data f2 cs( 16) / -0.0000000000 024731020 e0/ data f2 cs( 17) / 0.0000000000 009226360 e0/ data f2 cs( 18) / -0.0000000000 003552364 e0/ data f2 cs( 19) / 0.0000000000 001407606 e0/ data f2 cs( 20) / -0.0000000000 000572623 e0/ data f2 cs( 21) / 0.0000000000 000238654 e0/ data f2 cs( 22) / -0.0000000000 000101714 e0/ data f2 cs( 23) / 0.0000000000 000044259 e0/ data f2 cs( 24) / -0.0000000000 000019634 e0/ data f2 cs( 25) / 0.0000000000 000008868 e0/ data f2 cs( 26) / -0.0000000000 000004074 e0/ data f2 cs( 27) / 0.0000000000 000001901 e0/ data f2 cs( 28) / -0.0000000000 000000900 e0/ data f2 cs( 29) / 0.0000000000 000000432 e0/ c c series for g1 on the interval 2.00000e-02 to 6.25000e-02 c with weighted error 5.48e-17 c log weighted error 16.26 c significant figures required 15.47 c decimal places required 16.92 c data g1 cs( 1) / -0.3040578798 253495954 e0/ data g1 cs( 2) / -0.0566890984 597120588 e0/ data g1 cs( 3) / 0.0039046158 173275644 e0/ data g1 cs( 4) / -0.0003746075 959202261 e0/ data g1 cs( 5) / 0.0000435431 556559844 e0/ data g1 cs( 6) / -0.0000057417 294453025 e0/ data g1 cs( 7) / 0.0000008282 552104503 e0/ data g1 cs( 8) / -0.0000001278 245892595 e0/ data g1 cs( 9) / 0.0000000207 978352949 e0/ data g1 cs( 10) / -0.0000000035 313205922 e0/ data g1 cs( 11) / 0.0000000006 210824236 e0/ data g1 cs( 12) / -0.0000000001 125215474 e0/ data g1 cs( 13) / 0.0000000000 209088918 e0/ data g1 cs( 14) / -0.0000000000 039715832 e0/ data g1 cs( 15) / 0.0000000000 007690431 e0/ data g1 cs( 16) / -0.0000000000 001514697 e0/ data g1 cs( 17) / 0.0000000000 000302892 e0/ data g1 cs( 18) / -0.0000000000 000061400 e0/ data g1 cs( 19) / 0.0000000000 000012601 e0/ data g1 cs( 20) / -0.0000000000 000002615 e0/ data g1 cs( 21) / 0.0000000000 000000548 e0/ c c series for g2 on the interval 0.00000e+00 to 2.00000e-02 c with weighted error 5.01e-17 c log weighted error 16.30 c significant figures required 15.12 c decimal places required 17.07 c data g2 cs( 1) / -0.0967329367 532432218 e0/ data g2 cs( 2) / -0.0452077907 957459871 e0/ data g2 cs( 3) / 0.0028190005 352706523 e0/ data g2 cs( 4) / -0.0002899167 740759160 e0/ data g2 cs( 5) / 0.0000407444 664601121 e0/ data g2 cs( 6) / -0.0000071056 382192354 e0/ data g2 cs( 7) / 0.0000014534 723163019 e0/ data g2 cs( 8) / -0.0000003364 116512503 e0/ data g2 cs( 9) / 0.0000000859 774367886 e0/ data g2 cs( 10) / -0.0000000238 437656302 e0/ data g2 cs( 11) / 0.0000000070 831906340 e0/ data g2 cs( 12) / -0.0000000022 318068154 e0/ data g2 cs( 13) / 0.0000000007 401087359 e0/ data g2 cs( 14) / -0.0000000002 567171162 e0/ data g2 cs( 15) / 0.0000000000 926707021 e0/ data g2 cs( 16) / -0.0000000000 346693311 e0/ data g2 cs( 17) / 0.0000000000 133950573 e0/ data g2 cs( 18) / -0.0000000000 053290754 e0/ data g2 cs( 19) / 0.0000000000 021775312 e0/ data g2 cs( 20) / -0.0000000000 009118621 e0/ data g2 cs( 21) / 0.0000000000 003905864 e0/ data g2 cs( 22) / -0.0000000000 001708459 e0/ data g2 cs( 23) / 0.0000000000 000762015 e0/ data g2 cs( 24) / -0.0000000000 000346151 e0/ data g2 cs( 25) / 0.0000000000 000159996 e0/ data g2 cs( 26) / -0.0000000000 000075213 e0/ data g2 cs( 27) / 0.0000000000 000035970 e0/ data g2 cs( 28) / -0.0000000000 000017530 e0/ data g2 cs( 29) / 0.0000000000 000008738 e0/ data g2 cs( 30) / -0.0000000000 000004487 e0/ data g2 cs( 31) / 0.0000000000 000002397 e0/ data g2 cs( 32) / -0.0000000000 000001347 e0/ data g2 cs( 33) / 0.0000000000 000000801 e0/ data g2 cs( 34) / -0.0000000000 000000501 e0/ c data nf1, nf2, ng1, ng2 / 4*0 / data xbnd, xbig, xmaxf, xmaxg / 4*0.0 / c if (nf1.ne.0) go to 10 tol = 0.1*r1mach(3) nf1 = inits (f1cs, 20, tol) nf2 = inits (f2cs, 29, tol) ng1 = inits (g1cs, 21, tol) ng2 = inits (g2cs, 34, tol) c xbig = sqrt (1.0/r1mach(3)) xmaxf = exp (amin1(-alog(r1mach(1)), alog(r1mach(2))) - 0.01) xmaxg = 1.0/sqrt(r1mach(1)) xbnd = sqrt(50.0) c 10 if (x.lt.4.0) call seteru ( 1 34hr9sifg approxs invalid for x lt 4, 34, 1, 2) c if (x.gt.xbnd) go to 20 f = (1.0 + csevl ((1.0/x**2-0.04125)/0.02125, f1cs, nf1))/x g = (1.0 + csevl ((1.0/x**2-0.04125)/0.02125, g1cs, ng1))/x**2 return c 20 if (x.gt.xbig) go to 30 f = (1.0 + csevl (100./x**2-1., f2cs, nf2))/x g = (1.0 + csevl (100./x**2-1., g2cs, ng2))/x**2 return c 30 f = 0.0 if (x.lt.xmaxf) f = 1.0/x g = 0.0 if (x.lt.xmaxg) g = 1.0/x**2 return c end