242 lines · c
1//===-- Common utilities for half-precision exponential 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_MATH_GENERIC_EXPXF16_H10#define LLVM_LIBC_SRC_MATH_GENERIC_EXPXF16_H11 12#include "hdr/stdint_proxy.h"13#include "src/__support/FPUtil/FPBits.h"14#include "src/__support/FPUtil/cast.h"15#include "src/__support/FPUtil/multiply_add.h"16#include "src/__support/FPUtil/nearest_integer.h"17#include "src/__support/macros/attributes.h"18#include "src/__support/macros/config.h"19#include "src/__support/math/exp10_float16_constants.h"20#include "src/__support/math/expf16_utils.h"21 22namespace LIBC_NAMESPACE_DECL {23 24namespace math {25 26namespace expxf16_internal {27 28LIBC_INLINE static ExpRangeReduction exp2_range_reduction(float16 x) {29 // For -25 < x < 16, to compute 2^x, we perform the following range reduction:30 // find hi, mid, lo, such that:31 // x = hi + mid + lo, in which32 // hi is an integer,33 // mid * 2^3 is an integer,34 // -2^(-4) <= lo < 2^(-4).35 // In particular,36 // hi + mid = round(x * 2^3) * 2^(-3).37 // Then,38 // 2^x = 2^(hi + mid + lo) = 2^hi * 2^mid * 2^lo.39 // We store 2^mid in the lookup table EXP2_MID_BITS, and compute 2^hi * 2^mid40 // by adding hi to the exponent field of 2^mid. 2^lo is computed using a41 // degree-3 minimax polynomial generated by Sollya.42 43 float xf = x;44 float kf = fputil::nearest_integer(xf * 0x1.0p+3f);45 int x_hi_mid = static_cast<int>(kf);46 unsigned x_hi = static_cast<unsigned>(x_hi_mid) >> 3;47 unsigned x_mid = static_cast<unsigned>(x_hi_mid) & 0x7;48 // lo = x - (hi + mid) = round(x * 2^3) * (-2^(-3)) + x49 float lo = fputil::multiply_add(kf, -0x1.0p-3f, xf);50 51 uint32_t exp2_hi_mid_bits =52 EXP2_MID_BITS[x_mid] +53 static_cast<uint32_t>(x_hi << fputil::FPBits<float>::FRACTION_LEN);54 float exp2_hi_mid = fputil::FPBits<float>(exp2_hi_mid_bits).get_val();55 // Degree-3 minimax polynomial generated by Sollya with the following56 // commands:57 // > display = hexadecimal;58 // > P = fpminimax((2^x - 1)/x, 2, [|SG...|], [-2^-4, 2^-4]);59 // > 1 + x * P;60 float exp2_lo = fputil::polyeval(lo, 0x1p+0f, 0x1.62e43p-1f, 0x1.ec0aa6p-3f,61 0x1.c6b4a6p-5f);62 return {exp2_hi_mid, exp2_lo};63}64 65// Generated by Sollya with the following commands:66// > display = hexadecimal;67// > round(log2(exp(1)), SG, RN);68static constexpr float LOG2F_E = 0x1.715476p+0f;69 70// Generated by Sollya with the following commands:71// > display = hexadecimal;72// > round(log(2), SG, RN);73static constexpr float LOGF_2 = 0x1.62e43p-1f;74 75// Generated by Sollya with the following commands:76// > display = hexadecimal;77// > for i from 0 to 31 do printsingle(round(2^(i * 2^-5), SG, RN));78static constexpr cpp::array<uint32_t, 32> EXP2_MID_5_BITS = {79 0x3f80'0000U, 0x3f82'cd87U, 0x3f85'aac3U, 0x3f88'980fU, 0x3f8b'95c2U,80 0x3f8e'a43aU, 0x3f91'c3d3U, 0x3f94'f4f0U, 0x3f98'37f0U, 0x3f9b'8d3aU,81 0x3f9e'f532U, 0x3fa2'7043U, 0x3fa5'fed7U, 0x3fa9'a15bU, 0x3fad'583fU,82 0x3fb1'23f6U, 0x3fb5'04f3U, 0x3fb8'fbafU, 0x3fbd'08a4U, 0x3fc1'2c4dU,83 0x3fc5'672aU, 0x3fc9'b9beU, 0x3fce'248cU, 0x3fd2'a81eU, 0x3fd7'44fdU,84 0x3fdb'fbb8U, 0x3fe0'ccdfU, 0x3fe5'b907U, 0x3fea'c0c7U, 0x3fef'e4baU,85 0x3ff5'257dU, 0x3ffa'83b3U,86};87 88// This function correctly calculates sinh(x) and cosh(x) by calculating exp(x)89// and exp(-x) simultaneously.90// To compute e^x, we perform the following range reduction:91// find hi, mid, lo such that:92// x = (hi + mid) * log(2) + lo, in which93// hi is an integer,94// 0 <= mid * 2^5 < 32 is an integer95// -2^(-5) <= lo * log2(e) <= 2^-5.96// In particular,97// hi + mid = round(x * log2(e) * 2^5) * 2^(-5).98// Then,99// e^x = 2^(hi + mid) * e^lo = 2^hi * 2^mid * e^lo.100// We store 2^mid in the lookup table EXP2_MID_5_BITS, and compute 2^hi * 2^mid101// by adding hi to the exponent field of 2^mid.102// e^lo is computed using a degree-3 minimax polynomial generated by Sollya:103// e^lo ~ P(lo)104// = 1 + lo + c2 * lo^2 + ... + c5 * lo^5105// = (1 + c2*lo^2 + c4*lo^4) + lo * (1 + c3*lo^2 + c5*lo^4)106// = P_even + lo * P_odd107// To compute e^(-x), notice that:108// e^(-x) = 2^(-(hi + mid)) * e^(-lo)109// ~ 2^(-(hi + mid)) * P(-lo)110// = 2^(-(hi + mid)) * (P_even - lo * P_odd)111// So:112// sinh(x) = (e^x - e^(-x)) / 2113// ~ 0.5 * (2^(hi + mid) * (P_even + lo * P_odd) -114// 2^(-(hi + mid)) * (P_even - lo * P_odd))115// = 0.5 * (P_even * (2^(hi + mid) - 2^(-(hi + mid))) +116// lo * P_odd * (2^(hi + mid) + 2^(-(hi + mid))))117// And similarly:118// cosh(x) = (e^x + e^(-x)) / 2119// ~ 0.5 * (P_even * (2^(hi + mid) + 2^(-(hi + mid))) +120// lo * P_odd * (2^(hi + mid) - 2^(-(hi + mid))))121// The main point of these formulas is that the expensive part of calculating122// the polynomials approximating lower parts of e^x and e^(-x) is shared and123// only done once.124template <bool IsSinh>125LIBC_INLINE static constexpr float16 eval_sinh_or_cosh(float16 x) {126 float xf = x;127 float kf = fputil::nearest_integer(xf * (LOG2F_E * 0x1.0p+5f));128 int x_hi_mid_p = static_cast<int>(kf);129 int x_hi_mid_m = -x_hi_mid_p;130 131 unsigned x_hi_p = static_cast<unsigned>(x_hi_mid_p) >> 5;132 unsigned x_hi_m = static_cast<unsigned>(x_hi_mid_m) >> 5;133 unsigned x_mid_p = static_cast<unsigned>(x_hi_mid_p) & 0x1f;134 unsigned x_mid_m = static_cast<unsigned>(x_hi_mid_m) & 0x1f;135 136 uint32_t exp2_hi_mid_bits_p =137 EXP2_MID_5_BITS[x_mid_p] +138 static_cast<uint32_t>(x_hi_p << fputil::FPBits<float>::FRACTION_LEN);139 uint32_t exp2_hi_mid_bits_m =140 EXP2_MID_5_BITS[x_mid_m] +141 static_cast<uint32_t>(x_hi_m << fputil::FPBits<float>::FRACTION_LEN);142 // exp2_hi_mid_p = 2^(hi + mid)143 float exp2_hi_mid_p = fputil::FPBits<float>(exp2_hi_mid_bits_p).get_val();144 // exp2_hi_mid_m = 2^(-(hi + mid))145 float exp2_hi_mid_m = fputil::FPBits<float>(exp2_hi_mid_bits_m).get_val();146 147 // exp2_hi_mid_sum = 2^(hi + mid) + 2^(-(hi + mid))148 float exp2_hi_mid_sum = exp2_hi_mid_p + exp2_hi_mid_m;149 // exp2_hi_mid_diff = 2^(hi + mid) - 2^(-(hi + mid))150 float exp2_hi_mid_diff = exp2_hi_mid_p - exp2_hi_mid_m;151 152 // lo = x - (hi + mid) = round(x * log2(e) * 2^5) * log(2) * (-2^(-5)) + x153 float lo = fputil::multiply_add(kf, LOGF_2 * -0x1.0p-5f, xf);154 float lo_sq = lo * lo;155 156 // Degree-3 minimax polynomial generated by Sollya with the following157 // commands:158 // > display = hexadecimal;159 // > P = fpminimax(expm1(x)/x, 2, [|SG...|], [-2^-5, 2^-5]);160 // > 1 + x * P;161 constexpr cpp::array<float, 4> COEFFS = {0x1p+0f, 0x1p+0f, 0x1.0004p-1f,162 0x1.555778p-3f};163 float half_p_odd =164 fputil::polyeval(lo_sq, COEFFS[1] * 0.5f, COEFFS[3] * 0.5f);165 float half_p_even =166 fputil::polyeval(lo_sq, COEFFS[0] * 0.5f, COEFFS[2] * 0.5f);167 168 // sinh(x) = lo * (0.5 * P_odd * (2^(hi + mid) + 2^(-(hi + mid)))) +169 // (0.5 * P_even * (2^(hi + mid) - 2^(-(hi + mid))))170 if constexpr (IsSinh)171 return fputil::cast<float16>(fputil::multiply_add(172 lo, half_p_odd * exp2_hi_mid_sum, half_p_even * exp2_hi_mid_diff));173 // cosh(x) = lo * (0.5 * P_odd * (2^(hi + mid) - 2^(-(hi + mid)))) +174 // (0.5 * P_even * (2^(hi + mid) + 2^(-(hi + mid))))175 return fputil::cast<float16>(fputil::multiply_add(176 lo, half_p_odd * exp2_hi_mid_diff, half_p_even * exp2_hi_mid_sum));177}178 179// Generated by Sollya with the following commands:180// > display = hexadecimal;181// > for i from 0 to 31 do print(round(log(1 + i * 2^-5), SG, RN));182static constexpr cpp::array<float, 32> LOGF_F = {183 0x0p+0f, 0x1.f829bp-6f, 0x1.f0a30cp-5f, 0x1.6f0d28p-4f,184 0x1.e27076p-4f, 0x1.29553p-3f, 0x1.5ff308p-3f, 0x1.9525aap-3f,185 0x1.c8ff7cp-3f, 0x1.fb9186p-3f, 0x1.1675cap-2f, 0x1.2e8e2cp-2f,186 0x1.4618bcp-2f, 0x1.5d1bdcp-2f, 0x1.739d8p-2f, 0x1.89a338p-2f,187 0x1.9f323ep-2f, 0x1.b44f78p-2f, 0x1.c8ff7cp-2f, 0x1.dd46ap-2f,188 0x1.f128f6p-2f, 0x1.02552ap-1f, 0x1.0be72ep-1f, 0x1.154c3ep-1f,189 0x1.1e85f6p-1f, 0x1.2795e2p-1f, 0x1.307d74p-1f, 0x1.393e0ep-1f,190 0x1.41d8fep-1f, 0x1.4a4f86p-1f, 0x1.52a2d2p-1f, 0x1.5ad404p-1f,191};192 193// Generated by Sollya with the following commands:194// > display = hexadecimal;195// > for i from 0 to 31 do print(round(log2(1 + i * 2^-5), SG, RN));196static constexpr cpp::array<float, 32> LOG2F_F = {197 0x0p+0f, 0x1.6bad38p-5f, 0x1.663f7p-4f, 0x1.08c588p-3f,198 0x1.5c01a4p-3f, 0x1.acf5e2p-3f, 0x1.fbc16cp-3f, 0x1.24407ap-2f,199 0x1.49a784p-2f, 0x1.6e221cp-2f, 0x1.91bba8p-2f, 0x1.b47ecp-2f,200 0x1.d6753ep-2f, 0x1.f7a856p-2f, 0x1.0c105p-1f, 0x1.1bf312p-1f,201 0x1.2b8034p-1f, 0x1.3abb4p-1f, 0x1.49a784p-1f, 0x1.584822p-1f,202 0x1.66a008p-1f, 0x1.74b1fep-1f, 0x1.82809ep-1f, 0x1.900e62p-1f,203 0x1.9d5dap-1f, 0x1.aa709p-1f, 0x1.b74948p-1f, 0x1.c3e9cap-1f,204 0x1.d053f6p-1f, 0x1.dc899ap-1f, 0x1.e88c6cp-1f, 0x1.f45e08p-1f,205};206 207// Generated by Sollya with the following commands:208// > display = hexadecimal;209// > for i from 0 to 31 do print(round(log10(1 + i * 2^-5), SG, RN));210static constexpr cpp::array<float, 32> LOG10F_F = {211 0x0p+0f, 0x1.b5e908p-7f, 0x1.af5f92p-6f, 0x1.3ed11ap-5f,212 0x1.a30a9ep-5f, 0x1.02428cp-4f, 0x1.31b306p-4f, 0x1.5fe804p-4f,213 0x1.8cf184p-4f, 0x1.b8de4ep-4f, 0x1.e3bc1ap-4f, 0x1.06cbd6p-3f,214 0x1.1b3e72p-3f, 0x1.2f3b6ap-3f, 0x1.42c7e8p-3f, 0x1.55e8c6p-3f,215 0x1.68a288p-3f, 0x1.7af974p-3f, 0x1.8cf184p-3f, 0x1.9e8e7cp-3f,216 0x1.afd3e4p-3f, 0x1.c0c514p-3f, 0x1.d1653p-3f, 0x1.e1b734p-3f,217 0x1.f1bdeep-3f, 0x1.00be06p-2f, 0x1.087a08p-2f, 0x1.101432p-2f,218 0x1.178da6p-2f, 0x1.1ee778p-2f, 0x1.2622bp-2f, 0x1.2d404cp-2f,219};220 221// Generated by Sollya with the following commands:222// > display = hexadecimal;223// > for i from 0 to 31 do print(round(1 / (1 + i * 2^-5), SG, RN));224static constexpr cpp::array<float, 32> ONE_OVER_F_F = {225 0x1p+0f, 0x1.f07c2p-1f, 0x1.e1e1e2p-1f, 0x1.d41d42p-1f,226 0x1.c71c72p-1f, 0x1.bacf92p-1f, 0x1.af286cp-1f, 0x1.a41a42p-1f,227 0x1.99999ap-1f, 0x1.8f9c18p-1f, 0x1.861862p-1f, 0x1.7d05f4p-1f,228 0x1.745d18p-1f, 0x1.6c16c2p-1f, 0x1.642c86p-1f, 0x1.5c9882p-1f,229 0x1.555556p-1f, 0x1.4e5e0ap-1f, 0x1.47ae14p-1f, 0x1.414142p-1f,230 0x1.3b13b2p-1f, 0x1.3521dp-1f, 0x1.2f684cp-1f, 0x1.29e412p-1f,231 0x1.24924ap-1f, 0x1.1f7048p-1f, 0x1.1a7b96p-1f, 0x1.15b1e6p-1f,232 0x1.111112p-1f, 0x1.0c9714p-1f, 0x1.08421p-1f, 0x1.041042p-1f,233};234 235} // namespace expxf16_internal236 237} // namespace math238 239} // namespace LIBC_NAMESPACE_DECL240 241#endif // LLVM_LIBC_SRC_MATH_GENERIC_EXPXF16_H242