1 /* ******************************************************
 2    *                                                    *
 3    *                                                    *
 4    * Copyright (c) 1972 by Massachusetts Institute of   *
 5    * Technology and Honeywell Information Systems, Inc. *
 6    *                                                    *
 7    *                                                    *
 8    ****************************************************** */
 9 
10 dexerfc_: procedure (number) returns (float binary (63));
11 
12 /*        compute the special error function of a double-precision floating-point number        */
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_;