1//===-- Single-precision tanhf function -----------------------------------===//
2//
3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4// See https://llvm.org/LICENSE.txt for license information.
5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6//
7//===----------------------------------------------------------------------===//
8
9#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_TANHF_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_TANHF_H
11
12#include "exp10f_utils.h"
13#include "src/__support/FPUtil/FPBits.h"
14#include "src/__support/FPUtil/PolyEval.h"
15#include "src/__support/FPUtil/except_value_utils.h"
16#include "src/__support/FPUtil/multiply_add.h"
17#include "src/__support/FPUtil/nearest_integer.h"
18#include "src/__support/macros/config.h"
19#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
20#include "src/__support/macros/properties/cpu_features.h"
21
22namespace LIBC_NAMESPACE_DECL {
23
24namespace math {
25
26LIBC_INLINE float tanhf(float x) {
27 // 2^6 * log2(e)
28 constexpr double LOG2_E_EXP2_6 = ExpBase::LOG2_B * 2.0;
29
30 using FPBits = typename fputil::FPBits<float>;
31 FPBits xbits(x);
32 uint32_t x_abs = xbits.abs().uintval();
33
34 // When |x| >= 15, or x is inf or nan, or |x| <= 0.078125
35 if (LIBC_UNLIKELY((x_abs >= 0x4170'0000U) || (x_abs <= 0x3da0'0000U))) {
36 if (x_abs <= 0x3da0'0000U) {
37 // |x| <= 0.078125
38 if (LIBC_UNLIKELY(x_abs <= 0x3280'0000U)) {
39 // |x| <= 2^-26
40 return (x_abs != 0)
41 ? static_cast<float>(x - 0x1.5555555555555p-2 * x * x * x)
42 : x;
43 }
44
45 const double TAYLOR[] = {-0x1.5555555555555p-2, 0x1.1111111111111p-3,
46 -0x1.ba1ba1ba1ba1cp-5, 0x1.664f4882c10fap-6,
47 -0x1.226e355e6c23dp-7};
48 double xdbl = x;
49 double x2 = xdbl * xdbl;
50 // Taylor polynomial.
51 double x4 = x2 * x2;
52 double c0 = x2 * TAYLOR[0];
53 double c1 = fputil::multiply_add(x: x2, y: TAYLOR[2], z: TAYLOR[1]);
54 double c2 = fputil::multiply_add(x: x2, y: TAYLOR[4], z: TAYLOR[3]);
55 double pe = fputil::polyeval(x: x4, a0: c0, a: c1, a: c2);
56
57 return static_cast<float>(fputil::multiply_add(x: xdbl, y: pe, z: xdbl));
58 }
59
60 // |x| >= 15
61 if (LIBC_UNLIKELY(xbits.is_nan()))
62 return x + 1.0f; // sNaN to qNaN + signal
63
64 if (LIBC_UNLIKELY(xbits.is_inf()))
65 return xbits.is_neg() ? -1.0f : 1.0f;
66
67 if (xbits.is_pos())
68 return fputil::round_result_slightly_down(value_rn: 1.0f);
69
70 return fputil::round_result_slightly_up(value_rn: -1.0f);
71 }
72
73 // Range reduction: e^(2x) = 2^(hi + mid) * e^lo
74 // Let k = round( x * 2^6 * log2(e)),
75 // So k = (hi + mid) * 2^5
76 // Then lo = 2x - (hi + mid) * log(2) = 2x - k * 2^-5 * log(2).
77
78 double xd = static_cast<double>(x);
79 // k = round( x* 2^6 * log2(e) )
80 double k = 0;
81 // mk = -k
82 int mk = 0;
83#ifdef LIBC_TARGET_CPU_HAS_NEAREST_INT
84 k = fputil::nearest_integer(xd * LOG2_E_EXP2_6);
85 mk = -static_cast<int>(k);
86#else
87 const double half_way = xbits.is_neg() ? 0.5 : -0.5;
88
89 mk = static_cast<int>(fputil::multiply_add(x: xd, y: -LOG2_E_EXP2_6, z: half_way));
90 k = static_cast<double>(-mk);
91#endif // LIBC_TARGET_CPU_HAS_NEAREST_INT
92 // -hi = floor(-k * 2^(-MID_BITS))
93 // exp_mhi = shift -hi to the exponent field of double precision.
94 int64_t exp_mhi = static_cast<int64_t>(mk >> ExpBase::MID_BITS)
95 << fputil::FPBits<double>::FRACTION_LEN;
96 // mh = 2^(-hi - mid)
97 int64_t mh_bits = ExpBase::EXP_2_MID[mk & ExpBase::MID_MASK] + exp_mhi;
98 double mh = fputil::FPBits<double>(uint64_t(mh_bits)).get_val();
99 // dx = lo/2 = x - (hi + mid) * log(2)/2 = x - k * 2^-6 * log(2)
100 double dx = fputil::multiply_add(
101 x: k, y: ExpBase::M_LOGB_2_LO * 0.5,
102 z: fputil::multiply_add(x: k, y: ExpBase::M_LOGB_2_HI * 0.5, z: xd));
103
104 // > P = fpminimax(expm1(2*x)/x, 4, [|D...|], [-log(2)/128, log(2)/128]);
105 constexpr double COEFFS[] = {0x1.ffffffffe5bc8p0, 0x1.555555555cd67p0,
106 0x1.5555c2a9b48b4p-1, 0x1.11112a0e34bdbp-2};
107
108 double dx2 = dx * dx;
109 double c0 = fputil::multiply_add(x: dx, y: 2.0, z: 1.0);
110 double c1 = fputil::multiply_add(x: dx, y: COEFFS[1], z: COEFFS[0]);
111 double c2 = fputil::multiply_add(x: dx, y: COEFFS[3], z: COEFFS[2]);
112 double r = fputil::polyeval(x: dx2, a0: c0, a: c1, a: c2);
113
114 // tanh(x) = sinh(x) / cosh(x)
115 // = (e^x - e^(-x)) / (e^x + e^(-x))
116 // = (e^(2x) - 1) / (e^(2x) + 1)
117 // = (2^(hi + mid) * e^lo - 1) / (2^(hi + mid) * e^lo + 1)
118 // = (e^lo - 2^(-hi - mid)) / (e^lo + 2^(-hi - mid))
119 // = (r - mh) / (r + mh)
120 return static_cast<float>((r - mh) / (r + mh));
121}
122
123} // namespace math
124} // namespace LIBC_NAMESPACE_DECL
125
126#endif // LLVM_LIBC_SRC___SUPPORT_MATH_TANHF_H
127