|
18 | 18 | *
|
19 | 19 | * ## Notice
|
20 | 20 | *
|
21 |
| -* The original C code, long comment, copyright, license, and constants are from [Cephes]{@link http://www.netlib.org/cephes}. The implementation follows the original, but has been modified for JavaScript. |
| 21 | +* The original C code, long comment, copyright, license, and constants are from [Cephes]{@link http://www.netlib.org/cephes}. The implementation follows the original, but has been modified according to project conventions. |
22 | 22 | *
|
23 | 23 | * ```text
|
24 | 24 | * Copyright 1984, 1995, 2000 by Stephen L. Moshier
|
|
38 | 38 | #include "stdlib/constants/float64/ln_two.h"
|
39 | 39 | #include <stdint.h>
|
40 | 40 |
|
| 41 | +// MAXLOG + LN2 = ln(2^1024) + LN2 |
41 | 42 | static const double POS_OVERFLOW = 7.09782712893383996843e2 + STDLIB_CONSTANT_FLOAT64_LN2;
|
| 43 | + |
| 44 | +// MINLOG - LN2 = ln(2^-1022) - LN2 |
42 | 45 | static const double NEG_OVERFLOW = -7.08396418532264106224e2 - STDLIB_CONSTANT_FLOAT64_LN2;
|
| 46 | + |
| 47 | +// MAXLOG - LN2 = ln(2^1024) - LN2 |
43 | 48 | static const double LARGE = 7.09782712893383996843e2 - STDLIB_CONSTANT_FLOAT64_LN2;
|
44 | 49 |
|
45 | 50 | /* Begin auto-generated functions. The following functions are auto-generated. Do not edit directly. */
|
@@ -88,40 +93,67 @@ static double rational_pq( const double x ) {
|
88 | 93 | /* End auto-generated functions. */
|
89 | 94 |
|
90 | 95 | /**
|
91 |
| -* Computes the hyperbolic sine of a number. |
| 96 | +* Computes the hyperbolic sine of a double-precision floating-point number. |
| 97 | +* |
| 98 | +* ## Method |
| 99 | +* |
| 100 | +* The range is partitioned into two segments. If \\( |x| \le 1 \\), we use a rational function of the form |
| 101 | +* |
| 102 | +* ```tex |
| 103 | +* x + x^3 \frac{\mathrm{P}(x)}{\mathrm{Q}(x)} |
| 104 | +* ``` |
| 105 | +* |
| 106 | +* Otherwise, the calculation is |
| 107 | +* |
| 108 | +* ```tex |
| 109 | +* \operatorname{sinh}(x) = \frac{ e^x - e^{-x} }{2}. |
| 110 | +* ``` |
| 111 | +* |
| 112 | +* ## Notes |
| 113 | +* |
| 114 | +* - Relative error: |
| 115 | +* |
| 116 | +* | arithmetic | domain | # trials | peak | rms | |
| 117 | +* |:----------:|:--------:|:--------:|:-------:|:-------:| |
| 118 | +* | DEC | +- 88 | 50000 | 4.0e-17 | 7.7e-18 | |
| 119 | +* | IEEE | +-MAXLOG | 30000 | 2.6e-16 | 5.7e-17 | |
92 | 120 | *
|
93 | 121 | * @param x input value
|
94 | 122 | * @return output value
|
95 | 123 | *
|
96 | 124 | * @example
|
97 | 125 | * double out = stdlib_base_sinh( 0.0 );
|
98 | 126 | * // returns 0.0
|
| 127 | +* |
| 128 | +* @example |
| 129 | +* double out = stdlib_base_sinh( 2.0 ); |
| 130 | +* // returns ~3.627 |
99 | 131 | */
|
100 | 132 | double stdlib_base_sinh( const double x ) {
|
101 |
| - double a; |
102 |
| - if ( x == 0.0 ) { |
103 |
| - return x; // handles `+-0` |
104 |
| - } |
105 |
| - a = stdlib_base_abs( x ); |
106 |
| - if ( x > POS_OVERFLOW || x < NEG_OVERFLOW ) { |
107 |
| - return ( x > 0.0 ) ? STDLIB_CONSTANT_FLOAT64_PINF : STDLIB_CONSTANT_FLOAT64_NINF; |
108 |
| - } |
109 |
| - if ( a > 1.0 ) { |
110 |
| - if ( a >= LARGE ) { |
111 |
| - a = stdlib_base_exp( 0.5 * a ); |
112 |
| - a *= 0.5 * a; |
113 |
| - if ( x < 0.0 ) { |
114 |
| - a = -a; |
115 |
| - } |
116 |
| - return a; |
117 |
| - } |
118 |
| - a = stdlib_base_exp( a ); |
119 |
| - a = ( 0.5 * a ) - ( 0.5 / a ); |
120 |
| - if ( x < 0.0 ) { |
121 |
| - a = -a; |
122 |
| - } |
123 |
| - return a; |
124 |
| - } |
125 |
| - a *= a; |
126 |
| - return x + ( x * a * rational_pq( a ) ); |
| 133 | + double a; |
| 134 | + if ( x == 0.0 ) { |
| 135 | + return x; // handles `+-0` |
| 136 | + } |
| 137 | + if ( x > POS_OVERFLOW || x < NEG_OVERFLOW ) { |
| 138 | + return ( x > 0.0 ) ? STDLIB_CONSTANT_FLOAT64_PINF : STDLIB_CONSTANT_FLOAT64_NINF; |
| 139 | + } |
| 140 | + a = stdlib_base_abs( x ); |
| 141 | + if ( a > 1.0 ) { |
| 142 | + if ( a >= LARGE ) { |
| 143 | + a = stdlib_base_exp( 0.5*a ); |
| 144 | + a *= 0.5 * a; |
| 145 | + if ( x < 0.0 ) { |
| 146 | + a = -a; |
| 147 | + } |
| 148 | + return a; |
| 149 | + } |
| 150 | + a = stdlib_base_exp( a ); |
| 151 | + a = ( 0.5*a ) - ( 0.5/a ); |
| 152 | + if ( x < 0.0 ) { |
| 153 | + a = -a; |
| 154 | + } |
| 155 | + return a; |
| 156 | + } |
| 157 | + a *= a; |
| 158 | + return x + ( x * a * rational_pq( a ) ); |
127 | 159 | }
|
0 commit comments