diff --git a/libc/shared/math.h b/libc/shared/math.h index 3012cbb938816..b37aa46820523 100644 --- a/libc/shared/math.h +++ b/libc/shared/math.h @@ -12,6 +12,7 @@ #include "libc_common.h" #include "math/exp.h" +#include "math/exp10.h" #include "math/expf.h" #include "math/expf16.h" #include "math/frexpf.h" diff --git a/libc/shared/math/exp10.h b/libc/shared/math/exp10.h new file mode 100644 index 0000000000000..3d36d9103705f --- /dev/null +++ b/libc/shared/math/exp10.h @@ -0,0 +1,23 @@ +//===-- Shared exp10 function -----------------------------------*- 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_SHARED_MATH_EXP10_H +#define LLVM_LIBC_SHARED_MATH_EXP10_H + +#include "shared/libc_common.h" +#include "src/__support/math/exp10.h" + +namespace LIBC_NAMESPACE_DECL { +namespace shared { + +using math::exp10; + +} // namespace shared +} // namespace LIBC_NAMESPACE_DECL + +#endif // LLVM_LIBC_SHARED_MATH_EXP10_H diff --git a/libc/src/__support/FPUtil/double_double.h b/libc/src/__support/FPUtil/double_double.h index c27885aadc028..8e54e845de493 100644 --- a/libc/src/__support/FPUtil/double_double.h +++ b/libc/src/__support/FPUtil/double_double.h @@ -151,8 +151,8 @@ LIBC_INLINE DoubleDouble quick_mult(double a, const DoubleDouble &b) { } template -LIBC_INLINE DoubleDouble quick_mult(const DoubleDouble &a, - const DoubleDouble &b) { +LIBC_INLINE constexpr DoubleDouble quick_mult(const DoubleDouble &a, + const DoubleDouble &b) { DoubleDouble r = exact_mult(a.hi, b.hi); double t1 = multiply_add(a.hi, b.lo, r.lo); double t2 = multiply_add(a.lo, b.hi, t1); diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt index f7ef9e7694fe6..0bfc996c44fc8 100644 --- a/libc/src/__support/math/CMakeLists.txt +++ b/libc/src/__support/math/CMakeLists.txt @@ -149,3 +149,24 @@ add_header_library( libc.src.__support.integer_literals libc.src.__support.macros.optimization ) + +add_header_library( + exp10 + HDRS + exp10.h + DEPENDS + .exp_constants + .exp_utils + libc.src.__support.CPP.bit + libc.src.__support.CPP.optional + libc.src.__support.FPUtil.dyadic_float + libc.src.__support.FPUtil.fenv_impl + libc.src.__support.FPUtil.fp_bits + libc.src.__support.FPUtil.multiply_add + libc.src.__support.FPUtil.nearest_integer + libc.src.__support.FPUtil.polyeval + libc.src.__support.FPUtil.rounding_mode + libc.src.__support.FPUtil.triple_double + libc.src.__support.integer_literals + libc.src.__support.macros.optimization +) diff --git a/libc/src/__support/math/exp10.h b/libc/src/__support/math/exp10.h new file mode 100644 index 0000000000000..da94281c0c745 --- /dev/null +++ b/libc/src/__support/math/exp10.h @@ -0,0 +1,505 @@ +//===-- Implementation header for exp10 ------------------------*- 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_EXP10_H +#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP10_H + +#include "exp_constants.h" // Lookup tables EXP2_MID1 and EXP_M2. +#include "exp_utils.h" // ziv_test_denorm. +#include "src/__support/CPP/bit.h" +#include "src/__support/CPP/optional.h" +#include "src/__support/FPUtil/FEnvImpl.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/__support/FPUtil/PolyEval.h" +#include "src/__support/FPUtil/double_double.h" +#include "src/__support/FPUtil/dyadic_float.h" +#include "src/__support/FPUtil/multiply_add.h" +#include "src/__support/FPUtil/nearest_integer.h" +#include "src/__support/FPUtil/rounding_mode.h" +#include "src/__support/FPUtil/triple_double.h" +#include "src/__support/common.h" +#include "src/__support/integer_literals.h" +#include "src/__support/macros/config.h" +#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY + +namespace LIBC_NAMESPACE_DECL { + +using fputil::DoubleDouble; +using fputil::TripleDouble; +using Float128 = typename fputil::DyadicFloat<128>; + +using LIBC_NAMESPACE::operator""_u128; + +// log2(10) +constexpr double LOG2_10 = 0x1.a934f0979a371p+1; + +// -2^-12 * log10(2) +// > a = -2^-12 * log10(2); +// > b = round(a, 32, RN); +// > c = round(a - b, 32, RN); +// > d = round(a - b - c, D, RN); +// Errors < 1.5 * 2^-144 +constexpr double MLOG10_2_EXP2_M12_HI = -0x1.3441350ap-14; +constexpr double MLOG10_2_EXP2_M12_MID = 0x1.0c0219dc1da99p-51; + +#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS +constexpr double MLOG10_2_EXP2_M12_MID_32 = 0x1.0c0219dcp-51; +constexpr double MLOG10_2_EXP2_M12_LO = 0x1.da994fd20dba2p-87; +#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS + +// Error bounds: +// Errors when using double precision. +constexpr double ERR_D = 0x1.8p-63; + +#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS +// Errors when using double-double precision. +constexpr double ERR_DD = 0x1.8p-99; +#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS + +namespace { + +// Polynomial approximations with double precision. Generated by Sollya with: +// > P = fpminimax((10^x - 1)/x, 3, [|D...|], [-2^-14, 2^-14]); +// > P; +// Error bounds: +// | output - (10^dx - 1) / dx | < 2^-52. +LIBC_INLINE static constexpr double poly_approx_d(double dx) { + // dx^2 + double dx2 = dx * dx; + double c0 = + fputil::multiply_add(dx, 0x1.53524c73cea6ap+1, 0x1.26bb1bbb55516p+1); + double c1 = + fputil::multiply_add(dx, 0x1.2bd75cc6afc65p+0, 0x1.0470587aa264cp+1); + double p = fputil::multiply_add(dx2, c1, c0); + return p; +} + +#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS +// Polynomial approximation with double-double precision. Generated by Solya +// with: +// > P = fpminimax((10^x - 1)/x, 5, [|DD...|], [-2^-14, 2^-14]); +// Error bounds: +// | output - 10^(dx) | < 2^-101 +static constexpr DoubleDouble poly_approx_dd(const DoubleDouble &dx) { + // Taylor polynomial. + constexpr DoubleDouble COEFFS[] = { + {0, 0x1p0}, + {-0x1.f48ad494e927bp-53, 0x1.26bb1bbb55516p1}, + {-0x1.e2bfab3191cd2p-53, 0x1.53524c73cea69p1}, + {0x1.80fb65ec3b503p-53, 0x1.0470591de2ca4p1}, + {0x1.338fc05e21e55p-54, 0x1.2bd7609fd98c4p0}, + {0x1.d4ea116818fbp-56, 0x1.1429ffd519865p-1}, + {-0x1.872a8ff352077p-57, 0x1.a7ed70847c8b3p-3}, + + }; + + DoubleDouble p = fputil::polyeval(dx, COEFFS[0], COEFFS[1], COEFFS[2], + COEFFS[3], COEFFS[4], COEFFS[5], COEFFS[6]); + return p; +} + +// Polynomial approximation with 128-bit precision: +// Return exp(dx) ~ 1 + a0 * dx + a1 * dx^2 + ... + a6 * dx^7 +// For |dx| < 2^-14: +// | output - 10^dx | < 1.5 * 2^-124. +static constexpr Float128 poly_approx_f128(const Float128 &dx) { + constexpr Float128 COEFFS_128[]{ + {Sign::POS, -127, 0x80000000'00000000'00000000'00000000_u128}, // 1.0 + {Sign::POS, -126, 0x935d8ddd'aaa8ac16'ea56d62b'82d30a2d_u128}, + {Sign::POS, -126, 0xa9a92639'e753443a'80a99ce7'5f4d5bdb_u128}, + {Sign::POS, -126, 0x82382c8e'f1652304'6a4f9d7d'bf6c9635_u128}, + {Sign::POS, -124, 0x12bd7609'fd98c44c'34578701'9216c7af_u128}, + {Sign::POS, -127, 0x450a7ff4'7535d889'cc41ed7e'0d27aee5_u128}, + {Sign::POS, -130, 0xd3f6b844'702d636b'8326bb91'a6e7601d_u128}, + {Sign::POS, -130, 0x45b937f0'd05bb1cd'fa7b46df'314112a9_u128}, + }; + + Float128 p = fputil::polyeval(dx, COEFFS_128[0], COEFFS_128[1], COEFFS_128[2], + COEFFS_128[3], COEFFS_128[4], COEFFS_128[5], + COEFFS_128[6], COEFFS_128[7]); + return p; +} + +// Compute 10^(x) using 128-bit precision. +// TODO(lntue): investigate triple-double precision implementation for this +// step. +static constexpr Float128 exp10_f128(double x, double kd, int idx1, int idx2) { + double t1 = fputil::multiply_add(kd, MLOG10_2_EXP2_M12_HI, x); // exact + double t2 = kd * MLOG10_2_EXP2_M12_MID_32; // exact + double t3 = kd * MLOG10_2_EXP2_M12_LO; // Error < 2^-144 + + Float128 dx = fputil::quick_add( + Float128(t1), fputil::quick_add(Float128(t2), Float128(t3))); + + // TODO: Skip recalculating exp_mid1 and exp_mid2. + Float128 exp_mid1 = + fputil::quick_add(Float128(EXP2_MID1[idx1].hi), + fputil::quick_add(Float128(EXP2_MID1[idx1].mid), + Float128(EXP2_MID1[idx1].lo))); + + Float128 exp_mid2 = + fputil::quick_add(Float128(EXP2_MID2[idx2].hi), + fputil::quick_add(Float128(EXP2_MID2[idx2].mid), + Float128(EXP2_MID2[idx2].lo))); + + Float128 exp_mid = fputil::quick_mul(exp_mid1, exp_mid2); + + Float128 p = poly_approx_f128(dx); + + Float128 r = fputil::quick_mul(exp_mid, p); + + r.exponent += static_cast(kd) >> 12; + + return r; +} + +// Compute 10^x with double-double precision. +static constexpr DoubleDouble exp10_double_double(double x, double kd, + const DoubleDouble &exp_mid) { + // Recalculate dx: + // dx = x - k * 2^-12 * log10(2) + double t1 = fputil::multiply_add(kd, MLOG10_2_EXP2_M12_HI, x); // exact + double t2 = kd * MLOG10_2_EXP2_M12_MID_32; // exact + double t3 = kd * MLOG10_2_EXP2_M12_LO; // Error < 2^-140 + + DoubleDouble dx = fputil::exact_add(t1, t2); + dx.lo += t3; + + // Degree-6 polynomial approximation in double-double precision. + // | p - 10^x | < 2^-103. + DoubleDouble p = poly_approx_dd(dx); + + // Error bounds: 2^-102. + DoubleDouble r = fputil::quick_mult(exp_mid, p); + + return r; +} +#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS + +// When output is denormal. +static constexpr double exp10_denorm(double x) { + // Range reduction. + double tmp = fputil::multiply_add(x, LOG2_10, 0x1.8000'0000'4p21); + int k = static_cast(cpp::bit_cast(tmp) >> 19); + double kd = static_cast(k); + + uint32_t idx1 = (k >> 6) & 0x3f; + uint32_t idx2 = k & 0x3f; + + int hi = k >> 12; + + DoubleDouble exp_mid1{EXP2_MID1[idx1].mid, EXP2_MID1[idx1].hi}; + DoubleDouble exp_mid2{EXP2_MID2[idx2].mid, EXP2_MID2[idx2].hi}; + DoubleDouble exp_mid = fputil::quick_mult(exp_mid1, exp_mid2); + + // |dx| < 1.5 * 2^-15 + 2^-31 < 2^-14 + double lo_h = fputil::multiply_add(kd, MLOG10_2_EXP2_M12_HI, x); // exact + double dx = fputil::multiply_add(kd, MLOG10_2_EXP2_M12_MID, lo_h); + + double mid_lo = dx * exp_mid.hi; + + // Approximate (10^dx - 1)/dx ~ 1 + a0*dx + a1*dx^2 + a2*dx^3 + a3*dx^4. + double p = poly_approx_d(dx); + + double lo = fputil::multiply_add(p, mid_lo, exp_mid.lo); + +#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS + return ziv_test_denorm(hi, exp_mid.hi, lo, ERR_D) + .value(); +#else + if (auto r = ziv_test_denorm(hi, exp_mid.hi, lo, ERR_D); + LIBC_LIKELY(r.has_value())) + return r.value(); + + // Use double-double + DoubleDouble r_dd = exp10_double_double(x, kd, exp_mid); + + if (auto r = ziv_test_denorm(hi, r_dd.hi, r_dd.lo, ERR_DD); + LIBC_LIKELY(r.has_value())) + return r.value(); + + // Use 128-bit precision + Float128 r_f128 = exp10_f128(x, kd, idx1, idx2); + + return static_cast(r_f128); +#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS +} + +// Check for exceptional cases when: +// * log10(1 - 2^-54) < x < log10(1 + 2^-53) +// * x >= log10(2^1024) +// * x <= log10(2^-1022) +// * x is inf or nan +static constexpr double set_exceptional(double x) { + using FPBits = typename fputil::FPBits; + FPBits xbits(x); + + uint64_t x_u = xbits.uintval(); + uint64_t x_abs = xbits.abs().uintval(); + + // |x| < log10(1 + 2^-53) + if (x_abs <= 0x3c8bcb7b1526e50e) { + // 10^(x) ~ 1 + x/2 + return fputil::multiply_add(x, 0.5, 1.0); + } + + // x <= log10(2^-1022) || x >= log10(2^1024) or inf/nan. + if (x_u >= 0xc0733a7146f72a42) { + // x <= log10(2^-1075) or -inf/nan + if (x_u > 0xc07439b746e36b52) { + // exp(-Inf) = 0 + if (xbits.is_inf()) + return 0.0; + + // exp(nan) = nan + if (xbits.is_nan()) + return x; + + if (fputil::quick_get_round() == FE_UPWARD) + return FPBits::min_subnormal().get_val(); + fputil::set_errno_if_required(ERANGE); + fputil::raise_except_if_required(FE_UNDERFLOW); + return 0.0; + } + + return exp10_denorm(x); + } + + // x >= log10(2^1024) or +inf/nan + // x is finite + if (x_u < 0x7ff0'0000'0000'0000ULL) { + int rounding = fputil::quick_get_round(); + if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO) + return FPBits::max_normal().get_val(); + + fputil::set_errno_if_required(ERANGE); + fputil::raise_except_if_required(FE_OVERFLOW); + } + // x is +inf or nan + return x + FPBits::inf().get_val(); +} + +} // namespace + +namespace math { + +static constexpr double exp10(double x) { + using FPBits = typename fputil::FPBits; + FPBits xbits(x); + + uint64_t x_u = xbits.uintval(); + + // x <= log10(2^-1022) or x >= log10(2^1024) or + // log10(1 - 2^-54) < x < log10(1 + 2^-53). + if (LIBC_UNLIKELY(x_u >= 0xc0733a7146f72a42 || + (x_u <= 0xbc7bcb7b1526e50e && x_u >= 0x40734413509f79ff) || + x_u < 0x3c8bcb7b1526e50e)) { + return set_exceptional(x); + } + + // Now log10(2^-1075) < x <= log10(1 - 2^-54) or + // log10(1 + 2^-53) < x < log10(2^1024) + + // Range reduction: + // Let x = log10(2) * (hi + mid1 + mid2) + lo + // in which: + // hi is an integer + // mid1 * 2^6 is an integer + // mid2 * 2^12 is an integer + // then: + // 10^(x) = 2^hi * 2^(mid1) * 2^(mid2) * 10^(lo). + // With this formula: + // - multiplying by 2^hi is exact and cheap, simply by adding the exponent + // field. + // - 2^(mid1) and 2^(mid2) are stored in 2 x 64-element tables. + // - 10^(lo) ~ 1 + a0*lo + a1 * lo^2 + ... + // + // We compute (hi + mid1 + mid2) together by perform the rounding on + // x * log2(10) * 2^12. + // Since |x| < |log10(2^-1075)| < 2^9, + // |x * 2^12| < 2^9 * 2^12 < 2^21, + // So we can fit the rounded result round(x * 2^12) in int32_t. + // Thus, the goal is to be able to use an additional addition and fixed width + // shift to get an int32_t representing round(x * 2^12). + // + // Assuming int32_t using 2-complement representation, since the mantissa part + // of a double precision is unsigned with the leading bit hidden, if we add an + // extra constant C = 2^e1 + 2^e2 with e1 > e2 >= 2^23 to the product, the + // part that are < 2^e2 in resulted mantissa of (x*2^12*L2E + C) can be + // considered as a proper 2-complement representations of x*2^12. + // + // One small problem with this approach is that the sum (x*2^12 + C) in + // double precision is rounded to the least significant bit of the dorminant + // factor C. In order to minimize the rounding errors from this addition, we + // want to minimize e1. Another constraint that we want is that after + // shifting the mantissa so that the least significant bit of int32_t + // corresponds to the unit bit of (x*2^12*L2E), the sign is correct without + // any adjustment. So combining these 2 requirements, we can choose + // C = 2^33 + 2^32, so that the sign bit corresponds to 2^31 bit, and hence + // after right shifting the mantissa, the resulting int32_t has correct sign. + // With this choice of C, the number of mantissa bits we need to shift to the + // right is: 52 - 33 = 19. + // + // Moreover, since the integer right shifts are equivalent to rounding down, + // we can add an extra 0.5 so that it will become round-to-nearest, tie-to- + // +infinity. So in particular, we can compute: + // hmm = x * 2^12 + C, + // where C = 2^33 + 2^32 + 2^-1, then if + // k = int32_t(lower 51 bits of double(x * 2^12 + C) >> 19), + // the reduced argument: + // lo = x - log10(2) * 2^-12 * k is bounded by: + // |lo| = |x - log10(2) * 2^-12 * k| + // = log10(2) * 2^-12 * | x * log2(10) * 2^12 - k | + // <= log10(2) * 2^-12 * (2^-1 + 2^-19) + // < 1.5 * 2^-2 * (2^-13 + 2^-31) + // = 1.5 * (2^-15 * 2^-31) + // + // Finally, notice that k only uses the mantissa of x * 2^12, so the + // exponent 2^12 is not needed. So we can simply define + // C = 2^(33 - 12) + 2^(32 - 12) + 2^(-13 - 12), and + // k = int32_t(lower 51 bits of double(x + C) >> 19). + + // Rounding errors <= 2^-31. + double tmp = fputil::multiply_add(x, LOG2_10, 0x1.8000'0000'4p21); + int k = static_cast(cpp::bit_cast(tmp) >> 19); + double kd = static_cast(k); + + uint32_t idx1 = (k >> 6) & 0x3f; + uint32_t idx2 = k & 0x3f; + + int hi = k >> 12; + + DoubleDouble exp_mid1{EXP2_MID1[idx1].mid, EXP2_MID1[idx1].hi}; + DoubleDouble exp_mid2{EXP2_MID2[idx2].mid, EXP2_MID2[idx2].hi}; + DoubleDouble exp_mid = fputil::quick_mult(exp_mid1, exp_mid2); + + // |dx| < 1.5 * 2^-15 + 2^-31 < 2^-14 + double lo_h = fputil::multiply_add(kd, MLOG10_2_EXP2_M12_HI, x); // exact + double dx = fputil::multiply_add(kd, MLOG10_2_EXP2_M12_MID, lo_h); + + // We use the degree-4 polynomial to approximate 10^(lo): + // 10^(lo) ~ 1 + a0 * lo + a1 * lo^2 + a2 * lo^3 + a3 * lo^4 + // = 1 + lo * P(lo) + // So that the errors are bounded by: + // |P(lo) - (10^lo - 1)/lo| < |lo|^4 / 64 < 2^(-13 * 4) / 64 = 2^-58 + // Let P_ be an evaluation of P where all intermediate computations are in + // double precision. Using either Horner's or Estrin's schemes, the evaluated + // errors can be bounded by: + // |P_(lo) - P(lo)| < 2^-51 + // => |lo * P_(lo) - (2^lo - 1) | < 2^-65 + // => 2^(mid1 + mid2) * |lo * P_(lo) - expm1(lo)| < 2^-64. + // Since we approximate + // 2^(mid1 + mid2) ~ exp_mid.hi + exp_mid.lo, + // We use the expression: + // (exp_mid.hi + exp_mid.lo) * (1 + dx * P_(dx)) ~ + // ~ exp_mid.hi + (exp_mid.hi * dx * P_(dx) + exp_mid.lo) + // with errors bounded by 2^-64. + + double mid_lo = dx * exp_mid.hi; + + // Approximate (10^dx - 1)/dx ~ 1 + a0*dx + a1*dx^2 + a2*dx^3 + a3*dx^4. + double p = poly_approx_d(dx); + + double lo = fputil::multiply_add(p, mid_lo, exp_mid.lo); + +#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS + int64_t exp_hi = static_cast(hi) << FPBits::FRACTION_LEN; + double r = + cpp::bit_cast(exp_hi + cpp::bit_cast(exp_mid.hi + lo)); + return r; +#else + double upper = exp_mid.hi + (lo + ERR_D); + double lower = exp_mid.hi + (lo - ERR_D); + + if (LIBC_LIKELY(upper == lower)) { + // To multiply by 2^hi, a fast way is to simply add hi to the exponent + // field. + int64_t exp_hi = static_cast(hi) << FPBits::FRACTION_LEN; + double r = cpp::bit_cast(exp_hi + cpp::bit_cast(upper)); + return r; + } + + // Exact outputs when x = 1, 2, ..., 22 + hard to round with x = 23. + // Quick check mask: 0x800f'ffffU = ~(bits of 1.0 | ... | bits of 23.0) + if (LIBC_UNLIKELY((x_u & 0x8000'ffff'ffff'ffffULL) == 0ULL)) { + switch (x_u) { + case 0x3ff0000000000000: // x = 1.0 + return 10.0; + case 0x4000000000000000: // x = 2.0 + return 100.0; + case 0x4008000000000000: // x = 3.0 + return 1'000.0; + case 0x4010000000000000: // x = 4.0 + return 10'000.0; + case 0x4014000000000000: // x = 5.0 + return 100'000.0; + case 0x4018000000000000: // x = 6.0 + return 1'000'000.0; + case 0x401c000000000000: // x = 7.0 + return 10'000'000.0; + case 0x4020000000000000: // x = 8.0 + return 100'000'000.0; + case 0x4022000000000000: // x = 9.0 + return 1'000'000'000.0; + case 0x4024000000000000: // x = 10.0 + return 10'000'000'000.0; + case 0x4026000000000000: // x = 11.0 + return 100'000'000'000.0; + case 0x4028000000000000: // x = 12.0 + return 1'000'000'000'000.0; + case 0x402a000000000000: // x = 13.0 + return 10'000'000'000'000.0; + case 0x402c000000000000: // x = 14.0 + return 100'000'000'000'000.0; + case 0x402e000000000000: // x = 15.0 + return 1'000'000'000'000'000.0; + case 0x4030000000000000: // x = 16.0 + return 10'000'000'000'000'000.0; + case 0x4031000000000000: // x = 17.0 + return 100'000'000'000'000'000.0; + case 0x4032000000000000: // x = 18.0 + return 1'000'000'000'000'000'000.0; + case 0x4033000000000000: // x = 19.0 + return 10'000'000'000'000'000'000.0; + case 0x4034000000000000: // x = 20.0 + return 100'000'000'000'000'000'000.0; + case 0x4035000000000000: // x = 21.0 + return 1'000'000'000'000'000'000'000.0; + case 0x4036000000000000: // x = 22.0 + return 10'000'000'000'000'000'000'000.0; + case 0x4037000000000000: // x = 23.0 + return 0x1.52d02c7e14af6p76 + x; + } + } + + // Use double-double + DoubleDouble r_dd = exp10_double_double(x, kd, exp_mid); + + double upper_dd = r_dd.hi + (r_dd.lo + ERR_DD); + double lower_dd = r_dd.hi + (r_dd.lo - ERR_DD); + + if (LIBC_LIKELY(upper_dd == lower_dd)) { + // To multiply by 2^hi, a fast way is to simply add hi to the exponent + // field. + int64_t exp_hi = static_cast(hi) << FPBits::FRACTION_LEN; + double r = cpp::bit_cast(exp_hi + cpp::bit_cast(upper_dd)); + return r; + } + + // Use 128-bit precision + Float128 r_f128 = exp10_f128(x, kd, idx1, idx2); + + return static_cast(r_f128); +#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS +} + +} // namespace math + +} // namespace LIBC_NAMESPACE_DECL + +#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP10_H diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index b59beacd94143..352c2ad4ab22a 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -1457,20 +1457,7 @@ add_entrypoint_object( HDRS ../exp10.h DEPENDS - .common_constants - .explogxf - libc.src.__support.CPP.bit - libc.src.__support.CPP.optional - libc.src.__support.FPUtil.dyadic_float - libc.src.__support.FPUtil.fenv_impl - libc.src.__support.FPUtil.fp_bits - libc.src.__support.FPUtil.multiply_add - libc.src.__support.FPUtil.nearest_integer - libc.src.__support.FPUtil.polyeval - libc.src.__support.FPUtil.rounding_mode - libc.src.__support.FPUtil.triple_double - libc.src.__support.integer_literals - libc.src.__support.macros.optimization + libc.src.__support.math.exp10 libc.src.errno.errno ) diff --git a/libc/src/math/generic/exp10.cpp b/libc/src/math/generic/exp10.cpp index c464979b092c3..5c36d28c166ae 100644 --- a/libc/src/math/generic/exp10.cpp +++ b/libc/src/math/generic/exp10.cpp @@ -7,491 +7,10 @@ //===----------------------------------------------------------------------===// #include "src/math/exp10.h" -#include "common_constants.h" // Lookup tables EXP2_MID1 and EXP_M2. -#include "explogxf.h" // ziv_test_denorm. -#include "src/__support/CPP/bit.h" -#include "src/__support/CPP/optional.h" -#include "src/__support/FPUtil/FEnvImpl.h" -#include "src/__support/FPUtil/FPBits.h" -#include "src/__support/FPUtil/PolyEval.h" -#include "src/__support/FPUtil/double_double.h" -#include "src/__support/FPUtil/dyadic_float.h" -#include "src/__support/FPUtil/multiply_add.h" -#include "src/__support/FPUtil/nearest_integer.h" -#include "src/__support/FPUtil/rounding_mode.h" -#include "src/__support/FPUtil/triple_double.h" -#include "src/__support/common.h" -#include "src/__support/integer_literals.h" -#include "src/__support/macros/config.h" -#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY +#include "src/__support/math/exp10.h" namespace LIBC_NAMESPACE_DECL { -using fputil::DoubleDouble; -using fputil::TripleDouble; -using Float128 = typename fputil::DyadicFloat<128>; - -using LIBC_NAMESPACE::operator""_u128; - -// log2(10) -constexpr double LOG2_10 = 0x1.a934f0979a371p+1; - -// -2^-12 * log10(2) -// > a = -2^-12 * log10(2); -// > b = round(a, 32, RN); -// > c = round(a - b, 32, RN); -// > d = round(a - b - c, D, RN); -// Errors < 1.5 * 2^-144 -constexpr double MLOG10_2_EXP2_M12_HI = -0x1.3441350ap-14; -constexpr double MLOG10_2_EXP2_M12_MID = 0x1.0c0219dc1da99p-51; - -#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS -constexpr double MLOG10_2_EXP2_M12_MID_32 = 0x1.0c0219dcp-51; -constexpr double MLOG10_2_EXP2_M12_LO = 0x1.da994fd20dba2p-87; -#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS - -// Error bounds: -// Errors when using double precision. -constexpr double ERR_D = 0x1.8p-63; - -#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS -// Errors when using double-double precision. -constexpr double ERR_DD = 0x1.8p-99; -#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS - -namespace { - -// Polynomial approximations with double precision. Generated by Sollya with: -// > P = fpminimax((10^x - 1)/x, 3, [|D...|], [-2^-14, 2^-14]); -// > P; -// Error bounds: -// | output - (10^dx - 1) / dx | < 2^-52. -LIBC_INLINE double poly_approx_d(double dx) { - // dx^2 - double dx2 = dx * dx; - double c0 = - fputil::multiply_add(dx, 0x1.53524c73cea6ap+1, 0x1.26bb1bbb55516p+1); - double c1 = - fputil::multiply_add(dx, 0x1.2bd75cc6afc65p+0, 0x1.0470587aa264cp+1); - double p = fputil::multiply_add(dx2, c1, c0); - return p; -} - -#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS -// Polynomial approximation with double-double precision. Generated by Solya -// with: -// > P = fpminimax((10^x - 1)/x, 5, [|DD...|], [-2^-14, 2^-14]); -// Error bounds: -// | output - 10^(dx) | < 2^-101 -DoubleDouble poly_approx_dd(const DoubleDouble &dx) { - // Taylor polynomial. - constexpr DoubleDouble COEFFS[] = { - {0, 0x1p0}, - {-0x1.f48ad494e927bp-53, 0x1.26bb1bbb55516p1}, - {-0x1.e2bfab3191cd2p-53, 0x1.53524c73cea69p1}, - {0x1.80fb65ec3b503p-53, 0x1.0470591de2ca4p1}, - {0x1.338fc05e21e55p-54, 0x1.2bd7609fd98c4p0}, - {0x1.d4ea116818fbp-56, 0x1.1429ffd519865p-1}, - {-0x1.872a8ff352077p-57, 0x1.a7ed70847c8b3p-3}, - - }; - - DoubleDouble p = fputil::polyeval(dx, COEFFS[0], COEFFS[1], COEFFS[2], - COEFFS[3], COEFFS[4], COEFFS[5], COEFFS[6]); - return p; -} - -// Polynomial approximation with 128-bit precision: -// Return exp(dx) ~ 1 + a0 * dx + a1 * dx^2 + ... + a6 * dx^7 -// For |dx| < 2^-14: -// | output - 10^dx | < 1.5 * 2^-124. -Float128 poly_approx_f128(const Float128 &dx) { - constexpr Float128 COEFFS_128[]{ - {Sign::POS, -127, 0x80000000'00000000'00000000'00000000_u128}, // 1.0 - {Sign::POS, -126, 0x935d8ddd'aaa8ac16'ea56d62b'82d30a2d_u128}, - {Sign::POS, -126, 0xa9a92639'e753443a'80a99ce7'5f4d5bdb_u128}, - {Sign::POS, -126, 0x82382c8e'f1652304'6a4f9d7d'bf6c9635_u128}, - {Sign::POS, -124, 0x12bd7609'fd98c44c'34578701'9216c7af_u128}, - {Sign::POS, -127, 0x450a7ff4'7535d889'cc41ed7e'0d27aee5_u128}, - {Sign::POS, -130, 0xd3f6b844'702d636b'8326bb91'a6e7601d_u128}, - {Sign::POS, -130, 0x45b937f0'd05bb1cd'fa7b46df'314112a9_u128}, - }; - - Float128 p = fputil::polyeval(dx, COEFFS_128[0], COEFFS_128[1], COEFFS_128[2], - COEFFS_128[3], COEFFS_128[4], COEFFS_128[5], - COEFFS_128[6], COEFFS_128[7]); - return p; -} - -// Compute 10^(x) using 128-bit precision. -// TODO(lntue): investigate triple-double precision implementation for this -// step. -Float128 exp10_f128(double x, double kd, int idx1, int idx2) { - double t1 = fputil::multiply_add(kd, MLOG10_2_EXP2_M12_HI, x); // exact - double t2 = kd * MLOG10_2_EXP2_M12_MID_32; // exact - double t3 = kd * MLOG10_2_EXP2_M12_LO; // Error < 2^-144 - - Float128 dx = fputil::quick_add( - Float128(t1), fputil::quick_add(Float128(t2), Float128(t3))); - - // TODO: Skip recalculating exp_mid1 and exp_mid2. - Float128 exp_mid1 = - fputil::quick_add(Float128(EXP2_MID1[idx1].hi), - fputil::quick_add(Float128(EXP2_MID1[idx1].mid), - Float128(EXP2_MID1[idx1].lo))); - - Float128 exp_mid2 = - fputil::quick_add(Float128(EXP2_MID2[idx2].hi), - fputil::quick_add(Float128(EXP2_MID2[idx2].mid), - Float128(EXP2_MID2[idx2].lo))); - - Float128 exp_mid = fputil::quick_mul(exp_mid1, exp_mid2); - - Float128 p = poly_approx_f128(dx); - - Float128 r = fputil::quick_mul(exp_mid, p); - - r.exponent += static_cast(kd) >> 12; - - return r; -} - -// Compute 10^x with double-double precision. -DoubleDouble exp10_double_double(double x, double kd, - const DoubleDouble &exp_mid) { - // Recalculate dx: - // dx = x - k * 2^-12 * log10(2) - double t1 = fputil::multiply_add(kd, MLOG10_2_EXP2_M12_HI, x); // exact - double t2 = kd * MLOG10_2_EXP2_M12_MID_32; // exact - double t3 = kd * MLOG10_2_EXP2_M12_LO; // Error < 2^-140 - - DoubleDouble dx = fputil::exact_add(t1, t2); - dx.lo += t3; - - // Degree-6 polynomial approximation in double-double precision. - // | p - 10^x | < 2^-103. - DoubleDouble p = poly_approx_dd(dx); - - // Error bounds: 2^-102. - DoubleDouble r = fputil::quick_mult(exp_mid, p); - - return r; -} -#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS - -// When output is denormal. -double exp10_denorm(double x) { - // Range reduction. - double tmp = fputil::multiply_add(x, LOG2_10, 0x1.8000'0000'4p21); - int k = static_cast(cpp::bit_cast(tmp) >> 19); - double kd = static_cast(k); - - uint32_t idx1 = (k >> 6) & 0x3f; - uint32_t idx2 = k & 0x3f; - - int hi = k >> 12; - - DoubleDouble exp_mid1{EXP2_MID1[idx1].mid, EXP2_MID1[idx1].hi}; - DoubleDouble exp_mid2{EXP2_MID2[idx2].mid, EXP2_MID2[idx2].hi}; - DoubleDouble exp_mid = fputil::quick_mult(exp_mid1, exp_mid2); - - // |dx| < 1.5 * 2^-15 + 2^-31 < 2^-14 - double lo_h = fputil::multiply_add(kd, MLOG10_2_EXP2_M12_HI, x); // exact - double dx = fputil::multiply_add(kd, MLOG10_2_EXP2_M12_MID, lo_h); - - double mid_lo = dx * exp_mid.hi; - - // Approximate (10^dx - 1)/dx ~ 1 + a0*dx + a1*dx^2 + a2*dx^3 + a3*dx^4. - double p = poly_approx_d(dx); - - double lo = fputil::multiply_add(p, mid_lo, exp_mid.lo); - -#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS - return ziv_test_denorm(hi, exp_mid.hi, lo, ERR_D) - .value(); -#else - if (auto r = ziv_test_denorm(hi, exp_mid.hi, lo, ERR_D); - LIBC_LIKELY(r.has_value())) - return r.value(); - - // Use double-double - DoubleDouble r_dd = exp10_double_double(x, kd, exp_mid); - - if (auto r = ziv_test_denorm(hi, r_dd.hi, r_dd.lo, ERR_DD); - LIBC_LIKELY(r.has_value())) - return r.value(); - - // Use 128-bit precision - Float128 r_f128 = exp10_f128(x, kd, idx1, idx2); - - return static_cast(r_f128); -#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS -} - -// Check for exceptional cases when: -// * log10(1 - 2^-54) < x < log10(1 + 2^-53) -// * x >= log10(2^1024) -// * x <= log10(2^-1022) -// * x is inf or nan -double set_exceptional(double x) { - using FPBits = typename fputil::FPBits; - FPBits xbits(x); - - uint64_t x_u = xbits.uintval(); - uint64_t x_abs = xbits.abs().uintval(); - - // |x| < log10(1 + 2^-53) - if (x_abs <= 0x3c8bcb7b1526e50e) { - // 10^(x) ~ 1 + x/2 - return fputil::multiply_add(x, 0.5, 1.0); - } - - // x <= log10(2^-1022) || x >= log10(2^1024) or inf/nan. - if (x_u >= 0xc0733a7146f72a42) { - // x <= log10(2^-1075) or -inf/nan - if (x_u > 0xc07439b746e36b52) { - // exp(-Inf) = 0 - if (xbits.is_inf()) - return 0.0; - - // exp(nan) = nan - if (xbits.is_nan()) - return x; - - if (fputil::quick_get_round() == FE_UPWARD) - return FPBits::min_subnormal().get_val(); - fputil::set_errno_if_required(ERANGE); - fputil::raise_except_if_required(FE_UNDERFLOW); - return 0.0; - } - - return exp10_denorm(x); - } - - // x >= log10(2^1024) or +inf/nan - // x is finite - if (x_u < 0x7ff0'0000'0000'0000ULL) { - int rounding = fputil::quick_get_round(); - if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO) - return FPBits::max_normal().get_val(); - - fputil::set_errno_if_required(ERANGE); - fputil::raise_except_if_required(FE_OVERFLOW); - } - // x is +inf or nan - return x + FPBits::inf().get_val(); -} - -} // namespace - -LLVM_LIBC_FUNCTION(double, exp10, (double x)) { - using FPBits = typename fputil::FPBits; - FPBits xbits(x); - - uint64_t x_u = xbits.uintval(); - - // x <= log10(2^-1022) or x >= log10(2^1024) or - // log10(1 - 2^-54) < x < log10(1 + 2^-53). - if (LIBC_UNLIKELY(x_u >= 0xc0733a7146f72a42 || - (x_u <= 0xbc7bcb7b1526e50e && x_u >= 0x40734413509f79ff) || - x_u < 0x3c8bcb7b1526e50e)) { - return set_exceptional(x); - } - - // Now log10(2^-1075) < x <= log10(1 - 2^-54) or - // log10(1 + 2^-53) < x < log10(2^1024) - - // Range reduction: - // Let x = log10(2) * (hi + mid1 + mid2) + lo - // in which: - // hi is an integer - // mid1 * 2^6 is an integer - // mid2 * 2^12 is an integer - // then: - // 10^(x) = 2^hi * 2^(mid1) * 2^(mid2) * 10^(lo). - // With this formula: - // - multiplying by 2^hi is exact and cheap, simply by adding the exponent - // field. - // - 2^(mid1) and 2^(mid2) are stored in 2 x 64-element tables. - // - 10^(lo) ~ 1 + a0*lo + a1 * lo^2 + ... - // - // We compute (hi + mid1 + mid2) together by perform the rounding on - // x * log2(10) * 2^12. - // Since |x| < |log10(2^-1075)| < 2^9, - // |x * 2^12| < 2^9 * 2^12 < 2^21, - // So we can fit the rounded result round(x * 2^12) in int32_t. - // Thus, the goal is to be able to use an additional addition and fixed width - // shift to get an int32_t representing round(x * 2^12). - // - // Assuming int32_t using 2-complement representation, since the mantissa part - // of a double precision is unsigned with the leading bit hidden, if we add an - // extra constant C = 2^e1 + 2^e2 with e1 > e2 >= 2^23 to the product, the - // part that are < 2^e2 in resulted mantissa of (x*2^12*L2E + C) can be - // considered as a proper 2-complement representations of x*2^12. - // - // One small problem with this approach is that the sum (x*2^12 + C) in - // double precision is rounded to the least significant bit of the dorminant - // factor C. In order to minimize the rounding errors from this addition, we - // want to minimize e1. Another constraint that we want is that after - // shifting the mantissa so that the least significant bit of int32_t - // corresponds to the unit bit of (x*2^12*L2E), the sign is correct without - // any adjustment. So combining these 2 requirements, we can choose - // C = 2^33 + 2^32, so that the sign bit corresponds to 2^31 bit, and hence - // after right shifting the mantissa, the resulting int32_t has correct sign. - // With this choice of C, the number of mantissa bits we need to shift to the - // right is: 52 - 33 = 19. - // - // Moreover, since the integer right shifts are equivalent to rounding down, - // we can add an extra 0.5 so that it will become round-to-nearest, tie-to- - // +infinity. So in particular, we can compute: - // hmm = x * 2^12 + C, - // where C = 2^33 + 2^32 + 2^-1, then if - // k = int32_t(lower 51 bits of double(x * 2^12 + C) >> 19), - // the reduced argument: - // lo = x - log10(2) * 2^-12 * k is bounded by: - // |lo| = |x - log10(2) * 2^-12 * k| - // = log10(2) * 2^-12 * | x * log2(10) * 2^12 - k | - // <= log10(2) * 2^-12 * (2^-1 + 2^-19) - // < 1.5 * 2^-2 * (2^-13 + 2^-31) - // = 1.5 * (2^-15 * 2^-31) - // - // Finally, notice that k only uses the mantissa of x * 2^12, so the - // exponent 2^12 is not needed. So we can simply define - // C = 2^(33 - 12) + 2^(32 - 12) + 2^(-13 - 12), and - // k = int32_t(lower 51 bits of double(x + C) >> 19). - - // Rounding errors <= 2^-31. - double tmp = fputil::multiply_add(x, LOG2_10, 0x1.8000'0000'4p21); - int k = static_cast(cpp::bit_cast(tmp) >> 19); - double kd = static_cast(k); - - uint32_t idx1 = (k >> 6) & 0x3f; - uint32_t idx2 = k & 0x3f; - - int hi = k >> 12; - - DoubleDouble exp_mid1{EXP2_MID1[idx1].mid, EXP2_MID1[idx1].hi}; - DoubleDouble exp_mid2{EXP2_MID2[idx2].mid, EXP2_MID2[idx2].hi}; - DoubleDouble exp_mid = fputil::quick_mult(exp_mid1, exp_mid2); - - // |dx| < 1.5 * 2^-15 + 2^-31 < 2^-14 - double lo_h = fputil::multiply_add(kd, MLOG10_2_EXP2_M12_HI, x); // exact - double dx = fputil::multiply_add(kd, MLOG10_2_EXP2_M12_MID, lo_h); - - // We use the degree-4 polynomial to approximate 10^(lo): - // 10^(lo) ~ 1 + a0 * lo + a1 * lo^2 + a2 * lo^3 + a3 * lo^4 - // = 1 + lo * P(lo) - // So that the errors are bounded by: - // |P(lo) - (10^lo - 1)/lo| < |lo|^4 / 64 < 2^(-13 * 4) / 64 = 2^-58 - // Let P_ be an evaluation of P where all intermediate computations are in - // double precision. Using either Horner's or Estrin's schemes, the evaluated - // errors can be bounded by: - // |P_(lo) - P(lo)| < 2^-51 - // => |lo * P_(lo) - (2^lo - 1) | < 2^-65 - // => 2^(mid1 + mid2) * |lo * P_(lo) - expm1(lo)| < 2^-64. - // Since we approximate - // 2^(mid1 + mid2) ~ exp_mid.hi + exp_mid.lo, - // We use the expression: - // (exp_mid.hi + exp_mid.lo) * (1 + dx * P_(dx)) ~ - // ~ exp_mid.hi + (exp_mid.hi * dx * P_(dx) + exp_mid.lo) - // with errors bounded by 2^-64. - - double mid_lo = dx * exp_mid.hi; - - // Approximate (10^dx - 1)/dx ~ 1 + a0*dx + a1*dx^2 + a2*dx^3 + a3*dx^4. - double p = poly_approx_d(dx); - - double lo = fputil::multiply_add(p, mid_lo, exp_mid.lo); - -#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS - int64_t exp_hi = static_cast(hi) << FPBits::FRACTION_LEN; - double r = - cpp::bit_cast(exp_hi + cpp::bit_cast(exp_mid.hi + lo)); - return r; -#else - double upper = exp_mid.hi + (lo + ERR_D); - double lower = exp_mid.hi + (lo - ERR_D); - - if (LIBC_LIKELY(upper == lower)) { - // To multiply by 2^hi, a fast way is to simply add hi to the exponent - // field. - int64_t exp_hi = static_cast(hi) << FPBits::FRACTION_LEN; - double r = cpp::bit_cast(exp_hi + cpp::bit_cast(upper)); - return r; - } - - // Exact outputs when x = 1, 2, ..., 22 + hard to round with x = 23. - // Quick check mask: 0x800f'ffffU = ~(bits of 1.0 | ... | bits of 23.0) - if (LIBC_UNLIKELY((x_u & 0x8000'ffff'ffff'ffffULL) == 0ULL)) { - switch (x_u) { - case 0x3ff0000000000000: // x = 1.0 - return 10.0; - case 0x4000000000000000: // x = 2.0 - return 100.0; - case 0x4008000000000000: // x = 3.0 - return 1'000.0; - case 0x4010000000000000: // x = 4.0 - return 10'000.0; - case 0x4014000000000000: // x = 5.0 - return 100'000.0; - case 0x4018000000000000: // x = 6.0 - return 1'000'000.0; - case 0x401c000000000000: // x = 7.0 - return 10'000'000.0; - case 0x4020000000000000: // x = 8.0 - return 100'000'000.0; - case 0x4022000000000000: // x = 9.0 - return 1'000'000'000.0; - case 0x4024000000000000: // x = 10.0 - return 10'000'000'000.0; - case 0x4026000000000000: // x = 11.0 - return 100'000'000'000.0; - case 0x4028000000000000: // x = 12.0 - return 1'000'000'000'000.0; - case 0x402a000000000000: // x = 13.0 - return 10'000'000'000'000.0; - case 0x402c000000000000: // x = 14.0 - return 100'000'000'000'000.0; - case 0x402e000000000000: // x = 15.0 - return 1'000'000'000'000'000.0; - case 0x4030000000000000: // x = 16.0 - return 10'000'000'000'000'000.0; - case 0x4031000000000000: // x = 17.0 - return 100'000'000'000'000'000.0; - case 0x4032000000000000: // x = 18.0 - return 1'000'000'000'000'000'000.0; - case 0x4033000000000000: // x = 19.0 - return 10'000'000'000'000'000'000.0; - case 0x4034000000000000: // x = 20.0 - return 100'000'000'000'000'000'000.0; - case 0x4035000000000000: // x = 21.0 - return 1'000'000'000'000'000'000'000.0; - case 0x4036000000000000: // x = 22.0 - return 10'000'000'000'000'000'000'000.0; - case 0x4037000000000000: // x = 23.0 - return 0x1.52d02c7e14af6p76 + x; - } - } - - // Use double-double - DoubleDouble r_dd = exp10_double_double(x, kd, exp_mid); - - double upper_dd = r_dd.hi + (r_dd.lo + ERR_DD); - double lower_dd = r_dd.hi + (r_dd.lo - ERR_DD); - - if (LIBC_LIKELY(upper_dd == lower_dd)) { - // To multiply by 2^hi, a fast way is to simply add hi to the exponent - // field. - int64_t exp_hi = static_cast(hi) << FPBits::FRACTION_LEN; - double r = cpp::bit_cast(exp_hi + cpp::bit_cast(upper_dd)); - return r; - } - - // Use 128-bit precision - Float128 r_f128 = exp10_f128(x, kd, idx1, idx2); - - return static_cast(r_f128); -#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS -} +LLVM_LIBC_FUNCTION(double, exp10, (double x)) { return math::exp10(x); } } // namespace LIBC_NAMESPACE_DECL diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel index 4ab0126291276..26fc8b4cf6543 100644 --- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel +++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel @@ -2245,6 +2245,24 @@ libc_support_library( ], ) +libc_support_library( + name = "__support_math_exp10", + hdrs = ["src/__support/math/exp10.h"], + deps = [ + ":__support_math_exp_constants", + ":__support_math_exp_utils", + ":__support_fputil_double_double", + ":__support_fputil_dyadic_float", + ":__support_fputil_multiply_add", + ":__support_fputil_nearest_integer", + ":__support_fputil_polyeval", + ":__support_fputil_rounding_mode", + ":__support_fputil_triple_double", + ":__support_integer_literals", + ":__support_macros_optimization", + ], +) + ############################### complex targets ################################ libc_function( @@ -2849,17 +2867,8 @@ libc_math_function( libc_math_function( name = "exp10", additional_deps = [ - ":__support_fputil_double_double", - ":__support_fputil_dyadic_float", - ":__support_fputil_multiply_add", - ":__support_fputil_nearest_integer", - ":__support_fputil_polyeval", - ":__support_fputil_rounding_mode", - ":__support_fputil_triple_double", - ":__support_integer_literals", - ":__support_macros_optimization", - ":common_constants", - ":explogxf", + ":__support_math_exp10", + ":errno", ], )