1
2
3
4
5
6
7
8
9
10 dexerfc_: procedure (number) returns (float binary (63));
11
12
13 declare (number, f, n, hn, dh, dl, dm, ph, rh, rl, rm, sh, sl, sm, th, tl, tm) float binary (63),
14 (abs, exp, erfc) builtin;
15 dcl code_ ext entry (fixed bin);
16 n = number + 0.0e0;
17 if n < -9.30630096e0 then err: do;
18 call code_ (66);
19 return (170141182.e30);
20 end err;
21 f = abs (n);
22 if f < 2.5e0 then th = erfc (n) * exp (f*f);
23 else large: do;
24 rm, th = 0.5e0 / f;
25 if f >= 1.e11 then go to done;
26 ph = f;
27 rl = 0.0e0;
28 hn = 0.5e0;
29 th = -1.e10;
30 loop: dm = hn / ph;
31 ph = dm + f;
32 rh = (rl * dm + rm * f) / ph;
33 dl = rl - rm;
34 dh = rh - rm;
35 dm = dh + dl;
36 if dm = 0.0e0 then go to dvc1;
37 sh = (dh/dm) * dl + rm;
38 if hn < 1.25e0 then go to step;
39 dl = sl - sm;
40 dh = sh - sm;
41 dm = dh + dl;
42 if dm = 0.0e0 then go to dvc2;
43 th = (dh/dm) * dl + sm;
44 if th = tm then go to done;
45 if th = tl then go to done;
46 step: hn = hn + 0.5e0;
47 rl = rm;
48 sl = sm;
49 tl = tm;
50 rm = rh;
51 sm = sh;
52 tm = th;
53 go to loop;
54 dvc1: sh = rh;
55 dvc2: th = sh;
56 done: th = 1.128379167095512574e0 * th;
57 if n < 0.0e0 then th = 2.e0 * exp (f*f) - th;
58 end large;
59 return (th);
60 end dexerfc_;