function r9ln2r (x) c april 1978 edition. w. fullerton c3, los alamos scientific lab. c c evaluate alog(1+x) from 2-nd order with relative error accuracy so c that alog(1+x) = x - x**2/2 + x**3*r9ln2r(x) c real ln21cs(26), ln22cs(20) external csevl, inits, r1mach c c series for ln21 on the interval -6.25000d-01 to 0. c with weighted error 2.49e-17 c log weighted error 16.60 c significant figures required 15.87 c decimal places required 17.31 c data ln21cs( 1) / .1811196251 3478810e0 / data ln21cs( 2) / -.1562712319 2872463e0 / data ln21cs( 3) / .0286763053 61557275e0 / data ln21cs( 4) / -.0055586996 55948139e0 / data ln21cs( 5) / .0011178976 65229983e0 / data ln21cs( 6) / -.0002308050 89823279e0 / data ln21cs( 7) / .0000485988 53341100e0 / data ln21cs( 8) / -.0000103901 27388903e0 / data ln21cs( 9) / .0000022484 56370739e0 / data ln21cs(10) / -.0000004914 05927392e0 / data ln21cs(11) / .0000001082 82565070e0 / data ln21cs(12) / -.0000000240 25872763e0 / data ln21cs(13) / .0000000053 62460047e0 / data ln21cs(14) / -.0000000012 02995136e0 / data ln21cs(15) / .0000000002 71078892e0 / data ln21cs(16) / -.0000000000 61323562e0 / data ln21cs(17) / .0000000000 13920858e0 / data ln21cs(18) / -.0000000000 03169930e0 / data ln21cs(19) / .0000000000 00723837e0 / data ln21cs(20) / -.0000000000 00165700e0 / data ln21cs(21) / .0000000000 00038018e0 / data ln21cs(22) / -.0000000000 00008741e0 / data ln21cs(23) / .0000000000 00002013e0 / data ln21cs(24) / -.0000000000 00000464e0 / data ln21cs(25) / .0000000000 00000107e0 / data ln21cs(26) / -.0000000000 00000024e0 / c c series for ln22 on the interval 0. to 8.12500d-01 c with weighted error 1.42e-17 c log weighted error 16.85 c significant figures required 15.95 c decimal places required 17.50 c data ln22cs( 1) / -.2224253253 5020461e0 / data ln22cs( 2) / -.0610471001 08078624e0 / data ln22cs( 3) / .0074272350 09750394e0 / data ln22cs( 4) / -.0009335018 26163697e0 / data ln22cs( 5) / .0001200499 07687260e0 / data ln22cs( 6) / -.0000157047 22952820e0 / data ln22cs( 7) / .0000020818 74781051e0 / data ln22cs( 8) / -.0000002789 19557764e0 / data ln22cs( 9) / .0000000376 93558237e0 / data ln22cs(10) / -.0000000051 30902896e0 / data ln22cs(11) / .0000000007 02714117e0 / data ln22cs(12) / -.0000000000 96748595e0 / data ln22cs(13) / .0000000000 13381046e0 / data ln22cs(14) / -.0000000000 01858102e0 / data ln22cs(15) / .0000000000 00258929e0 / data ln22cs(16) / -.0000000000 00036195e0 / data ln22cs(17) / .0000000000 00005074e0 / data ln22cs(18) / -.0000000000 00000713e0 / data ln22cs(19) / .0000000000 00000100e0 / data ln22cs(20) / -.0000000000 00000014e0 / c data ntln21, ntln22, xmin, xbig, xmax / 2*0, 3*0.0 / c if (ntln21.ne.0) go to 10 eps = r1mach(3) ntln21 = inits (ln21cs, 26, 0.1*eps) ntln22 = inits (ln22cs, 20, 0.1*eps) c xmin = -1.0 + sqrt(r1mach(4)) sqeps = sqrt(eps) xmax = 6.0/sqeps xmax = xmax - (eps*xmax**2 - 2.0*alog(xmax)) / (2.0*eps*xmax) xbig = 4.0/sqrt(sqeps) xbig = xbig - (sqeps*xbig**2 - 2.0*alog(xbig)) / (2.*sqeps*xbig) c 10 if (x.lt.(-0.625) .or. x.gt.0.8125) go to 20 c if (x.lt.0.0) r9ln2r = 0.375 + csevl (16.*x/5.+1.0, ln21cs, 1 ntln21) if (x.ge.0.0) r9ln2r = 0.375 + csevl (32.*x/13.-1.0, ln22cs, 1 ntln22) return c 20 if (x.lt.xmin) call seteru ( 1 57hr9ln2r answer lt half precision because x is too near -1, 2 57, 1, 1) if (x.gt.xmax) call seteru ( 1 51hr9ln2r no precision in answer because x is too big, 51,3,2) if (x.gt.xbig) call seteru ( 1 53hr9ln2r answer lt half precision because x is too big, 53, 2 2, 1) c r9ln2r = (alog(1.0+x) - x*(1.0-0.5*x) ) / x**3 return c end