[libc-commits] [libc] 76bb278 - [libc][math] Implement double precision exp10 function correctly rounded for all rounding modes.
Tue Ly via libc-commits
libc-commits at lists.llvm.org
Wed Aug 30 05:43:58 PDT 2023
Author: Tue Ly
Date: 2023-08-30T08:43:50-04:00
New Revision: 76bb278ebb3e7207723ff6735dd9b3bf4653b85c
URL: https://github.com/llvm/llvm-project/commit/76bb278ebb3e7207723ff6735dd9b3bf4653b85c
DIFF: https://github.com/llvm/llvm-project/commit/76bb278ebb3e7207723ff6735dd9b3bf4653b85c.diff
LOG: [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
Added:
libc/src/math/exp10.h
libc/src/math/generic/exp10.cpp
libc/test/src/math/exp10_test.cpp
Modified:
libc/config/darwin/arm/entrypoints.txt
libc/config/linux/aarch64/entrypoints.txt
libc/config/linux/riscv64/entrypoints.txt
libc/config/linux/x86_64/entrypoints.txt
libc/config/windows/entrypoints.txt
libc/docs/math/index.rst
libc/spec/gnu_ext.td
libc/src/math/CMakeLists.txt
libc/src/math/generic/CMakeLists.txt
libc/src/math/generic/exp2.cpp
libc/test/src/math/CMakeLists.txt
utils/bazel/llvm-project-overlay/libc/BUILD.bazel
utils/bazel/llvm-project-overlay/libc/test/src/math/BUILD.bazel
Removed:
################################################################################
diff --git a/libc/config/darwin/arm/entrypoints.txt b/libc/config/darwin/arm/entrypoints.txt
index 10058865e0378e..f3d2cabb2c141d 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 7cd16df79a21c5..e0bf9800ec881c 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 4655744ebdcd01..2b2f2629f78ce6 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 4df1bf334c0c83..dcb8c6231432d3 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 15554b91eaf1c5..26552fe7b07b51 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 30c29252809b57..4e2ef031ab2444 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 a358caf4313b49..362add7283d6e1 100644
--- a/libc/spec/gnu_ext.td
+++ b/libc/spec/gnu_ext.td
@@ -31,6 +31,7 @@ def GnuExtensions : StandardSpec<"GNUExtensions"> {
RetValSpec<VoidType>,
[ArgSpec<FloatType>, ArgSpec<FloatPtr>, ArgSpec<FloatPtr>]
>,
+ FunctionSpec<"exp10", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
FunctionSpec<"exp10f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
]
>;
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index 79056fcb64b387..8b2021cac8239f 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 00000000000000..5afc95db2ce9c3
--- /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 fa7a5ed0ad9d52..9b092c3d3e7e65 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 00000000000000..4a43259b3307d3
--- /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 <errno.h>
+
+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<int>(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<int>(cpp::bit_cast<uint64_t>(tmp) >> 19);
+ double kd = static_cast<double>(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<double>(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<double>;
+ using FloatProp = typename fputil::FloatProperties<double>;
+ 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<double>(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<double>(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<double>(FPBits::inf());
+}
+
+LLVM_LIBC_FUNCTION(double, exp10, (double x)) {
+ using FPBits = typename fputil::FPBits<double>;
+ using FloatProp = typename fputil::FloatProperties<double>;
+ 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<int>(cpp::bit_cast<uint64_t>(tmp) >> 19);
+ double kd = static_cast<double>(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<int64_t>(hi) << FloatProp::MANTISSA_WIDTH;
+ double r = cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(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<int64_t>(hi) << FloatProp::MANTISSA_WIDTH;
+ double r = cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(upper_dd));
+ return r;
+ }
+
+ // Use 128-bit precision
+ Float128 r_f128 = exp10_f128(x, kd, idx1, idx2);
+
+ return static_cast<double>(r_f128);
+}
+
+} // namespace __llvm_libc
diff --git a/libc/src/math/generic/exp2.cpp b/libc/src/math/generic/exp2.cpp
index 6b66e95aef5903..96710c7d537069 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<double>;
@@ -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 ca55496b3244d2..04b1256e09b462 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 00000000000000..2b4c5ccc13f5d8
--- /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 <math.h>
+
+#include <errno.h>
+#include <stdint.h>
+
+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<double>(0.25).uintval();
+ uint64_t STOP = __llvm_libc::fputil::FPBits<double>(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<uint64_t>(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 54b7dce3e78843..a25e83691c98c1 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 bfada564199947..76b884714cfa95 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"],
More information about the libc-commits
mailing list