From 76bb278ebb3e7207723ff6735dd9b3bf4653b85c Mon Sep 17 00:00:00 2001 From: Tue Ly Date: Tue, 29 Aug 2023 18:08:41 -0400 Subject: [PATCH] [libc][math] Implement double precision exp10 function correctly rounded for all rounding modes. Implement double precision exp10 function correctly rounded for all rounding modes. Using the same algorithm as double precision exp (https://reviews.llvm.org/D158551) and exp2 (https://reviews.llvm.org/D158812) functions. Reviewed By: zimmermann6 Differential Revision: https://reviews.llvm.org/D159143 --- libc/config/darwin/arm/entrypoints.txt | 1 + libc/config/linux/aarch64/entrypoints.txt | 1 + libc/config/linux/riscv64/entrypoints.txt | 1 + libc/config/linux/x86_64/entrypoints.txt | 1 + libc/config/windows/entrypoints.txt | 1 + libc/docs/math/index.rst | 2 +- libc/spec/gnu_ext.td | 1 + libc/src/math/CMakeLists.txt | 1 + libc/src/math/exp10.h | 18 + libc/src/math/generic/CMakeLists.txt | 27 + libc/src/math/generic/exp10.cpp | 476 ++++++++++++++++++ libc/src/math/generic/exp2.cpp | 8 +- libc/test/src/math/CMakeLists.txt | 14 + libc/test/src/math/exp10_test.cpp | 150 ++++++ .../llvm-project-overlay/libc/BUILD.bazel | 16 + .../libc/test/src/math/BUILD.bazel | 7 + 16 files changed, 720 insertions(+), 5 deletions(-) create mode 100644 libc/src/math/exp10.h create mode 100644 libc/src/math/generic/exp10.cpp create mode 100644 libc/test/src/math/exp10_test.cpp diff --git a/libc/config/darwin/arm/entrypoints.txt b/libc/config/darwin/arm/entrypoints.txt index 10058865e037..f3d2cabb2c14 100644 --- a/libc/config/darwin/arm/entrypoints.txt +++ b/libc/config/darwin/arm/entrypoints.txt @@ -131,6 +131,7 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.erff libc.src.math.exp libc.src.math.expf + libc.src.math.exp10 libc.src.math.exp10f libc.src.math.exp2 libc.src.math.exp2f diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt index 7cd16df79a21..e0bf9800ec88 100644 --- a/libc/config/linux/aarch64/entrypoints.txt +++ b/libc/config/linux/aarch64/entrypoints.txt @@ -245,6 +245,7 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.erff libc.src.math.exp libc.src.math.expf + libc.src.math.exp10 libc.src.math.exp10f libc.src.math.exp2 libc.src.math.exp2f diff --git a/libc/config/linux/riscv64/entrypoints.txt b/libc/config/linux/riscv64/entrypoints.txt index 4655744ebdcd..2b2f2629f78c 100644 --- a/libc/config/linux/riscv64/entrypoints.txt +++ b/libc/config/linux/riscv64/entrypoints.txt @@ -254,6 +254,7 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.erff libc.src.math.exp libc.src.math.expf + libc.src.math.exp10 libc.src.math.exp10f libc.src.math.exp2 libc.src.math.exp2f diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt index 4df1bf334c0c..dcb8c6231432 100644 --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -258,6 +258,7 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.erff libc.src.math.exp libc.src.math.expf + libc.src.math.exp10 libc.src.math.exp10f libc.src.math.exp2 libc.src.math.exp2f diff --git a/libc/config/windows/entrypoints.txt b/libc/config/windows/entrypoints.txt index 15554b91eaf1..26552fe7b07b 100644 --- a/libc/config/windows/entrypoints.txt +++ b/libc/config/windows/entrypoints.txt @@ -130,6 +130,7 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.erff libc.src.math.exp libc.src.math.expf + libc.src.math.exp10 libc.src.math.exp10f libc.src.math.exp2 libc.src.math.exp2f diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst index 30c29252809b..4e2ef031ab24 100644 --- a/libc/docs/math/index.rst +++ b/libc/docs/math/index.rst @@ -358,7 +358,7 @@ Higher Math Functions +------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+ | expl | | | | | | | | | | | | | +------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+ -| exp10 | | | | | | | | | | | | | +| exp10 | |check| | |check| | | |check| | |check| | | | |check| | | | | | +------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+ | exp10f | |check| | |check| | | |check| | |check| | | | |check| | | | | | +------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+ diff --git a/libc/spec/gnu_ext.td b/libc/spec/gnu_ext.td index a358caf4313b..362add7283d6 100644 --- a/libc/spec/gnu_ext.td +++ b/libc/spec/gnu_ext.td @@ -31,6 +31,7 @@ def GnuExtensions : StandardSpec<"GNUExtensions"> { RetValSpec, [ArgSpec, ArgSpec, ArgSpec] >, + FunctionSpec<"exp10", RetValSpec, [ArgSpec]>, FunctionSpec<"exp10f", RetValSpec, [ArgSpec]>, ] >; diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt index 79056fcb64b3..8b2021cac823 100644 --- a/libc/src/math/CMakeLists.txt +++ b/libc/src/math/CMakeLists.txt @@ -85,6 +85,7 @@ add_math_entrypoint_object(expf) add_math_entrypoint_object(exp2) add_math_entrypoint_object(exp2f) +add_math_entrypoint_object(exp10) add_math_entrypoint_object(exp10f) add_math_entrypoint_object(expm1f) diff --git a/libc/src/math/exp10.h b/libc/src/math/exp10.h new file mode 100644 index 000000000000..5afc95db2ce9 --- /dev/null +++ b/libc/src/math/exp10.h @@ -0,0 +1,18 @@ +//===-- 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_MATH_EXP10_H +#define LLVM_LIBC_SRC_MATH_EXP10_H + +namespace __llvm_libc { + +double exp10(double x); + +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_MATH_EXP10_H diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index fa7a5ed0ad9d..9b092c3d3e7e 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -648,6 +648,33 @@ add_entrypoint_object( -O3 ) +add_entrypoint_object( + exp10 + SRCS + exp10.cpp + 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.macros.optimization + libc.include.errno + libc.src.errno.errno + libc.include.math + COMPILE_OPTIONS + -O3 +) + add_entrypoint_object( exp10f SRCS diff --git a/libc/src/math/generic/exp10.cpp b/libc/src/math/generic/exp10.cpp new file mode 100644 index 000000000000..4a43259b3307 --- /dev/null +++ b/libc/src/math/generic/exp10.cpp @@ -0,0 +1,476 @@ +//===-- Double-precision 10^x function ------------------------------------===// +// +// 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 +// +//===----------------------------------------------------------------------===// + +#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/macros/optimization.h" // LIBC_UNLIKELY + +#include + +namespace __llvm_libc { + +using fputil::DoubleDouble; +using fputil::TripleDouble; +using Float128 = typename fputil::DyadicFloat<128>; + +// 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; +constexpr double MLOG10_2_EXP2_M12_MID_32 = 0x1.0c0219dcp-51; +constexpr double MLOG10_2_EXP2_M12_LO = 0x1.da994fd20dba2p-87; + +// Error bounds: +// Errors when using double precision. +constexpr double ERR_D = 0x1.8p-63; + +// Errors when using double-double precision. +constexpr double ERR_DD = 0x1.8p-99; + +// 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; +} + +// 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) { + using MType = typename Float128::MantissaType; + + constexpr Float128 COEFFS_128[]{ + {false, -127, MType({0, 0x8000000000000000})}, // 1.0 + {false, -126, MType({0xea56d62b82d30a2d, 0x935d8dddaaa8ac16})}, + {false, -126, MType({0x80a99ce75f4d5bdb, 0xa9a92639e753443a})}, + {false, -126, MType({0x6a4f9d7dbf6c9635, 0x82382c8ef1652304})}, + {false, -124, MType({0x345787019216c7af, 0x12bd7609fd98c44c})}, + {false, -127, MType({0xcc41ed7e0d27aee5, 0x450a7ff47535d889})}, + {false, -130, MType({0x8326bb91a6e7601d, 0xd3f6b844702d636b})}, + {false, -130, MType({0xfa7b46df314112a9, 0x45b937f0d05bb1cd})}, + }; + + 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; +} + +// 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); + + 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); +} + +// 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; + using FloatProp = typename fputil::FloatProperties; + FPBits xbits(x); + + uint64_t x_u = xbits.uintval(); + uint64_t x_abs = x_u & FloatProp::EXP_MANT_MASK; + + // |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 static_cast(FPBits(FPBits::MIN_SUBNORMAL)); + 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 static_cast(FPBits(FPBits::MAX_NORMAL)); + + fputil::set_errno_if_required(ERANGE); + fputil::raise_except_if_required(FE_OVERFLOW); + } + // x is +inf or nan + return x + static_cast(FPBits::inf()); +} + +LLVM_LIBC_FUNCTION(double, exp10, (double x)) { + using FPBits = typename fputil::FPBits; + using FloatProp = typename fputil::FloatProperties; + 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); + + 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) << FloatProp::MANTISSA_WIDTH; + 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) << FloatProp::MANTISSA_WIDTH; + 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); +} + +} // namespace __llvm_libc diff --git a/libc/src/math/generic/exp2.cpp b/libc/src/math/generic/exp2.cpp index 6b66e95aef59..96710c7d5370 100644 --- a/libc/src/math/generic/exp2.cpp +++ b/libc/src/math/generic/exp2.cpp @@ -104,7 +104,7 @@ Float128 poly_approx_f128(const Float128 &dx) { return p; } -// Compute exp(x) using 128-bit precision. +// Compute 2^(x) using 128-bit precision. // TODO(lntue): investigate triple-double precision implementation for this // step. Float128 exp2_f128(double x, int hi, int idx1, int idx2) { @@ -192,7 +192,7 @@ double exp2_denorm(double x) { // Check for exceptional cases when: // * log2(1 - 2^-54) < x < log2(1 + 2^-53) // * x >= 1024 -// * x <= -1075 +// * x <= -1022 // * x is inf or nan double set_exceptional(double x) { using FPBits = typename fputil::FPBits; @@ -208,9 +208,9 @@ double set_exceptional(double x) { return fputil::multiply_add(x, 0.5, 1.0); } - // x <= 2^-1075 || x >= 1024 or inf/nan. + // x <= -1022 || x >= 1024 or inf/nan. if (x_u > 0xc08ff00000000000) { - // x <= 2^-1075 or -inf/nan + // x <= -1075 or -inf/nan if (x_u >= 0xc090cc0000000000) { // exp(-Inf) = 0 if (xbits.is_inf()) diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt index ca55496b3244..04b1256e09b4 100644 --- a/libc/test/src/math/CMakeLists.txt +++ b/libc/test/src/math/CMakeLists.txt @@ -647,6 +647,20 @@ add_fp_unittest( libc.src.__support.FPUtil.fp_bits ) +add_fp_unittest( + exp10_test + NEED_MPFR + SUITE + libc_math_unittests + SRCS + exp10_test.cpp + DEPENDS + libc.src.errno.errno + libc.include.math + libc.src.math.exp10 + libc.src.__support.FPUtil.fp_bits +) + add_fp_unittest( copysign_test SUITE diff --git a/libc/test/src/math/exp10_test.cpp b/libc/test/src/math/exp10_test.cpp new file mode 100644 index 000000000000..2b4c5ccc13f5 --- /dev/null +++ b/libc/test/src/math/exp10_test.cpp @@ -0,0 +1,150 @@ +//===-- Unittests for 10^x ------------------------------------------------===// +// +// 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 +// +//===----------------------------------------------------------------------===// + +#include "src/__support/FPUtil/FPBits.h" +#include "src/errno/libc_errno.h" +#include "src/math/exp10.h" +#include "test/UnitTest/FPMatcher.h" +#include "test/UnitTest/Test.h" +#include "utils/MPFRWrapper/MPFRUtils.h" +#include + +#include +#include + +namespace mpfr = __llvm_libc::testing::mpfr; +using __llvm_libc::testing::tlog; + +DECLARE_SPECIAL_CONSTANTS(double) + +TEST(LlvmLibcExp10Test, SpecialNumbers) { + EXPECT_FP_EQ(aNaN, __llvm_libc::exp10(aNaN)); + EXPECT_FP_EQ(inf, __llvm_libc::exp10(inf)); + EXPECT_FP_EQ_ALL_ROUNDING(zero, __llvm_libc::exp10(neg_inf)); + EXPECT_FP_EQ_WITH_EXCEPTION(zero, __llvm_libc::exp10(-0x1.0p20), + FE_UNDERFLOW); + EXPECT_FP_EQ_WITH_EXCEPTION(inf, __llvm_libc::exp10(0x1.0p20), FE_OVERFLOW); + EXPECT_FP_EQ_ALL_ROUNDING(1.0, __llvm_libc::exp10(0.0)); + EXPECT_FP_EQ_ALL_ROUNDING(1.0, __llvm_libc::exp10(-0.0)); +} + +TEST(LlvmLibcExp10Test, TrickyInputs) { + constexpr int N = 41; + constexpr uint64_t INPUTS[N] = { + 0x40033093317082F8, 0x3FD79289C6E6A5C0, + 0x3FD05DE80A173EA0, // 0x1.05de80a173eap-2 + 0xbf1eb7a4cb841fcc, // -0x1.eb7a4cb841fccp-14 + 0xbf19a61fb925970d, + 0x3fda7b764e2cf47a, // 0x1.a7b764e2cf47ap-2 + 0xc04757852a4b93aa, // -0x1.757852a4b93aap+5 + 0x4044c19e5712e377, // x=0x1.4c19e5712e377p+5 + 0xbf19a61fb925970d, // x=-0x1.9a61fb925970dp-14 + 0xc039a74cdab36c28, // x=-0x1.9a74cdab36c28p+4 + 0xc085b3e4e2e3bba9, // x=-0x1.5b3e4e2e3bba9p+9 + 0xc086960d591aec34, // x=-0x1.6960d591aec34p+9 + 0xc086232c09d58d91, // x=-0x1.6232c09d58d91p+9 + 0xc0874910d52d3051, // x=-0x1.74910d52d3051p9 + 0xc0867a172ceb0990, // x=-0x1.67a172ceb099p+9 + 0xc08ff80000000000, // x=-0x1.ff8p+9 + 0xbc971547652b82fe, // x=-0x1.71547652b82fep-54 + 0x0000000000000000, // x = 0 + 0x3ff0000000000000, // x = 1 + 0x4000000000000000, // x = 2 + 0x4008000000000000, // x = 3 + 0x4010000000000000, // x = 4 + 0x4014000000000000, // x = 5 + 0x4018000000000000, // x = 6 + 0x401c000000000000, // x = 7 + 0x4020000000000000, // x = 8 + 0x4022000000000000, // x = 9 + 0x4024000000000000, // x = 10 + 0x4026000000000000, // x = 11 + 0x4028000000000000, // x = 12 + 0x402a000000000000, // x = 13 + 0x402c000000000000, // x = 14 + 0x402e000000000000, // x = 15 + 0x4030000000000000, // x = 16 + 0x4031000000000000, // x = 17 + 0x4032000000000000, // x = 18 + 0x4033000000000000, // x = 19 + 0x4034000000000000, // x = 20 + 0x4035000000000000, // x = 21 + 0x4036000000000000, // x = 22 + 0x4037000000000000, // x = 23 + }; + for (int i = 0; i < N; ++i) { + double x = double(FPBits(INPUTS[i])); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp10, x, + __llvm_libc::exp10(x), 0.5); + } +} + +TEST(LlvmLibcExp10Test, InDoubleRange) { + constexpr uint64_t COUNT = 1'231; + uint64_t START = __llvm_libc::fputil::FPBits(0.25).uintval(); + uint64_t STOP = __llvm_libc::fputil::FPBits(4.0).uintval(); + uint64_t STEP = (STOP - START) / COUNT; + + auto test = [&](mpfr::RoundingMode rounding_mode) { + mpfr::ForceRoundingMode __r(rounding_mode); + if (!__r.success) + return; + + uint64_t fails = 0; + uint64_t count = 0; + uint64_t cc = 0; + double mx, mr = 0.0; + double tol = 0.5; + + for (uint64_t i = 0, v = START; i <= COUNT; ++i, v += STEP) { + double x = FPBits(v).get_val(); + if (isnan(x) || isinf(x) || x < 0.0) + continue; + libc_errno = 0; + double result = __llvm_libc::exp10(x); + ++cc; + if (isnan(result) || isinf(result)) + continue; + + ++count; + + if (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Exp10, x, result, + 0.5, rounding_mode)) { + ++fails; + while (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Exp10, x, + result, tol, rounding_mode)) { + mx = x; + mr = result; + + if (tol > 1000.0) + break; + + tol *= 2.0; + } + } + } + tlog << " Exp10 failed: " << fails << "/" << count << "/" << cc + << " tests.\n"; + tlog << " Max ULPs is at most: " << static_cast(tol) << ".\n"; + if (fails) { + EXPECT_MPFR_MATCH(mpfr::Operation::Exp10, mx, mr, 0.5, rounding_mode); + } + }; + + tlog << " Test Rounding To Nearest...\n"; + test(mpfr::RoundingMode::Nearest); + + tlog << " Test Rounding Downward...\n"; + test(mpfr::RoundingMode::Downward); + + tlog << " Test Rounding Upward...\n"; + test(mpfr::RoundingMode::Upward); + + tlog << " Test Rounding Toward Zero...\n"; + test(mpfr::RoundingMode::TowardZero); +} diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel index 54b7dce3e788..a25e83691c98 100644 --- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel +++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel @@ -1246,6 +1246,22 @@ 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_macros_optimization", + ":common_constants", + ":explogxf", + ], +) + libc_math_function( name = "exp10f", additional_deps = [ diff --git a/utils/bazel/llvm-project-overlay/libc/test/src/math/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/test/src/math/BUILD.bazel index bfada5641999..76b884714cfa 100644 --- a/utils/bazel/llvm-project-overlay/libc/test/src/math/BUILD.bazel +++ b/utils/bazel/llvm-project-overlay/libc/test/src/math/BUILD.bazel @@ -755,6 +755,13 @@ math_test( ], ) +math_test( + name = "exp10", + deps = [ + "//libc/utils/MPFRWrapper:mpfr_wrapper", + ], +) + math_test( name = "fmod", hdrs = ["FModTest.h"],