|
40 | 40 | #include "stdlib/constants/float64/ninf.h" |
41 | 41 | #include <stdint.h> |
42 | 42 |
|
43 | | -static const double one = 1.0; |
44 | | - |
45 | | -static const double zero = 0.0; |
46 | | - |
47 | | -static const double huge = 1e300; |
48 | | - |
49 | | -// 0x3ff00000 = 1072693248 => 0 01111111111 00000000000000000000 => biased exponent: 1023 = 0+1023 => 2^0 = 1 |
50 | | -static const int32_t HIGH_BIASED_EXP_0 = 0x3ff00000; |
51 | | - |
52 | | -// 0x3e300000 = 1040187392 => 0 01111100011 00000000000000000000 => biased exponent: 100 = -23+127 |
53 | | -static const int32_t HIGH_BIASED_EXP_NEG_7 = 0x3e300000; |
54 | | - |
55 | | -// 0x3fe00000 = 1071644672 => 0 01111111110 00000000000000000000 => biased exponent: 1022 = -1+1023 => 2^-1 |
56 | | -static const int32_t HIGH_BIASED_EXP_NEG_1 = 0x3fe00000; |
| 43 | +static const double NEAR_ZERO = 1.0 / (1 << 28); // 2**-28 |
57 | 44 |
|
58 | 45 | /** |
59 | 46 | * Computes the hyperbolic arctangent of a double-precision floating-point number. |
@@ -95,37 +82,42 @@ static const int32_t HIGH_BIASED_EXP_NEG_1 = 0x3fe00000; |
95 | 82 | * // returns ~1.472 |
96 | 83 | */ |
97 | 84 | double stdlib_base_atanh( const double x ) { |
| 85 | + uint32_t ahx; |
98 | 86 | uint32_t lx; |
99 | 87 | uint32_t hx; |
100 | | - uint32_t ahx; |
101 | 88 | double ax; |
102 | 89 | double t; |
| 90 | + |
103 | 91 | if ( stdlib_base_is_nan( x ) || x < -1.0 || x > 1.0 ) { |
104 | 92 | return 0.0 / 0.0; // NaN |
105 | 93 | } |
| 94 | + |
106 | 95 | // Split `x` into high and low words: |
107 | 96 | stdlib_base_float64_to_words( x, &hx, &lx ); |
108 | 97 |
|
| 98 | + ax = x; |
109 | 99 | ahx = hx & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK; |
110 | 100 |
|
111 | | - if( ( ahx | ( ( lx | (-lx) )>>31 ))> HIGH_BIASED_EXP_0 ) { |
112 | | - return ( x - x ) / ( x - x ); // |x|>1 |
| 101 | + // special cases |
| 102 | + if ( ax == 1.0 ) { |
| 103 | + return STDLIB_CONSTANT_FLOAT64_PINF; |
113 | 104 | } |
114 | | - if( ahx == HIGH_BIASED_EXP_0 ) { |
115 | | - return x / zero; |
| 105 | + if ( ax == -1.0 ) { |
| 106 | + return STDLIB_CONSTANT_FLOAT64_NINF; |
116 | 107 | } |
117 | | - if( ahx < HIGH_BIASED_EXP_NEG_7 && ( huge+x ) > zero) { |
118 | | - return x; // x<2**-28 |
| 108 | + |
| 109 | + // Case: |x| < 2**-28 |
| 110 | + if ( ax < NEAR_ZERO ) { |
| 111 | + return ( hx == 1 ) ? -ax : ax; |
119 | 112 | } |
120 | 113 |
|
121 | | - ax = x; |
122 | 114 | stdlib_base_float64_set_high_word( ahx, &ax ); |
123 | 115 |
|
124 | | - if( ahx < HIGH_BIASED_EXP_NEG_1 ) { |
| 116 | + if( ahx < 0.5 ) { |
125 | 117 | t = ax + ax; |
126 | | - t = 0.5 * stdlib_base_log1p( t + ( t * ax / ( one - ax ) ) ); |
| 118 | + t = 0.5 * stdlib_base_log1p( t + ( t * ax / ( 1.0 - ax ) ) ); |
127 | 119 | } else { |
128 | | - t = 0.5 * stdlib_base_log1p( ( ax + ax ) / ( one - ax ) ); |
| 120 | + t = 0.5 * stdlib_base_log1p( ( ax + ax ) / ( 1.0 - ax ) ); |
129 | 121 | } |
130 | 122 | return ( hx == 1 ) ? -t : t; |
131 | 123 | } |
0 commit comments