erf.f
1.77 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
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