1
2
3
4
5
6
7
8
9
10 derf_: procedure (number) returns (float binary (63));
11
12
13 declare (number, r) float binary (63),
14 (f, n, p, q) float binary (63),
15 (abs, exp) builtin,
16 dexerfc_ entry (float binary (63)) returns (float binary (63));
17 r = 0.0e0;
18 erfs: n = number + 0.0e0;
19 f = abs (n);
20 if f < 0.47693628e0 then small: do;
21 q = 1.128379167095512574e0 * n;
22 if f < 1.e-10 then go to comp;
23 p = f*f;
24 n = (((( 0.8350702795147239592e-10 *p - 0.1089222103714857338e-8)*p + 0.1312253296380280507e-7)*p
25 - 0.1450385222315046876e-6)*p + 0.1458916900093370682e-5)*p;
26 q = ((((((((n - 0.1322751322751322751e-4)*p + 0.1068376068376068376e-3)*p - 0.7575757575757575758e-3)*p
27 + 0.4629629629629629630e-2)*p - 0.2380952380952380952e-1)*p + 0.1e0)*p
28 - 0.3333333333333333333e0)*p + 1.e0)*q;
29 comp: if r ^= 0.0e0 then q = 1.e0 - q;
30 go to finis;
31 end small;
32 if f >= 2.5e0 then large: do;
33 if f >= 9.30630096e0 then q = 0.0e0;
34 else q = exp (-f*f) * dexerfc_ (f);
35 end large; else middle: do;
36 if f < 1.5e0 then lower: do;
37 p = f - 1.e0;
38 q = ((((( - 0.8445354974538876839e-12 *p + 0.2035190183702655211e-11)*p + 0.1016478011836914822e-10)*p
39 - 0.3914544150228919332e-10)*p - 0.9687156253803097249e-10)*p + 0.6116760848521522320e-9)*p;
40 q = ((((((( q + 0.5758182413135021396e-9)*p - 0.7972478608404360084e-8)*p + 0.1720401851653628375e-8)*p
41 + 0.8630591951220226579e-7)*p - 0.1092604901414462467e-6)*p - 0.7524484361200326476e-6)*p
42 + 0.1844279900355114423e-5)*p;
43 q = (((((((((q + 0.4854967260442840621e-5)*p - 0.2100986406780402429e-4)*p - 0.1658718964594168655e-4)*p
44 + 0.1772942579639506779e-3)*p - 0.7579366392581423493e-4)*p - 0.1086769072243678816e-2)*p
45 + 0.1670704693150730120e-2)*p + 0.4233533251576121955e-2)*p - 0.1343051928086217999e-1)*p;
46 q = ((((((( q - 0.4087549346349359129e-2)*p + 0.6131324019524038693e-1)*p - 0.6131324019524038693e-1)*p
47 - 0.1226264803904807739e0)*p + 0.3678794411714423216e0)*p - 0.3678794411714423216e0)*p
48 + 0.1394027926403309882e0) * 1.128379167095512574e0;
49 end lower; else upper: do;
50 p = f - 2.e0;
51 q = (( - 0.1856500126366457284e-12 *p - 0.1944964938417663150e-12)*p + 0.3103020133929132565e-11)*p;
52 q = (((((((((q - 0.3723550879099818846e-11)*p - 0.3426266649193350739e-10)*p + 0.1200727503356961138e-9)*p
53 + 0.1787795198978365218e-9)*p - 0.1821176414744498095e-8)*p + 0.1759567016164738904e-8)*p
54 + 0.1642444033947200479e-7)*p - 0.5324702588728536810e-7)*p - 0.5245213918278798165e-7)*p;
55 q = (((((((((q + 0.6206354808105919169e-6)*p - 0.8484562971386512658e-6)*p - 0.3501612055936534974e-5)*p
56 + 0.1439484990506010484e-4)*p - 0.4634950036778170257e-5)*p - 0.9195995379200109923e-4)*p
57 + 0.2329025685851383420e-3)*p + 0.4441623187303262417e-4)*p - 0.1410013470005726578e-2)*p;
58 q = ((((((((q + 0.2994461596094635826e-2)*p - 0.4070141975274262287e-3)*p - 0.1159990462953164752e-1)*p
59 + 0.3052606481455696716e-1)*p - 0.4273649074037975402e-1)*p + 0.3663127777746836059e-1)*p
60 - 0.1831563888873418029e-1)*p + 0.4145534690336333682e-2) * 1.128379167095512574e0;
61 end upper; end middle;
62 if r = 0.0e0 then q = 1.e0 - q;
63 if n < 0.0e0 then q = r - q;
64 finis: return (q);
65 derfc_: entry (number) returns (float binary (63));
66 r = 2.e0;
67 go to erfs;
68 end derf_;