function erf (x) c april 1977 version. w. fullerton, c3, los alamos scientific lab. dimension erfcs(13) external csevl, erfc, inits, r1mach c c series for erf on the interval 0. to 1.00000d+00 c with weighted error 7.10e-18 c log weighted error 17.15 c significant figures required 16.31 c decimal places required 17.71 c data erf cs( 1) / -.0490461212 34691808e0 / data erf cs( 2) / -.1422612051 0371364e0 / data erf cs( 3) / .0100355821 87599796e0 / data erf cs( 4) / -.0005768764 69976748e0 / data erf cs( 5) / .0000274199 31252196e0 / data erf cs( 6) / -.0000011043 17550734e0 / data erf cs( 7) / .0000000384 88755420e0 / data erf cs( 8) / -.0000000011 80858253e0 / data erf cs( 9) / .0000000000 32334215e0 / data erf cs(10) / -.0000000000 00799101e0 / data erf cs(11) / .0000000000 00017990e0 / data erf cs(12) / -.0000000000 00000371e0 / data erf cs(13) / .0000000000 00000007e0 / c data sqrtpi /1.772453850 9055160e0/ data nterf, xbig, sqeps / 0, 0., 0./ c if (nterf.ne.0) go to 10 nterf = inits (erfcs, 13, 0.1*r1mach(3)) xbig = sqrt(-alog(sqrtpi*r1mach(3))) sqeps = sqrt(2.0*r1mach(3)) c 10 y = abs(x) if (y.gt.1.) go to 20 c c erf(x) = 1. - erfc(x) for -1. .le. x .le. 1. c if (y.le.sqeps) erf = 2.0*x/sqrtpi if (y.gt.sqeps) erf = x*(1.0 + csevl(2.*x**2-1., erfcs, nterf)) return c c erf(x) = 1. - erfc(x) for abs(x) .gt. 1. c 20 if (y.le.xbig) erf = sign (1.0-erfc(y), x) if (y.gt.xbig) erf = sign (1.0, x) c return end