1 /* ******************************************************
 2    *                                                    *
 3    *                                                    *
 4    * Copyright (c) 1972 by Massachusetts Institute of   *
 5    * Technology and Honeywell Information Systems, Inc. *
 6    *                                                    *
 7    *                                                    *
 8    ****************************************************** */
 9 
10 /* calls to round_ removed 07/16/73 by A. Downing */
11 erf_: procedure (number) returns (float binary (27));
12 
13 /*    compute the error (or complementary error) function of a single-precision floating-point number    */
14 declare (number, s) float binary (27),
15         (f, n, p, q, r) float binary (63),
16          exerfc_ entry (float binary (27)) returns (float binary (27)),
17           (abs, exp, round) builtin;
18           r = 0.0e0;
19 erfs:     n = number;
20           f = abs (n);
21 if f < 0.5e0 then small: do;
22                q = 1.128379167095512574e0 * n;
23                if f < 5.e-5 then go to comp;
24                p = f*f;
25                q = ((((((( - 0.1322751322751322751e-4 *p + 0.1068376068376068376e-3)*p - 0.7575757575757575758e-3)*p
26                + 0.4629629629629629630e-2)*p - 0.2380952380952380952e-1)*p + 0.1e0)*p
27                - 0.3333333333333333333e0)*p + 1.e0)*q;
28 comp:          if r ^= 0.0e0 then q = 1.e0 - q;
29                go to finis;
30           end small;
31 if f >= 2.5e0 then large: do;
32                s = f;
33                if f >= 9.30630096e0 then q = 0.0e0;
34                else q = exp (-f*f) * exerfc_ (s);
35 end large; else middle: do;
36 if f < 1.5e0 then lower: do;
37                     p = f - 1.e0;
38                     q = (((((((( 0.4854967260442840621e-5 *p - 0.2100986406780402429e-4)*p - 0.1658718964594168655e-4)*p
39                     + 0.1772942579639506779e-3)*p - 0.7579366392581423493e-4)*p - 0.1086769072243678816e-2)*p
40                     + 0.1670704693150730120e-2)*p + 0.4233533251576121955e-2)*p - 0.1343051928086217999e-1)*p;
41                     q = ((((((( q - 0.4087549346349359129e-2)*p + 0.6131324019524038693e-1)*p - 0.6131324019524038693e-1)*p
42                     - 0.1226264803904807739e0)*p + 0.3678794411714423216e0)*p - 0.3678794411714423216e0)*p
43                     + 0.1394027926403309882e0) * 1.128379167095512574e0;
44 end lower; else upper: do;
45                     p = f - 2.e0;
46                     q = (((((((( 0.6206354808105919169e-6 *p - 0.8484562971386512658e-6)*p - 0.3501612055936534974e-5)*p
47                     + 0.1439484990506010484e-4)*p - 0.4634950036778170257e-5)*p - 0.9195995379200109923e-4)*p
48                     + 0.2329025685851383420e-3)*p + 0.4441623187303262417e-4)*p - 0.1410013470005726578e-2)*p;
49                     q = ((((((((q + 0.2994461596094635826e-2)*p - 0.4070141975274262287e-3)*p - 0.1159990462953164752e-1)*p
50                     + 0.3052606481455696716e-1)*p - 0.4273649074037975402e-1)*p + 0.3663127777746836059e-1)*p
51                     - 0.1831563888873418029e-1)*p + 0.4145534690336333682e-2) * 1.128379167095512574e0;
52           end upper; end middle;
53           if r = 0.0e0 then q = 1.e0 - q;
54           if n < 0.0e0 then q = r - q;
55 finis:              s = round (q, 28);
56           return (s);
57 erfc_: entry (number) returns (float binary (27));
58           r = 2.e0;
59           go to erfs;
60      end erf_;