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