1 /* 2 * Double-precision tanh(x) function. 3 * 4 * Copyright (c) 2023, Arm Limited. 5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 6 */ 7 #include "math_config.h" 8 #include "poly_scalar_f64.h" 9 #include "pl_sig.h" 10 #include "pl_test.h" 11 12 #define AbsMask 0x7fffffffffffffff 13 #define InvLn2 0x1.71547652b82fep0 14 #define Ln2hi 0x1.62e42fefa39efp-1 15 #define Ln2lo 0x1.abc9e3b39803fp-56 16 #define Shift 0x1.8p52 17 18 #define BoringBound 0x403241bf835f9d5f /* asuint64 (0x1.241bf835f9d5fp+4). */ 19 #define TinyBound 0x3e40000000000000 /* asuint64 (0x1p-27). */ 20 #define One 0x3ff0000000000000 21 22 static inline double 23 expm1_inline (double x) 24 { 25 /* Helper routine for calculating exp(x) - 1. Copied from expm1_2u5.c, with 26 several simplifications: 27 - No special-case handling for tiny or special values. 28 - Simpler combination of p and t in final stage of the algorithm. 29 - Use shift-and-add instead of ldexp to calculate t. */ 30 31 /* Reduce argument: f in [-ln2/2, ln2/2], i is exact. */ 32 double j = fma (InvLn2, x, Shift) - Shift; 33 int64_t i = j; 34 double f = fma (j, -Ln2hi, x); 35 f = fma (j, -Ln2lo, f); 36 37 /* Approximate expm1(f) using polynomial. */ 38 double f2 = f * f; 39 double f4 = f2 * f2; 40 double p = fma (f2, estrin_10_f64 (f, f2, f4, f4 * f4, __expm1_poly), f); 41 42 /* t = 2 ^ i. */ 43 double t = asdouble ((uint64_t) (i + 1023) << 52); 44 /* expm1(x) = p * t + (t - 1). */ 45 return fma (p, t, t - 1); 46 } 47 48 /* Approximation for double-precision tanh(x), using a simplified version of 49 expm1. The greatest observed error is 2.77 ULP: 50 tanh(-0x1.c4a4ca0f9f3b7p-3) got -0x1.bd6a21a163627p-3 51 want -0x1.bd6a21a163624p-3. */ 52 double 53 tanh (double x) 54 { 55 uint64_t ix = asuint64 (x); 56 uint64_t ia = ix & AbsMask; 57 uint64_t sign = ix & ~AbsMask; 58 59 if (unlikely (ia > BoringBound)) 60 { 61 if (ia > 0x7ff0000000000000) 62 return __math_invalid (x); 63 return asdouble (One | sign); 64 } 65 66 if (unlikely (ia < TinyBound)) 67 return x; 68 69 /* tanh(x) = (e^2x - 1) / (e^2x + 1). */ 70 double q = expm1_inline (2 * x); 71 return q / (q + 2); 72 } 73 74 PL_SIG (S, D, 1, tanh, -10.0, 10.0) 75 PL_TEST_ULP (tanh, 2.27) 76 PL_TEST_SYM_INTERVAL (tanh, 0, TinyBound, 1000) 77 PL_TEST_SYM_INTERVAL (tanh, TinyBound, BoringBound, 100000) 78 PL_TEST_SYM_INTERVAL (tanh, BoringBound, inf, 1000) 79