Blame view
fvn_fnlib/derfc.f
11.8 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 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 |
double precision function derfc (x)
c july 1977 edition. w. fullerton, c3, los alamos scientific lab.
double precision x, erfcs(21), erfccs(59), erc2cs(49), sqeps,
1 sqrtpi, xmax, xsml, y, d1mach, dcsevl, dexp, dlog, dsqrt
external d1mach, dcsevl, 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
c series for erc2 on the interval 2.50000e-01 to 1.00000e+00
c with weighted error 2.67e-32
c log weighted error 31.57
c significant figures required 30.31
c decimal places required 32.42
c
data erc2cs( 1) / -.6960134660 2309501127 3915082619 7 d-1 /
data erc2cs( 2) / -.4110133936 2620893489 8221208466 6 d-1 /
data erc2cs( 3) / +.3914495866 6896268815 6114370524 4 d-2 /
data erc2cs( 4) / -.4906395650 5489791612 8093545077 4 d-3 /
data erc2cs( 5) / +.7157479001 3770363807 6089414182 5 d-4 /
data erc2cs( 6) / -.1153071634 1312328338 0823284791 2 d-4 /
data erc2cs( 7) / +.1994670590 2019976350 5231486770 9 d-5 /
data erc2cs( 8) / -.3642666471 5992228739 3611843071 1 d-6 /
data erc2cs( 9) / +.6944372610 0050125899 3127721463 3 d-7 /
data erc2cs( 10) / -.1371220902 1043660195 3460514121 0 d-7 /
data erc2cs( 11) / +.2788389661 0071371319 6386034808 7 d-8 /
data erc2cs( 12) / -.5814164724 3311615518 6479105031 6 d-9 /
data erc2cs( 13) / +.1238920491 7527531811 8016881795 0 d-9 /
data erc2cs( 14) / -.2690639145 3067434323 9042493788 9 d-10 /
data erc2cs( 15) / +.5942614350 8479109824 4470968384 0 d-11 /
data erc2cs( 16) / -.1332386735 7581195792 8775442057 0 d-11 /
data erc2cs( 17) / +.3028046806 1771320171 7369724330 4 d-12 /
data erc2cs( 18) / -.6966648814 9410325887 9586758895 4 d-13 /
data erc2cs( 19) / +.1620854541 0539229698 1289322762 8 d-13 /
data erc2cs( 20) / -.3809934465 2504919998 7691305772 9 d-14 /
data erc2cs( 21) / +.9040487815 9788311493 6897101297 5 d-15 /
data erc2cs( 22) / -.2164006195 0896073478 0981204700 3 d-15 /
data erc2cs( 23) / +.5222102233 9958549846 0798024417 2 d-16 /
data erc2cs( 24) / -.1269729602 3645553363 7241552778 0 d-16 /
data erc2cs( 25) / +.3109145504 2761975838 3622741295 1 d-17 /
data erc2cs( 26) / -.7663762920 3203855240 0956671481 1 d-18 /
data erc2cs( 27) / +.1900819251 3627452025 3692973329 0 d-18 /
data erc2cs( 28) / -.4742207279 0690395452 2565599996 5 d-19 /
data erc2cs( 29) / +.1189649200 0765283828 8068307845 1 d-19 /
data erc2cs( 30) / -.3000035590 3257802568 4527131306 6 d-20 /
data erc2cs( 31) / +.7602993453 0432461730 1938527709 8 d-21 /
data erc2cs( 32) / -.1935909447 6068728815 6981104913 0 d-21 /
data erc2cs( 33) / +.4951399124 7733378810 0004238677 3 d-22 /
data erc2cs( 34) / -.1271807481 3363718796 0862198988 8 d-22 /
data erc2cs( 35) / +.3280049600 4695130433 1584165205 3 d-23 /
data erc2cs( 36) / -.8492320176 8228965689 2479242239 9 d-24 /
data erc2cs( 37) / +.2206917892 8075602235 1987998719 9 d-24 /
data erc2cs( 38) / -.5755617245 6965284983 1281950719 9 d-25 /
data erc2cs( 39) / +.1506191533 6392342503 5414405119 9 d-25 /
data erc2cs( 40) / -.3954502959 0187969531 0428569599 9 d-26 /
data erc2cs( 41) / +.1041529704 1515009799 8464505173 3 d-26 /
data erc2cs( 42) / -.2751487795 2787650794 5017890133 3 d-27 /
data erc2cs( 43) / +.7290058205 4975574089 9770368000 0 d-28 /
data erc2cs( 44) / -.1936939645 9159478040 7750109866 6 d-28 /
data erc2cs( 45) / +.5160357112 0514872983 7005482666 6 d-29 /
data erc2cs( 46) / -.1378419322 1930940993 8964480000 0 d-29 /
data erc2cs( 47) / +.3691326793 1070690422 5109333333 3 d-30 /
data erc2cs( 48) / -.9909389590 6243654206 5322666666 6 d-31 /
data erc2cs( 49) / +.2666491705 1953884133 2394666666 6 d-31 /
c
c series for erfc on the interval 0. to 2.50000e-01
c with weighted error 1.53e-31
c log weighted error 30.82
c significant figures required 29.47
c decimal places required 31.70
c
data erfccs( 1) / +.7151793102 0292477450 3697709496 d-1 /
data erfccs( 2) / -.2653243433 7606715755 8893386681 d-1 /
data erfccs( 3) / +.1711153977 9208558833 2699194606 d-2 /
data erfccs( 4) / -.1637516634 5851788416 3746404749 d-3 /
data erfccs( 5) / +.1987129350 0552036499 5974806758 d-4 /
data erfccs( 6) / -.2843712412 7665550875 0175183152 d-5 /
data erfccs( 7) / +.4606161308 9631303696 9379968464 d-6 /
data erfccs( 8) / -.8227753025 8792084205 7766536366 d-7 /
data erfccs( 9) / +.1592141872 7709011298 9358340826 d-7 /
data erfccs( 10) / -.3295071362 2528432148 6631665072 d-8 /
data erfccs( 11) / +.7223439760 4005554658 1261153890 d-9 /
data erfccs( 12) / -.1664855813 3987295934 4695966886 d-9 /
data erfccs( 13) / +.4010392588 2376648207 7671768814 d-10 /
data erfccs( 14) / -.1004816214 4257311327 2170176283 d-10 /
data erfccs( 15) / +.2608275913 3003338085 9341009439 d-11 /
data erfccs( 16) / -.6991110560 4040248655 7697812476 d-12 /
data erfccs( 17) / +.1929492333 2617070862 4205749803 d-12 /
data erfccs( 18) / -.5470131188 7543310649 0125085271 d-13 /
data erfccs( 19) / +.1589663309 7626974483 9084032762 d-13 /
data erfccs( 20) / -.4726893980 1975548392 0369584290 d-14 /
data erfccs( 21) / +.1435873376 7849847867 2873997840 d-14 /
data erfccs( 22) / -.4449510561 8173583941 7250062829 d-15 /
data erfccs( 23) / +.1404810884 7682334373 7305537466 d-15 /
data erfccs( 24) / -.4513818387 7642108962 5963281623 d-16 /
data erfccs( 25) / +.1474521541 0451330778 7018713262 d-16 /
data erfccs( 26) / -.4892621406 9457761543 6841552532 d-17 /
data erfccs( 27) / +.1647612141 4106467389 5301522827 d-17 /
data erfccs( 28) / -.5626817176 3294080929 9928521323 d-18 /
data erfccs( 29) / +.1947443382 2320785142 9197867821 d-18 /
data erfccs( 30) / -.6826305642 9484207295 6664144723 d-19 /
data erfccs( 31) / +.2421988887 2986492401 8301125438 d-19 /
data erfccs( 32) / -.8693414133 5030704256 3800861857 d-20 /
data erfccs( 33) / +.3155180346 2280855712 2363401262 d-20 /
data erfccs( 34) / -.1157372324 0496087426 1239486742 d-20 /
data erfccs( 35) / +.4288947161 6056539462 3737097442 d-21 /
data erfccs( 36) / -.1605030742 0576168500 5737770964 d-21 /
data erfccs( 37) / +.6063298757 4538026449 5069923027 d-22 /
data erfccs( 38) / -.2311404251 6979584909 8840801367 d-22 /
data erfccs( 39) / +.8888778540 6618855255 4702955697 d-23 /
data erfccs( 40) / -.3447260576 6513765223 0718495566 d-23 /
data erfccs( 41) / +.1347865460 2069650682 7582774181 d-23 /
data erfccs( 42) / -.5311794071 1250217364 5873201807 d-24 /
data erfccs( 43) / +.2109341058 6197831682 8954734537 d-24 /
data erfccs( 44) / -.8438365587 9237891159 8133256738 d-25 /
data erfccs( 45) / +.3399982524 9452089062 7359576337 d-25 /
data erfccs( 46) / -.1379452388 0732420900 2238377110 d-25 /
data erfccs( 47) / +.5634490311 8332526151 3392634811 d-26 /
data erfccs( 48) / -.2316490434 4770654482 3427752700 d-26 /
data erfccs( 49) / +.9584462844 6018101526 3158381226 d-27 /
data erfccs( 50) / -.3990722880 3301097262 4224850193 d-27 /
data erfccs( 51) / +.1672129225 9444773601 7228709669 d-27 /
data erfccs( 52) / -.7045991522 7660138563 8803782587 d-28 /
data erfccs( 53) / +.2979768402 8642063541 2357989444 d-28 /
data erfccs( 54) / -.1262522466 4606192972 2422632994 d-28 /
data erfccs( 55) / +.5395438704 5424879398 5299653154 d-29 /
data erfccs( 56) / -.2380992882 5314591867 5346190062 d-29 /
data erfccs( 57) / +.1099052830 1027615735 9726683750 d-29 /
data erfccs( 58) / -.4867713741 6449657273 2518677435 d-30 /
data erfccs( 59) / +.1525877264 1103575676 3200828211 d-30 /
c
data sqrtpi / 1.772453850 9055160272 9816748334 115d0 /
data nterf, nterfc, nterc2, xsml, xmax, sqeps / 3*0, 3*0.d0 /
c
if (nterf.ne.0) go to 10
eta = 0.1*sngl(d1mach(3))
nterf = initds (erfcs, 21, eta)
nterfc = initds (erfccs, 59, eta)
nterc2 = initds (erc2cs, 49, eta)
c
xsml = -dsqrt (-dlog(sqrtpi*d1mach(3)))
xmax = dsqrt (-dlog(sqrtpi*d1mach(1)) )
xmax = xmax - 0.5d0*dlog(xmax)/xmax - 0.01d0
sqeps = dsqrt (2.0d0*d1mach(3))
c
10 if (x.gt.xsml) go to 20
c
c erfc(x) = 1.0 - erf(x) for x .lt. xsml
c
derfc = 2.0d0
return
c
20 if (x.gt.xmax) go to 40
y = dabs(x)
if (y.gt.1.0d0) go to 30
c
c erfc(x) = 1.0 - erf(x) for abs(x) .le. 1.0
c
if (y.lt.sqeps) derfc = 1.0d0 - 2.0d0*x/sqrtpi
if (y.ge.sqeps) derfc = 1.0d0 - x*(1.0d0 + dcsevl (2.d0*x*x-1.d0,
1 erfcs, nterf))
return
c
c erfc(x) = 1.0 - erf(x) for 1.0 .lt. abs(x) .le. xmax
c
30 y = y*y
if (y.le.4.d0) derfc = dexp(-y)/dabs(x) * (0.5d0 + dcsevl (
1 (8.d0/y-5.d0)/3.d0, erc2cs, nterc2) )
if (y.gt.4.d0) derfc = dexp(-y)/dabs(x) * (0.5d0 + dcsevl (
1 8.d0/y-1.d0, erfccs, nterfc) )
if (x.lt.0.d0) derfc = 2.0d0 - derfc
return
c
40 call seteru (32hderfc x so big erfc underflows, 32, 1, 0)
derfc = 0.d0
return
c
end
|