1 /* 2 * Double-precision atanh(x) function. 3 * 4 * Copyright (c) 2022-2023, Arm Limited. 5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 6 */ 7 8 #include "math_config.h" 9 #include "estrin.h" 10 #include "pl_sig.h" 11 #include "pl_test.h" 12 13 #define AbsMask 0x7fffffffffffffff 14 #define Half 0x3fe0000000000000 15 #define One 0x3ff0000000000000 16 #define Ln2Hi 0x1.62e42fefa3800p-1 17 #define Ln2Lo 0x1.ef35793c76730p-45 18 #define OneMHfRt2Top \ 19 0x00095f62 /* top32(asuint64(1)) - top32(asuint64(sqrt(2)/2)). */ 20 #define OneTop12 0x3ff 21 #define HfRt2Top 0x3fe6a09e /* top32(asuint64(sqrt(2)/2)). */ 22 #define BottomMask 0xffffffff 23 #define C(i) __log1p_data.coeffs[i] 24 25 static inline double 26 log1p_inline (double x) 27 { 28 /* Helper for calculating log(1 + x) using order-18 polynomial on a reduced 29 interval. Copied from log1p_2u.c, with no special-case handling. See that 30 file for details of the algorithm. */ 31 double m = x + 1; 32 uint64_t mi = asuint64 (m); 33 34 /* Decompose x + 1 into (f + 1) * 2^k, with k chosen such that f is in 35 [sqrt(2)/2, sqrt(2)]. */ 36 uint32_t u = (mi >> 32) + OneMHfRt2Top; 37 int32_t k = (int32_t) (u >> 20) - OneTop12; 38 uint32_t utop = (u & 0x000fffff) + HfRt2Top; 39 uint64_t u_red = ((uint64_t) utop << 32) | (mi & BottomMask); 40 double f = asdouble (u_red) - 1; 41 42 /* Correction term for round-off in f. */ 43 double cm = (x - (m - 1)) / m; 44 45 /* Approximate log1p(f) with polynomial. */ 46 double f2 = f * f; 47 double f4 = f2 * f2; 48 double f8 = f4 * f4; 49 double p = fma (f, ESTRIN_18 (f, f2, f4, f8, f8 * f8, C) * f, f); 50 51 /* Recombine log1p(x) = k*log2 + log1p(f) + c/m. */ 52 double kd = k; 53 double y = fma (Ln2Lo, kd, cm); 54 return y + fma (Ln2Hi, kd, p); 55 } 56 57 /* Approximation for double-precision inverse tanh(x), using a simplified 58 version of log1p. Greatest observed error is 3.00 ULP: 59 atanh(0x1.e58f3c108d714p-4) got 0x1.e7da77672a647p-4 60 want 0x1.e7da77672a64ap-4. */ 61 double 62 atanh (double x) 63 { 64 uint64_t ix = asuint64 (x); 65 uint64_t sign = ix & ~AbsMask; 66 uint64_t ia = ix & AbsMask; 67 68 if (unlikely (ia == One)) 69 return __math_divzero (sign >> 32); 70 71 if (unlikely (ia > One)) 72 return __math_invalid (x); 73 74 double halfsign = asdouble (Half | sign); 75 double ax = asdouble (ia); 76 return halfsign * log1p_inline ((2 * ax) / (1 - ax)); 77 } 78 79 PL_SIG (S, D, 1, atanh, -1.0, 1.0) 80 PL_TEST_ULP (atanh, 3.00) 81 PL_TEST_INTERVAL (atanh, 0, 0x1p-23, 10000) 82 PL_TEST_INTERVAL (atanh, -0, -0x1p-23, 10000) 83 PL_TEST_INTERVAL (atanh, 0x1p-23, 1, 90000) 84 PL_TEST_INTERVAL (atanh, -0x1p-23, -1, 90000) 85 PL_TEST_INTERVAL (atanh, 1, inf, 100) 86 PL_TEST_INTERVAL (atanh, -1, -inf, 100) 87