Blame view
fvn_fnlib/erfc.f
5.47 KB
38581db0c
git-svn-id: https...
|
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
|