127 lines · c
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-exception6//7//===----------------------------------------------------------------------===//8 9#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_SINHFCOSHF_UTILS_H10#define LLVM_LIBC_SRC___SUPPORT_MATH_SINHFCOSHF_UTILS_H11 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 range24// reduction: find hi, mid, lo such that:25// x = (hi + mid) * log(2) + lo, in which26// hi is an integer,27// 0 <= mid * 2^5 < 32 is an integer28// -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 polynomial35// generated by Sollya:36// e^lo ~ P(lo) = 1 + lo + c2 * lo^2 + ... + c5 * lo^537// = (1 + c2*lo^2 + c4*lo^4) + lo * (1 + c3*lo^2 + c5*lo^4)38// = P_even + lo * P_odd39// We perform 2^hi * 2^mid by simply add hi to the exponent field40// 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)) / 247// ~ 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)) / 253// ~ 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 calculating56// the polynomials approximating lower parts of e^(x) and e^(-x) are shared57// 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_INT68 kd = fputil::nearest_integer(ExpBase::LOG2_B * xd);69 k_p = static_cast<int>(kd);70 k_m = -k_p;71#else72 constexpr double HALF_WAY[2] = {0.5, -0.5};73 74 k_p = static_cast<int>(75 fputil::multiply_add(xd, ExpBase::LOG2_B, HALF_WAY[x < 0.0f]));76 k_m = -k_p;77 kd = static_cast<double>(k_p);78#endif // LIBC_TARGET_CPU_HAS_NEAREST_INT79 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(kd, ExpBase::M_LOGB_2_LO,101 fputil::multiply_add(kd, ExpBase::M_LOGB_2_HI, xd));102 double dx2 = dx * dx;103 104 // c0 = 1 + COEFFS[0] * lo^2105 // P_even = (1 + COEFFS[0] * lo^2 + COEFFS[2] * lo^4) / 2106 double p_even = fputil::polyeval(dx2, 0.5, ExpBase::COEFFS[0] * 0.5,107 ExpBase::COEFFS[2] * 0.5);108 // P_odd = (1 + COEFFS[1] * lo^2 + COEFFS[3] * lo^4) / 2109 double p_odd = fputil::polyeval(dx2, 0.5, ExpBase::COEFFS[1] * 0.5,110 ExpBase::COEFFS[3] * 0.5);111 112 double r;113 if constexpr (is_sinh)114 r = fputil::multiply_add(dx * mh_sum, p_odd, p_even * mh_diff);115 else116 r = fputil::multiply_add(dx * mh_diff, p_odd, p_even * mh_sum);117 return r;118}119 120} // namespace sinhfcoshf_internal121 122} // namespace math123 124} // namespace LIBC_NAMESPACE_DECL125 126#endif // LLVM_LIBC_SRC___SUPPORT_MATH_SINHFCOSHF_UTILS_H127