1//===-- Single-precision general sinhf/coshf functions --------------------===//
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_SINHFCOSHF_UTILS_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_SINHFCOSHF_UTILS_H
11
12#include "exp10f_utils.h"
13#include "src/__support/FPUtil/multiply_add.h"
14
15namespace LIBC_NAMESPACE_DECL {
16
17namespace math {
18
19namespace sinhfcoshf_internal {
20
21// The function correctly calculates sinh(x) and cosh(x) by calculating exp(x)
22// and exp(-x) simultaneously.
23// To compute e^x, we perform the following range
24// reduction: find hi, mid, lo such that:
25// x = (hi + mid) * log(2) + lo, in which
26// hi is an integer,
27// 0 <= mid * 2^5 < 32 is an integer
28// -2^(-6) <= lo * log2(e) <= 2^-6.
29// In particular,
30// hi + mid = round(x * log2(e) * 2^5) * 2^(-5).
31// Then,
32// e^x = 2^(hi + mid) * e^lo = 2^hi * 2^mid * e^lo.
33// 2^mid is stored in the lookup table of 32 elements.
34// e^lo is computed using a degree-5 minimax polynomial
35// generated by Sollya:
36// e^lo ~ P(lo) = 1 + lo + c2 * lo^2 + ... + c5 * lo^5
37// = (1 + c2*lo^2 + c4*lo^4) + lo * (1 + c3*lo^2 + c5*lo^4)
38// = P_even + lo * P_odd
39// We perform 2^hi * 2^mid by simply add hi to the exponent field
40// of 2^mid.
41// To compute e^(-x), notice that:
42// e^(-x) = 2^(-(hi + mid)) * e^(-lo)
43// ~ 2^(-(hi + mid)) * P(-lo)
44// = 2^(-(hi + mid)) * (P_even - lo * P_odd)
45// So:
46// sinh(x) = (e^x - e^(-x)) / 2
47// ~ 0.5 * (2^(hi + mid) * (P_even + lo * P_odd) -
48// 2^(-(hi + mid)) * (P_even - lo * P_odd))
49// = 0.5 * (P_even * (2^(hi + mid) - 2^(-(hi + mid))) +
50// lo * P_odd * (2^(hi + mid) + 2^(-(hi + mid))))
51// And similarly:
52// cosh(x) = (e^x + e^(-x)) / 2
53// ~ 0.5 * (P_even * (2^(hi + mid) + 2^(-(hi + mid))) +
54// lo * P_odd * (2^(hi + mid) - 2^(-(hi + mid))))
55// The main point of these formulas is that the expensive part of calculating
56// the polynomials approximating lower parts of e^(x) and e^(-x) are shared
57// and only done once.
58template <bool is_sinh> LIBC_INLINE double exp_pm_eval(float x) {
59 double xd = static_cast<double>(x);
60
61 // kd = round(x * log2(e) * 2^5)
62 // k_p = round(x * log2(e) * 2^5)
63 // k_m = round(-x * log2(e) * 2^5)
64 double kd;
65 int k_p, k_m;
66
67#ifdef LIBC_TARGET_CPU_HAS_NEAREST_INT
68 kd = fputil::nearest_integer(ExpBase::LOG2_B * xd);
69 k_p = static_cast<int>(kd);
70 k_m = -k_p;
71#else
72 constexpr double HALF_WAY[2] = {0.5, -0.5};
73
74 k_p = static_cast<int>(
75 fputil::multiply_add(x: xd, y: ExpBase::LOG2_B, z: HALF_WAY[x < 0.0f]));
76 k_m = -k_p;
77 kd = static_cast<double>(k_p);
78#endif // LIBC_TARGET_CPU_HAS_NEAREST_INT
79
80 // hi = floor(kf * 2^(-5))
81 // exp_hi = shift hi to the exponent field of double precision.
82 int64_t exp_hi_p = static_cast<int64_t>((k_p >> ExpBase::MID_BITS))
83 << fputil::FPBits<double>::FRACTION_LEN;
84 int64_t exp_hi_m = static_cast<int64_t>((k_m >> ExpBase::MID_BITS))
85 << fputil::FPBits<double>::FRACTION_LEN;
86 // mh_p = 2^(hi + mid)
87 // mh_m = 2^(-(hi + mid))
88 // mh_bits_* = bit field of mh_*
89 int64_t mh_bits_p = ExpBase::EXP_2_MID[k_p & ExpBase::MID_MASK] + exp_hi_p;
90 int64_t mh_bits_m = ExpBase::EXP_2_MID[k_m & ExpBase::MID_MASK] + exp_hi_m;
91 double mh_p = fputil::FPBits<double>(uint64_t(mh_bits_p)).get_val();
92 double mh_m = fputil::FPBits<double>(uint64_t(mh_bits_m)).get_val();
93 // mh_sum = 2^(hi + mid) + 2^(-(hi + mid))
94 double mh_sum = mh_p + mh_m;
95 // mh_diff = 2^(hi + mid) - 2^(-(hi + mid))
96 double mh_diff = mh_p - mh_m;
97
98 // dx = lo = x - (hi + mid) * log(2)
99 double dx =
100 fputil::multiply_add(x: kd, y: ExpBase::M_LOGB_2_LO,
101 z: fputil::multiply_add(x: kd, y: ExpBase::M_LOGB_2_HI, z: xd));
102 double dx2 = dx * dx;
103
104 // c0 = 1 + COEFFS[0] * lo^2
105 // P_even = (1 + COEFFS[0] * lo^2 + COEFFS[2] * lo^4) / 2
106 double p_even = fputil::polyeval(x: dx2, a0: 0.5, a: ExpBase::COEFFS[0] * 0.5,
107 a: ExpBase::COEFFS[2] * 0.5);
108 // P_odd = (1 + COEFFS[1] * lo^2 + COEFFS[3] * lo^4) / 2
109 double p_odd = fputil::polyeval(x: dx2, a0: 0.5, a: ExpBase::COEFFS[1] * 0.5,
110 a: ExpBase::COEFFS[3] * 0.5);
111
112 double r;
113 if constexpr (is_sinh)
114 r = fputil::multiply_add(x: dx * mh_sum, y: p_odd, z: p_even * mh_diff);
115 else
116 r = fputil::multiply_add(x: dx * mh_diff, y: p_odd, z: p_even * mh_sum);
117 return r;
118}
119
120} // namespace sinhfcoshf_internal
121
122} // namespace math
123
124} // namespace LIBC_NAMESPACE_DECL
125
126#endif // LLVM_LIBC_SRC___SUPPORT_MATH_SINHFCOSHF_UTILS_H
127