1*f3087befSAndrew Turner /* 2*f3087befSAndrew Turner * Single-precision atanh(x) function. 3*f3087befSAndrew Turner * 4*f3087befSAndrew Turner * Copyright (c) 2022-2024, Arm Limited. 5*f3087befSAndrew Turner * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 6*f3087befSAndrew Turner */ 7*f3087befSAndrew Turner 8*f3087befSAndrew Turner #include "math_config.h" 9*f3087befSAndrew Turner #include "mathlib.h" 10*f3087befSAndrew Turner #include "test_sig.h" 11*f3087befSAndrew Turner #include "test_defs.h" 12*f3087befSAndrew Turner 13*f3087befSAndrew Turner #define AbsMask 0x7fffffff 14*f3087befSAndrew Turner #define Half 0x3f000000 15*f3087befSAndrew Turner #define One 0x3f800000 16*f3087befSAndrew Turner #define Four 0x40800000 17*f3087befSAndrew Turner #define Ln2 0x1.62e43p-1f 18*f3087befSAndrew Turner /* asuint(0x1p-12), below which atanhf(x) rounds to x. */ 19*f3087befSAndrew Turner #define TinyBound 0x39800000 20*f3087befSAndrew Turner 21*f3087befSAndrew Turner #define C(i) __log1pf_data.coeffs[i] 22*f3087befSAndrew Turner 23*f3087befSAndrew Turner static inline float 24*f3087befSAndrew Turner eval_poly (float m) 25*f3087befSAndrew Turner { 26*f3087befSAndrew Turner /* Approximate log(1+m) on [-0.25, 0.5] using Estrin scheme. */ 27*f3087befSAndrew Turner float p_12 = fmaf (m, C (1), C (0)); 28*f3087befSAndrew Turner float p_34 = fmaf (m, C (3), C (2)); 29*f3087befSAndrew Turner float p_56 = fmaf (m, C (5), C (4)); 30*f3087befSAndrew Turner float p_78 = fmaf (m, C (7), C (6)); 31*f3087befSAndrew Turner 32*f3087befSAndrew Turner float m2 = m * m; 33*f3087befSAndrew Turner float p_02 = fmaf (m2, p_12, m); 34*f3087befSAndrew Turner float p_36 = fmaf (m2, p_56, p_34); 35*f3087befSAndrew Turner float p_79 = fmaf (m2, C (8), p_78); 36*f3087befSAndrew Turner 37*f3087befSAndrew Turner float m4 = m2 * m2; 38*f3087befSAndrew Turner float p_06 = fmaf (m4, p_36, p_02); 39*f3087befSAndrew Turner 40*f3087befSAndrew Turner return fmaf (m4 * p_79, m4, p_06); 41*f3087befSAndrew Turner } 42*f3087befSAndrew Turner 43*f3087befSAndrew Turner static inline float 44*f3087befSAndrew Turner log1pf_inline (float x) 45*f3087befSAndrew Turner { 46*f3087befSAndrew Turner /* Helper for calculating log(x + 1). Copied from log1pf_2u1.c, with no 47*f3087befSAndrew Turner special-case handling. See that file for details of the algorithm. */ 48*f3087befSAndrew Turner float m = x + 1.0f; 49*f3087befSAndrew Turner int k = (asuint (m) - 0x3f400000) & 0xff800000; 50*f3087befSAndrew Turner float s = asfloat (Four - k); 51*f3087befSAndrew Turner float m_scale = asfloat (asuint (x) - k) + fmaf (0.25f, s, -1.0f); 52*f3087befSAndrew Turner float p = eval_poly (m_scale); 53*f3087befSAndrew Turner float scale_back = (float) k * 0x1.0p-23f; 54*f3087befSAndrew Turner return fmaf (scale_back, Ln2, p); 55*f3087befSAndrew Turner } 56*f3087befSAndrew Turner 57*f3087befSAndrew Turner /* Approximation for single-precision inverse tanh(x), using a simplified 58*f3087befSAndrew Turner version of log1p. Maximum error is 3.08 ULP: 59*f3087befSAndrew Turner atanhf(0x1.ff0d5p-5) got 0x1.ffb768p-5 60*f3087befSAndrew Turner want 0x1.ffb76ep-5. */ 61*f3087befSAndrew Turner float 62*f3087befSAndrew Turner atanhf (float x) 63*f3087befSAndrew Turner { 64*f3087befSAndrew Turner uint32_t ix = asuint (x); 65*f3087befSAndrew Turner uint32_t iax = ix & AbsMask; 66*f3087befSAndrew Turner uint32_t sign = ix & ~AbsMask; 67*f3087befSAndrew Turner 68*f3087befSAndrew Turner if (unlikely (iax < TinyBound)) 69*f3087befSAndrew Turner return x; 70*f3087befSAndrew Turner 71*f3087befSAndrew Turner if (iax == One) 72*f3087befSAndrew Turner return __math_divzero (sign); 73*f3087befSAndrew Turner 74*f3087befSAndrew Turner if (unlikely (iax > One)) 75*f3087befSAndrew Turner return __math_invalidf (x); 76*f3087befSAndrew Turner 77*f3087befSAndrew Turner float halfsign = asfloat (Half | sign); 78*f3087befSAndrew Turner float ax = asfloat (iax); 79*f3087befSAndrew Turner return halfsign * log1pf_inline ((2 * ax) / (1 - ax)); 80*f3087befSAndrew Turner } 81*f3087befSAndrew Turner 82*f3087befSAndrew Turner TEST_SIG (S, F, 1, atanh, -1.0, 1.0) 83*f3087befSAndrew Turner TEST_ULP (atanhf, 2.59) 84*f3087befSAndrew Turner TEST_SYM_INTERVAL (atanhf, 0, 0x1p-12, 500) 85*f3087befSAndrew Turner TEST_SYM_INTERVAL (atanhf, 0x1p-12, 1, 200000) 86*f3087befSAndrew Turner TEST_SYM_INTERVAL (atanhf, 1, inf, 1000) 87