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