1 /* ******************************************************
 2    *                                                    *
 3    *                                                    *
 4    * Copyright (c) 1972 by Massachusetts Institute of   *
 5    * Technology and Honeywell Information Systems, Inc. *
 6    *                                                    *
 7    *                                                    *
 8    ****************************************************** */
 9 
10 dsinh_: procedure (number) returns (float binary (63));
11 
12 /*      compute the hyperbolic sine or cosine of a double-precision floating-point number      */
13 declare (number, b, n, p, r) float binary (63),
14           (abs, exp) builtin;
15 dcl  code_ ext entry (fixed bin);
16           n = number;
17           r = abs (n);
18 if r >= 0.7135236e0 then large: do;
19                p = -1.e0;
20                go to sinhs;
21           end large;
22           if r < 1.e-10 then go to negate;
23           b = r*r;
24           r = ((((((((0.2811457254345520763e-14 * b + 0.7647163731819816476e-12) * b + 0.1605904383682161460e-9) * b
25           + 0.2505210838544171878e-7) * b + 0.2755731922398589065e-5) * b + 0.1984126984126984127e-3) * b
26           + 0.8333333333333333333e-2) * b + 0.1666666666666666667e0) * b + 1.e0) * n;
27           go to finis;
28 dcosh_: entry (number) returns (float binary (63));
29           n, r = abs (number);
30           p = 1.e0;
31 sinhs: if r > 88.028e0 then err: do;
32                call code_ (50);
33                r = 170141182.e30;
34                go to negate;
35           end err;
36           r = exp (r);
37           r = (p/r + r) * 0.5e0;
38 negate:   if n < 0.0e0 then r = -r;
39 finis:    return (r);
40      end dsinh_;