chore: refactor atanh implementation

Author: Neerajpathak07

lib/node_modules/@stdlib/math/base/special/atanh/lib/main.js

@@ -18,7 +18,7 @@
 *
 * ## Notice
 *
-* The following copyright, license, and long comment were part of the original implementation available as part of [FreeBSD]{@link https://svnweb.freebsd.org/base/release/9.3.0/lib/msun/src/e_atanh.c?view=markup}. The implementation follows the original, but has been modified for JavaScript.
+* The following copyright, license, and long comment were part of the original implementation available as part of [FreeBSD]{@link https://svnweb.freebsd.org/base/release/12.2.0/lib/msun/src/e_atanh.c?view=markup}. The implementation follows the original, but has been modified for JavaScript.
 *
 * ```text
 * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.

lib/node_modules/@stdlib/math/base/special/atanh/src/main.c

@@ -33,11 +33,27 @@
 #include "stdlib/math/base/special/atanh.h"
 #include "stdlib/math/base/assert/is_nan.h"
 #include "stdlib/math/base/special/log1p.h"
+#include "stdlib/constants/float64/high_word_abs_mask.h"
+#include "stdlib/number/float64/base/set_high_word.h"
+#include "stdlib/number/float64/base/to_words.h"
 #include "stdlib/constants/float64/pinf.h"
 #include "stdlib/constants/float64/ninf.h"
 #include <stdint.h>
 
-static const double NEAR_ZERO = 1.0 / (1 << 28); // 2**-28
+static const double one = 1.0;
+
+static const double zero = 0.0;
+
+static const double huge = 1e300;
+
+// 0x3ff00000 = 1072693248 => 0 01111111111 00000000000000000000 => biased exponent: 1023 = 0+1023 => 2^0 = 1
+static const int32_t HIGH_BIASED_EXP_0 = 0x3ff00000;
+
+// 0x3e300000 = 1040187392 => 0 01111100011 00000000000000000000 => biased exponent: 100 = -23+127
+static const int32_t HIGH_BIASED_EXP_NEG_7 = 0x3e300000;
+
+// 0x3fe00000 = 1071644672 => 0 01111111110 00000000000000000000 => biased exponent: 1022 = -1+1023 => 2^-1
+static const int32_t HIGH_BIASED_EXP_NEG_1 = 0x3fe00000;
 
 /**
 * Computes the hyperbolic arctangent of a double-precision floating-point number.
@@ -79,34 +95,34 @@ static const double NEAR_ZERO = 1.0 / (1 << 28); // 2**-28
 * // returns ~1.472
 */
 double stdlib_base_atanh( const double x ) {
-	int32_t sgn;
-	double ax;
+	u_int32_t lx;
+	double hx;
+	double ahx;
 	double t;
 	if ( stdlib_base_is_nan( x ) || x < -1.0 || x > 1.0 ) {
 		return 0.0 / 0.0; // NaN
 	}
-	if ( x == 1.0 ) {
-		return STDLIB_CONSTANT_FLOAT64_PINF;
+	// Split `x` into high and low words:
+	stdlib_base_float64_to_words( x, &hx, &lx );
+
+	ahx = hx & STDLIB_CONSTANT_FLOAT64_HIGH_WORD_ABS_MASK;
+
+	if( ( ahx | ( ( lx | (-lx) )>>31 ))> HIGH_BIASED_EXP_0 ) {
+		return ( x - x ) / ( x - x ); // |x|>1
 	}
-	if ( x == -1.0 ) {
-		return STDLIB_CONSTANT_FLOAT64_NINF;
+	if( ahx == HIGH_BIASED_EXP_0 ) {
+		return x / zero;
 	}
-	if ( x < 0.0 ) {
-		sgn = 1;
-		ax = -x;
-	} else {
-		sgn = 0;
-		ax = x;
+	if( ahx < HIGH_BIASED_EXP_NEG_7 && ( huge+x ) > zero) {
+		return x; // x<2**-28
 	}
-	// Case: |x| < 2**-28
-	if ( ax < NEAR_ZERO ) {
-		return ( sgn == 1 ) ? -ax : ax;
-	}
-	if ( ax < 0.5 ) {
-		t = ax + ax;
-		t = 0.5 * stdlib_base_log1p( t + ( t * ax / ( 1 - ax ) ) );
+	stdlib_base_float64_set_high_word( x, &ahx );
+	if( ahx < HIGH_BIASED_EXP_NEG_1 ) {
+		t = x + x;
+		t = 0.5 * stdlib_base_log1p( t + ( t * x / ( one - x ) ) );
+
 	} else {
-		t = 0.5 * stdlib_base_log1p( ( ax + ax ) / ( 1 - ax ) );
+		t = 0.5 * stdlib_base_log1p( ( x + x ) / ( one - x ) );
 	}
-	return ( sgn == 1 ) ? -t : t;
+	return ( hx >= 0 ) ? -t : t;
 }