mirror of
https://github.com/capstone-engine/llvm-capstone.git
synced 2024-11-23 13:50:11 +00:00
[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
This commit is contained in:
Tue Ly
parent
74f4daef04
commit
76bb278ebb
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -358,7 +358,7 @@ Higher Math Functions
|
||||
+------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
|
||||
| expl | | | | | | | | | | | | |
|
||||
+------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
|
||||
| exp10 | | | | | | | | | | | | |
|
||||
| exp10 | |check| | |check| | | |check| | |check| | | | |check| | | | | |
|
||||
+------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
|
||||
| exp10f | |check| | |check| | | |check| | |check| | | | |check| | | | | |
|
||||
+------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
|
||||
|
||||
@ -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>]>,
|
||||
]
|
||||
>;
|
||||
|
||||
@ -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)
|
||||
|
||||
18
libc/src/math/exp10.h
Normal file
18
libc/src/math/exp10.h
Normal file
@ -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
|
||||
@ -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
|
||||
|
||||
476
libc/src/math/generic/exp10.cpp
Normal file
476
libc/src/math/generic/exp10.cpp
Normal file
@ -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
|
||||
@ -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())
|
||||
|
||||
@ -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
|
||||
|
||||
150
libc/test/src/math/exp10_test.cpp
Normal file
150
libc/test/src/math/exp10_test.cpp
Normal file
@ -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);
|
||||
}
|
||||
@ -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 = [
|
||||
|
||||
@ -755,6 +755,13 @@ math_test(
|
||||
],
|
||||
)
|
||||
|
||||
math_test(
|
||||
name = "exp10",
|
||||
deps = [
|
||||
"//libc/utils/MPFRWrapper:mpfr_wrapper",
|
||||
],
|
||||
)
|
||||
|
||||
math_test(
|
||||
name = "fmod",
|
||||
hdrs = ["FModTest.h"],
|
||||
|
||||
Loading…
Reference in New Issue
Block a user