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 dlog_: procedure (number) returns (float binary (63));
12 
13 /*     compute the logarithm or hyperbolic arctangent of a double-precision floating-point number     */
14 declare (number, a, f, n, p, r) float binary (63),
15          i fixed binary (17);
16 dcl  code_ ext entry (fixed bin),
17      1 word aligned based,
18      2 exponent fixed bin (7) unal;
19 declare   (abs, addr, log, log10, log2, sign) builtin;
20 
21           return(log(number));                    /* Natural log of value. */
22 
23 long:     i = addr (f) -> word.exponent;
24           addr (f) -> word.exponent = 0;
25 if f < 0.7071067811865475244e0 then lower: do;
26                a = 0.5946035575013605334e0;
27                n = 0.75e0;
28 end lower; else upper: do;
29                a = 0.840896415253714543e0;
30                n = 0.25e0;
31           end upper;
32           r = f + a;
33           f = (f - a) / r;
34           n = i - n;
35 short:    if abs (f) < 1.e-19 then go to finis;
36           a = f*f;
37           f = ((((((((0.1176470588235294118e0 * a + 0.1333333333333333333e0) * a + 0.1538461538461538462e0) * a
38           + 0.1818181818181818182e0)* a + 0.2222222222222222222e0) * a + 0.2857142857142857143e0) * a
39           + 0.4e0) * a + 0.6666666666666666667e0) * a + 2.e0) * f;
40 finis:    return ((0.6931471805599453094e0 * n + f) * p);
41 
42 dlog2_: entry (number) returns (float binary (63));
43           return(log2(number));                             /* Log (2) of value. */
44 
45 dlog10_: entry (number) returns (float binary (63));
46           return(log10(number));                            /* Log (10) of value. */
47 
48 dlone_: entry (number) returns (float binary (63));
49           return(log(number+1.0e0));                        /* Natural log of x+1. */
50 
51 datanh_: entry (number) returns (float binary (63));
52           p = 0.5e0;
53           f = number;
54           a = abs (f);
55           if a < 0.1e0
56           then do;
57                n = 0.0e0;
58                go to short;
59                end;
60           a = a - 1.e0;
61           if a >= 0.0e0 then err2: do;
62                call code_ (sign (a) + 44);
63                if a = 0.0e0 then f = 170141182.e30 * f; else f = 0.0e0;
64                return (f);
65           end err2;
66           f = (1.e0 + f) / (1.e0 - f);
67           go to long;
68      end dlog_;