r9sifg.f 8.16 KB
subroutine r9sifg (x, f, g)
c december 1980 edition.  w. fullerton, bell labs
      dimension f1cs(20), f2cs(29), g1cs(21), g2cs(34)
      external csevl, inits, r1mach
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