erfc.f
5.47 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
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
function erfc (x)
c july 1977 edition. w. fullerton, c3, los alamos scientific lab.
dimension erfcs(13), erfccs(24), erc2cs(23)
external csevl, 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
c series for erc2 on the interval 2.50000d-01 to 1.00000d+00
c with weighted error 5.22e-17
c log weighted error 16.28
c approx significant figures required 15.0
c decimal places required 16.96
c
data erc2cs( 1) / -.0696013466 02309501e0 /
data erc2cs( 2) / -.0411013393 62620893e0 /
data erc2cs( 3) / .0039144958 66689626e0 /
data erc2cs( 4) / -.0004906395 65054897e0 /
data erc2cs( 5) / .0000715747 90013770e0 /
data erc2cs( 6) / -.0000115307 16341312e0 /
data erc2cs( 7) / .0000019946 70590201e0 /
data erc2cs( 8) / -.0000003642 66647159e0 /
data erc2cs( 9) / .0000000694 43726100e0 /
data erc2cs(10) / -.0000000137 12209021e0 /
data erc2cs(11) / .0000000027 88389661e0 /
data erc2cs(12) / -.0000000005 81416472e0 /
data erc2cs(13) / .0000000001 23892049e0 /
data erc2cs(14) / -.0000000000 26906391e0 /
data erc2cs(15) / .0000000000 05942614e0 /
data erc2cs(16) / -.0000000000 01332386e0 /
data erc2cs(17) / .0000000000 00302804e0 /
data erc2cs(18) / -.0000000000 00069666e0 /
data erc2cs(19) / .0000000000 00016208e0 /
data erc2cs(20) / -.0000000000 00003809e0 /
data erc2cs(21) / .0000000000 00000904e0 /
data erc2cs(22) / -.0000000000 00000216e0 /
data erc2cs(23) / .0000000000 00000052e0 /
c
c series for erfc on the interval 0. to 2.50000d-01
c with weighted error 4.81e-17
c log weighted error 16.32
c approx significant figures required 15.0
c decimal places required 17.01
c
data erfccs( 1) / 0.0715179310 202925e0 /
data erfccs( 2) / -.0265324343 37606719e0 /
data erfccs( 3) / .0017111539 77920853e0 /
data erfccs( 4) / -.0001637516 63458512e0 /
data erfccs( 5) / .0000198712 93500549e0 /
data erfccs( 6) / -.0000028437 12412769e0 /
data erfccs( 7) / .0000004606 16130901e0 /
data erfccs( 8) / -.0000000822 77530261e0 /
data erfccs( 9) / .0000000159 21418724e0 /
data erfccs(10) / -.0000000032 95071356e0 /
data erfccs(11) / .0000000007 22343973e0 /
data erfccs(12) / -.0000000001 66485584e0 /
data erfccs(13) / .0000000000 40103931e0 /
data erfccs(14) / -.0000000000 10048164e0 /
data erfccs(15) / .0000000000 02608272e0 /
data erfccs(16) / -.0000000000 00699105e0 /
data erfccs(17) / .0000000000 00192946e0 /
data erfccs(18) / -.0000000000 00054704e0 /
data erfccs(19) / .0000000000 00015901e0 /
data erfccs(20) / -.0000000000 00004729e0 /
data erfccs(21) / .0000000000 00001432e0 /
data erfccs(22) / -.0000000000 00000439e0 /
data erfccs(23) / .0000000000 00000138e0 /
data erfccs(24) / -.0000000000 00000048e0 /
c
data sqrtpi /1.772453850 9055160e0/
data nterf, nterfc, nterc2, xsml, xmax, sqeps /3*0, 3*0./
c
if (nterf.ne.0) go to 10
eta = 0.1*r1mach(3)
nterf = inits (erfcs, 13, eta)
nterfc = inits (erfccs, 24, eta)
nterc2 = inits (erc2cs, 23, eta)
c
xsml = -sqrt (-alog(sqrtpi*r1mach(3)))
xmax = sqrt (-alog(sqrtpi*r1mach(1)))
xmax = xmax - 0.5*alog(xmax)/xmax - 0.01
sqeps = sqrt (2.0*r1mach(3))
c
10 if (x.gt.xsml) go to 20
c
c erfc(x) = 1.0 - erf(x) for x .lt. xsml
c
erfc = 2.
return
c
20 if (x.gt.xmax) go to 40
y = abs(x)
if (y.gt.1.0) go to 30
c
c erfc(x) = 1.0 - erf(x) for -1. .le. x .le. 1.
c
if (y.lt.sqeps) erfc = 1.0 - 2.0*x/sqrtpi
if (y.ge.sqeps) erfc = 1.0 -
1 x*(1.0 + csevl (2.*x*x-1., erfcs, nterf) )
return
c
c erfc(x) = 1.0 - erf(x) for 1. .lt. abs(x) .le. xmax
c
30 y = y*y
if (y.le.4.) erfc = exp(-y)/abs(x) * (0.5 + csevl ((8./y-5.)/3.,
1 erc2cs, nterc2) )
if (y.gt.4.) erfc = exp(-y)/abs(x) * (0.5 + csevl (8./y-1.,
1 erfccs, nterfc) )
if (x.lt.0.) erfc = 2.0 - erfc
return
c
40 call seteru (32herfc x so big erfc underflows, 32, 1, 0)
erfc = 0.
return
c
end