1 /* 2 * Single-precision tanh(x) function. 3 * 4 * Copyright (c) 2022-2024, Arm Limited. 5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 6 */ 7 #include "math_config.h" 8 #include "test_sig.h" 9 #include "test_defs.h" 10 11 /* 0x1.205966p+3, above which tanhf rounds to 1 (or -1 for negative). */ 12 #define BoringBound 0x41102cb3 13 #define AbsMask 0x7fffffff 14 #define One 0x3f800000 15 16 #define Shift (0x1.8p23f) 17 #define InvLn2 (0x1.715476p+0f) 18 #define Ln2hi (0x1.62e4p-1f) 19 #define Ln2lo (0x1.7f7d1cp-20f) 20 21 #define C(i) __expm1f_poly[i] 22 23 static inline float 24 expm1f_inline (float x) 25 { 26 /* Helper routine for calculating exp(x) - 1. 27 Copied from expm1f_1u6.c, with several simplifications: 28 - No special-case handling for tiny or special values, instead return 29 early from the main routine. 30 - No special handling for large values: 31 - No early return for infinity. 32 - Simpler combination of p and t in final stage of algorithm. 33 - |i| < 27, so can calculate t by simpler shift-and-add, instead of 34 ldexpf (same as vector algorithm). */ 35 36 /* Reduce argument: f in [-ln2/2, ln2/2], i is exact. */ 37 float j = fmaf (InvLn2, x, Shift) - Shift; 38 int32_t i = j; 39 float f = fmaf (j, -Ln2hi, x); 40 f = fmaf (j, -Ln2lo, f); 41 42 /* Approximate expm1(f) with polynomial P, expm1(f) ~= f + f^2 * P(f). 43 Uses Estrin scheme, where the main expm1f routine uses Horner. */ 44 float f2 = f * f; 45 float p_01 = fmaf (f, C (1), C (0)); 46 float p_23 = fmaf (f, C (3), C (2)); 47 float p = fmaf (f2, p_23, p_01); 48 p = fmaf (f2 * f2, C (4), p); 49 p = fmaf (f2, p, f); 50 51 /* t = 2^i. */ 52 float t = asfloat ((uint32_t) (i + 127) << 23); 53 /* expm1(x) ~= p * t + (t - 1). */ 54 return fmaf (p, t, t - 1); 55 } 56 57 /* Approximation for single-precision tanh(x), using a simplified version of 58 expm1f. The maximum error is 2.58 ULP: 59 tanhf(0x1.fa5eep-5) got 0x1.f9ba02p-5 60 want 0x1.f9ba08p-5. */ 61 float 62 tanhf (float x) 63 { 64 uint32_t ix = asuint (x); 65 uint32_t iax = ix & AbsMask; 66 uint32_t sign = ix & ~AbsMask; 67 68 if (unlikely (iax > BoringBound)) 69 { 70 if (iax > 0x7f800000) 71 return __math_invalidf (x); 72 return asfloat (One | sign); 73 } 74 75 if (unlikely (iax < 0x34000000)) 76 return x; 77 78 /* tanh(x) = (e^2x - 1) / (e^2x + 1). */ 79 float q = expm1f_inline (2 * x); 80 return q / (q + 2); 81 } 82 83 TEST_SIG (S, F, 1, tanh, -10.0, 10.0) 84 TEST_ULP (tanhf, 2.09) 85 TEST_SYM_INTERVAL (tanhf, 0, 0x1p-23, 1000) 86 TEST_SYM_INTERVAL (tanhf, 0x1p-23, 0x1.205966p+3, 100000) 87 TEST_SYM_INTERVAL (tanhf, 0x1.205966p+3, inf, 100) 88