1 /* ******************************************************
 2    *                                                    *
 3    *                                                    *
 4    * Copyright (c) 1972 by Massachusetts Institute of   *
 5    * Technology and Honeywell Information Systems, Inc. *
 6    *                                                    *
 7    *                                                    *
 8    ****************************************************** */
 9 
10 exerfc_: procedure (number) returns (float binary (27));
11 
12 /*       compute the special error function of a single-precision floating-point number       */
13 declare (number, r, uh, ul, um) float binary (27),
14           (abs, exp, erfc) builtin,
15         (f, n, hn, dh, dl, dm, ph, rh, rl, rm, sh, sl, sm, th, tl, tm) float binary (63);
16 dcl  code_ ext entry (fixed bin),
17      round builtin;
18           n, r = number + 0.0e0;
19 if n < -9.30630096e0 then err: do;
20                call code_ (41);
21                return (170141182.e30);
22           end err;
23           f = abs (n);
24           if f < 2.5e0 then th = erfc (r) * exp (f*f);
25 else large: do;
26                rm, th = 0.5e0 / f;
27                if f >= 1.e5 then go to done;
28                ph = f;
29                rl = 0.0e0;
30                hn = 0.5e0;
31                uh = -1.e10;
32 loop:          dm = hn / ph;
33                ph = dm + f;
34                rh = (rl * dm + rm * f) / ph;
35                dl = rl - rm;
36                dh = rh - rm;
37                dm = dh + dl;
38                if dm = 0.0e0 then go to dvc1;
39                sh = (dh/dm) * dl + rm;
40                if hn < 1.25e0 then go to step;
41                dl = sl - sm;
42                dh = sh - sm;
43                dm = dh + dl;
44                if dm = 0.0e0 then go to dvc2;
45                th = (dh/dm) * dl + sm;
46                uh = th;
47                if uh = um then go to done;
48                if uh = ul then go to done;
49 step:          hn = hn + 0.5e0;
50                rl = rm;
51                sl = sm;
52                tl = tm;
53                ul = um;
54                rm = rh;
55                sm = sh;
56                tm = th;
57                um = uh;
58                go to loop;
59 dvc1:          sh = rh;
60 dvc2:          th = sh;
61 done:          th = 1.128379167095512574e0 * th;
62                if n < 0.0e0 then th = 2.e0 * exp (f*f) - th;
63           end large;
64           r = round (th, 28);
65           return (r);
66      end exerfc_;