1
2
3
4
5
6
7
8
9
10 exerfc_: procedure (number) returns (float binary (27));
11
12
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_;