derf.f
2.83 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
50
51
52
53
54
55
56
57
58
59
60
double precision function derf (x)
c july 1977 edition. w. fullerton, c3, los alamos scientific lab.
c last revised june 1983. w. fullerton, imsl, houston.
double precision x, erfcs(21), sqeps, sqrtpi, xbig, y, d1mach,
1 dcsevl, derfc, dlog, dsqrt
external d1mach, dcsevl, derfc, initds
c
c series for erf on the interval 0. to 1.00000e+00
c with weighted error 1.28e-32
c log weighted error 31.89
c significant figures required 31.05
c decimal places required 32.55
c
data erf cs( 1) / -.4904612123 4691808039 9845440333 76 d-1 /
data erf cs( 2) / -.1422612051 0371364237 8247418996 31 d+0 /
data erf cs( 3) / +.1003558218 7599795575 7546767129 33 d-1 /
data erf cs( 4) / -.5768764699 7674847650 8270255091 67 d-3 /
data erf cs( 5) / +.2741993125 2196061034 4221607914 71 d-4 /
data erf cs( 6) / -.1104317550 7344507604 1353812959 05 d-5 /
data erf cs( 7) / +.3848875542 0345036949 9613114981 74 d-7 /
data erf cs( 8) / -.1180858253 3875466969 6317518015 81 d-8 /
data erf cs( 9) / +.3233421582 6050909646 4029309533 54 d-10 /
data erf cs( 10) / -.7991015947 0045487581 6073747085 95 d-12 /
data erf cs( 11) / +.1799072511 3961455611 9672454866 34 d-13 /
data erf cs( 12) / -.3718635487 8186926382 3168282094 93 d-15 /
data erf cs( 13) / +.7103599003 7142529711 6899083946 66 d-17 /
data erf cs( 14) / -.1261245511 9155225832 4954248533 33 d-18 /
data erf cs( 15) / +.2091640694 1769294369 1705002666 66 d-20 /
data erf cs( 16) / -.3253973102 9314072982 3641600000 00 d-22 /
data erf cs( 17) / +.4766867209 7976748332 3733333333 33 d-24 /
data erf cs( 18) / -.6598012078 2851343155 1999999999 99 d-26 /
data erf cs( 19) / +.8655011469 9637626197 3333333333 33 d-28 /
data erf cs( 20) / -.1078892517 7498064213 3333333333 33 d-29 /
data erf cs( 21) / +.1281188399 3017002666 6666666666 66 d-31 /
c
data sqrtpi / 1.772453850 9055160272 9816748334 115d0 /
data nterf, xbig, sqeps / 0, 2*0.d0 /
c
if (nterf.ne.0) go to 10
nterf = initds (erfcs, 21, 0.1*sngl(d1mach(3)))
xbig = dsqrt (-dlog(sqrtpi*d1mach(3)))
sqeps = dsqrt (2.0d0*d1mach(3))
c
10 y = dabs(x)
if (y.gt.1.d0) go to 20
c
c erf(x) = 1.0 - erfc(x) for -1.0 .le. x .le. 1.0
c
if (y.le.sqeps) derf = 2.0d0*x/sqrtpi
if (y.gt.sqeps) derf = x*(1.0d0 + dcsevl (2.d0*x*x-1.d0,
1 erfcs, nterf))
return
c
c erf(x) = 1.0 - erfc(x) for abs(x) .gt. 1.0
c
20 if (y.le.xbig) derf = dsign (1.0d0-derfc(y), x)
if (y.gt.xbig) derf = dsign (1.0d0, x)
c
return
end