1 /* ******************************************************
 2    *                                                    *
 3    *                                                    *
 4    * Copyright (c) 1972 by Massachusetts Institute of   *
 5    * Technology and Honeywell Information Systems, Inc. *
 6    *                                                    *
 7    *                                                    *
 8    ****************************************************** */
 9 
10 /* modified by A. Downing to remove calls to round_ */
11 
12 log_: procedure (number) returns (float binary (27));
13 
14 /*    compute the logarithm or hyperbolic arctangent of a single-precision floating-point number     */
15 declare (number, r) float binary (27),
16         (a, f, n, p) float binary (63),
17          i fixed binary (17);
18 dcl  code_ ext entry (fixed bin),
19      (abs, addr, log, log10, log2, round, sign) builtin,
20      1 word aligned based,
21      2 exponent fixed bin (7) unal;
22           return(log(number));                              /* Natural log of value. */
23 
24 long:     i = addr (f) -> exponent;
25           addr (f) -> exponent = 0;
26 if f < 0.7071067811865475244e0 then lower: do;
27                a = 0.5946035575013605334e0;
28                n = 0.75e0;
29 end lower; else upper: do;
30                a = 0.840896415253714543e0;
31                n = 0.25e0;
32           end upper;
33           f = (f - a) / (f + a);
34           n = i - n;
35 short:    if abs (f) < 0.7450580597e-8 then go to finis;
36           a = f*f;
37           f = (((0.2857142857142857143e0 * a + 0.4e0) * a + 0.6666666666666666667e0) * a + 2.e0) * f;
38 finis:    f = (0.6931471805599453094e0 * n + f) * p;
39           r = round (f, 28);
40           return (r);
41 
42 log2_: entry (number) returns (float binary (27));
43           return(log2(number));                             /* log(2) of value. */
44 
45 log10_: entry (number) returns (float binary (27));
46           return(log10(number));                            /* log(10) of value. */
47 
48 lone_: entry (number) returns (float binary (27));
49           return(log(number+1.0e0));                        /* Natural log of x+1. */
50 
51 atanh_: entry (number) returns (float binary (27));
52           p = 0.5e0;
53           f = number;
54           a = abs (f);
55           if a < 0.1e0        then do;
56                n = 0.0e0;
57                go to short;
58                end;
59           a = a - 1.e0;
60           if a >= 0.0e0 then err2: do;
61                call code_ (sign (a) + 42);
62                if a = 0.0e0 then f = 170141182.e30 * f; else f = 0.0e0;
63                return (f);
64           end err2;
65           f = (1.e0 + f) / (1.e0 - f);
66           go to long;
67      end log_;