//===-- Common utils for exp10f ---------------------------------*- C++ -*-===// // // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. // See https://llvm.org/LICENSE.txt for license information. // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception // //===----------------------------------------------------------------------===// #ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F_UTILS_H #define LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F_UTILS_H #include "src/__support/FPUtil/FPBits.h" #include "src/__support/FPUtil/PolyEval.h" #include "src/__support/FPUtil/nearest_integer.h" #include "src/__support/macros/config.h" namespace LIBC_NAMESPACE_DECL { struct ExpBase { // Base = e static constexpr int MID_BITS = 5; static constexpr int MID_MASK = (1 << MID_BITS) - 1; // log2(e) * 2^5 static constexpr double LOG2_B = 0x1.71547652b82fep+0 * (1 << MID_BITS); // High and low parts of -log(2) * 2^(-5) static constexpr double M_LOGB_2_HI = -0x1.62e42fefa0000p-1 / (1 << MID_BITS); static constexpr double M_LOGB_2_LO = -0x1.cf79abc9e3b3ap-40 / (1 << MID_BITS); // Look up table for bit fields of 2^(i/32) for i = 0..31, generated by Sollya // with: // > for i from 0 to 31 do printdouble(round(2^(i/32), D, RN)); static constexpr int64_t EXP_2_MID[1 << MID_BITS] = { 0x3ff0000000000000, 0x3ff059b0d3158574, 0x3ff0b5586cf9890f, 0x3ff11301d0125b51, 0x3ff172b83c7d517b, 0x3ff1d4873168b9aa, 0x3ff2387a6e756238, 0x3ff29e9df51fdee1, 0x3ff306fe0a31b715, 0x3ff371a7373aa9cb, 0x3ff3dea64c123422, 0x3ff44e086061892d, 0x3ff4bfdad5362a27, 0x3ff5342b569d4f82, 0x3ff5ab07dd485429, 0x3ff6247eb03a5585, 0x3ff6a09e667f3bcd, 0x3ff71f75e8ec5f74, 0x3ff7a11473eb0187, 0x3ff82589994cce13, 0x3ff8ace5422aa0db, 0x3ff93737b0cdc5e5, 0x3ff9c49182a3f090, 0x3ffa5503b23e255d, 0x3ffae89f995ad3ad, 0x3ffb7f76f2fb5e47, 0x3ffc199bdd85529c, 0x3ffcb720dcef9069, 0x3ffd5818dcfba487, 0x3ffdfc97337b9b5f, 0x3ffea4afa2a490da, 0x3fff50765b6e4540, }; // Approximating e^dx with degree-5 minimax polynomial generated by Sollya: // > Q = fpminimax(expm1(x)/x, 4, [|1, D...|], [-log(2)/64, log(2)/64]); // Then: // e^dx ~ P(dx) = 1 + dx + COEFFS[0] * dx^2 + ... + COEFFS[3] * dx^5. static constexpr double COEFFS[4] = { 0x1.ffffffffe5bc8p-2, 0x1.555555555cd67p-3, 0x1.5555c2a9b48b4p-5, 0x1.11112a0e34bdbp-7}; LIBC_INLINE static double powb_lo(double dx) { using fputil::multiply_add; double dx2 = dx * dx; double c0 = 1.0 + dx; // c1 = COEFFS[0] + COEFFS[1] * dx double c1 = multiply_add(dx, ExpBase::COEFFS[1], ExpBase::COEFFS[0]); // c2 = COEFFS[2] + COEFFS[3] * dx double c2 = multiply_add(dx, ExpBase::COEFFS[3], ExpBase::COEFFS[2]); // r = c4 + c5 * dx^4 // = 1 + dx + COEFFS[0] * dx^2 + ... + COEFFS[5] * dx^7 return fputil::polyeval(dx2, c0, c1, c2); } }; struct Exp10Base : public ExpBase { // log2(10) * 2^5 static constexpr double LOG2_B = 0x1.a934f0979a371p1 * (1 << MID_BITS); // High and low parts of -log10(2) * 2^(-5). // Notice that since |x * log2(10)| < 150: // |k| = |round(x * log2(10) * 2^5)| < 2^8 * 2^5 = 2^13 // So when the FMA instructions are not available, in order for the product // k * M_LOGB_2_HI // to be exact, we only store the high part of log10(2) up to 38 bits // (= 53 - 15) of precision. // It is generated by Sollya with: // > round(log10(2), 44, RN); static constexpr double M_LOGB_2_HI = -0x1.34413509f8p-2 / (1 << MID_BITS); // > round(log10(2) - 0x1.34413509f8p-2, D, RN); static constexpr double M_LOGB_2_LO = 0x1.80433b83b532ap-44 / (1 << MID_BITS); // Approximating 10^dx with degree-5 minimax polynomial generated by Sollya: // > Q = fpminimax((10^x - 1)/x, 4, [|D...|], [-log10(2)/2^6, log10(2)/2^6]); // Then: // 10^dx ~ P(dx) = 1 + COEFFS[0] * dx + ... + COEFFS[4] * dx^5. static constexpr double COEFFS[5] = {0x1.26bb1bbb55515p1, 0x1.53524c73bd3eap1, 0x1.0470591dff149p1, 0x1.2bd7c0a9fbc4dp0, 0x1.1429e74a98f43p-1}; LIBC_INLINE static double powb_lo(double dx) { using fputil::multiply_add; double dx2 = dx * dx; // c0 = 1 + COEFFS[0] * dx double c0 = multiply_add(dx, Exp10Base::COEFFS[0], 1.0); // c1 = COEFFS[1] + COEFFS[2] * dx double c1 = multiply_add(dx, Exp10Base::COEFFS[2], Exp10Base::COEFFS[1]); // c2 = COEFFS[3] + COEFFS[4] * dx double c2 = multiply_add(dx, Exp10Base::COEFFS[4], Exp10Base::COEFFS[3]); // r = c0 + dx^2 * (c1 + c2 * dx^2) // = c0 + c1 * dx^2 + c2 * dx^4 // = 1 + COEFFS[0] * dx + ... + COEFFS[4] * dx^5. return fputil::polyeval(dx2, c0, c1, c2); } }; // Output of range reduction for exp_b: (2^(mid + hi), lo) // where: // b^x = 2^(mid + hi) * b^lo struct exp_b_reduc_t { double mh; // 2^(mid + hi) double lo; }; // The function correctly calculates b^x value with at least float precision // in a limited range. // Range reduction: // b^x = 2^(hi + mid) * b^lo // where: // x = (hi + mid) * log_b(2) + lo // hi is an integer, // 0 <= mid * 2^MID_BITS < 2^MID_BITS is an integer // -2^(-MID_BITS - 1) <= lo * log2(b) <= 2^(-MID_BITS - 1) // Base class needs to provide the following constants: // - MID_BITS : number of bits after decimal points used for mid // - MID_MASK : 2^MID_BITS - 1, mask to extract mid bits // - LOG2_B : log2(b) * 2^MID_BITS for scaling // - M_LOGB_2_HI : high part of -log_b(2) * 2^(-MID_BITS) // - M_LOGB_2_LO : low part of -log_b(2) * 2^(-MID_BITS) // - EXP_2_MID : look up table for bit fields of 2^mid // Return: // { 2^(hi + mid), lo } template LIBC_INLINE static constexpr exp_b_reduc_t exp_b_range_reduc(float x) { double xd = static_cast(x); // kd = round((hi + mid) * log2(b) * 2^MID_BITS) double kd = fputil::nearest_integer(Base::LOG2_B * xd); // k = round((hi + mid) * log2(b) * 2^MID_BITS) int k = static_cast(kd); // hi = floor(kd * 2^(-MID_BITS)) // exp_hi = shift hi to the exponent field of double precision. uint64_t exp_hi = static_cast(k >> Base::MID_BITS) << fputil::FPBits::FRACTION_LEN; // mh = 2^hi * 2^mid // mh_bits = bit field of mh uint64_t mh_bits = Base::EXP_2_MID[k & Base::MID_MASK] + exp_hi; double mh = fputil::FPBits(mh_bits).get_val(); // dx = lo = x - (hi + mid) * log(2) double dx = fputil::multiply_add( kd, Base::M_LOGB_2_LO, fputil::multiply_add(kd, Base::M_LOGB_2_HI, xd)); return {mh, dx}; } } // namespace LIBC_NAMESPACE_DECL #endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F_UTILS_H