diff --git a/libc/src/__support/FPUtil/except_value_utils.h b/libc/src/__support/FPUtil/except_value_utils.h --- a/libc/src/__support/FPUtil/except_value_utils.h +++ b/libc/src/__support/FPUtil/except_value_utils.h @@ -12,6 +12,7 @@ #include "FEnvImpl.h" #include "FPBits.h" #include "src/__support/CPP/optional.h" +#include "src/__support/common.h" namespace __llvm_libc { @@ -47,7 +48,7 @@ Mapping values[N]; - constexpr cpp::optional lookup(UIntType x_bits) const { + LIBC_INLINE constexpr cpp::optional lookup(UIntType x_bits) const { for (size_t i = 0; i < N; ++i) { if (LIBC_UNLIKELY(x_bits == values[i].input)) { UIntType out_bits = values[i].rnd_towardzero_result; @@ -68,7 +69,8 @@ return cpp::nullopt; } - constexpr cpp::optional lookup_odd(UIntType x_abs, bool sign) const { + LIBC_INLINE constexpr cpp::optional lookup_odd(UIntType x_abs, + bool sign) const { for (size_t i = 0; i < N; ++i) { if (LIBC_UNLIKELY(x_abs == values[i].input)) { UIntType out_bits = values[i].rnd_towardzero_result; @@ -96,6 +98,19 @@ } }; +// Helper functions to set results for exceptional cases. +LIBC_INLINE float round_result_slightly_down(float value_rn) { + volatile float tmp = value_rn; + tmp = tmp - 0x1.0p-100f; + return tmp; +} + +LIBC_INLINE float round_result_slightly_up(float value_rn) { + volatile float tmp = value_rn; + tmp = tmp + 0x1.0p-100f; + return tmp; +} + } // namespace fputil } // namespace __llvm_libc diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -782,7 +782,7 @@ ../log10f.h DEPENDS .common_constants - libc.src.__support.FPUtil.basic_operations + libc.src.__support.FPUtil.except_value_utils libc.src.__support.FPUtil.fenv_impl libc.src.__support.FPUtil.fp_bits libc.src.__support.FPUtil.fma @@ -799,7 +799,7 @@ ../log1pf.h DEPENDS .common_constants - libc.src.__support.FPUtil.basic_operations + libc.src.__support.FPUtil.except_value_utils libc.src.__support.FPUtil.fenv_impl libc.src.__support.FPUtil.fp_bits libc.src.__support.FPUtil.fma @@ -816,7 +816,7 @@ ../log2f.h DEPENDS .common_constants - libc.src.__support.FPUtil.basic_operations + libc.src.__support.FPUtil.except_value_utils libc.src.__support.FPUtil.fenv_impl libc.src.__support.FPUtil.fp_bits libc.src.__support.FPUtil.fma @@ -833,10 +833,10 @@ ../logf.h DEPENDS .common_constants - libc.src.__support.FPUtil.basic_operations + libc.src.__support.FPUtil.except_value_utils libc.src.__support.FPUtil.fenv_impl libc.src.__support.FPUtil.fp_bits - libc.src.__support.FPUtil.fma + libc.src.__support.FPUtil.multiply_add libc.src.__support.FPUtil.polyeval COMPILE_OPTIONS -O3 diff --git a/libc/src/math/generic/log10f.cpp b/libc/src/math/generic/log10f.cpp --- a/libc/src/math/generic/log10f.cpp +++ b/libc/src/math/generic/log10f.cpp @@ -8,12 +8,14 @@ #include "src/math/log10f.h" #include "common_constants.h" // Lookup table for (1/f) -#include "src/__support/FPUtil/BasicOperations.h" #include "src/__support/FPUtil/FEnvImpl.h" #include "src/__support/FPUtil/FMA.h" #include "src/__support/FPUtil/FPBits.h" #include "src/__support/FPUtil/PolyEval.h" +#include "src/__support/FPUtil/except_value_utils.h" +#include "src/__support/FPUtil/multiply_add.h" #include "src/__support/common.h" +#include "src/__support/macros/attributes.h" // LIBC_UNLIKELY // This is an algorithm for log10(x) in single precision which is // correctly rounded for all rounding modes, based on the implementation of @@ -106,10 +108,10 @@ using FPBits = typename fputil::FPBits; FPBits xbits(x); - double m = 0.0; + uint32_t x_u = xbits.uintval(); // Exact powers of 10 and other hard-to-round cases. - switch (xbits.uintval()) { + switch (x_u) { case 0x4120'0000U: // x = 10 return 1.0f; case 0x42c8'0000U: // x = 100 @@ -131,37 +133,44 @@ case 0x5015'02f9U: // x = 10,000,000,000 return 10.0f; case 0x4f13'4f83U: // x = 2471461632.0 - if (fputil::get_round() == FE_UPWARD) - return 0x1.2c9314p+3f; - break; - case 0x7956'ba5eU: { // x = 69683218960000541503257137270226944.0 - int round_mode = fputil::get_round(); - if (round_mode == FE_DOWNWARD || round_mode == FE_TOWARDZERO) - return 0x1.16bebap+5f; - break; - } + return fputil::round_result_slightly_down(0x1.2c9314p+3f); + case 0x7956'ba5eU: // x = 69683218960000541503257137270226944.0 + return fputil::round_result_slightly_up(0x1.16bebap+5f); +#ifndef LIBC_TARGET_HAS_FMA + case 0x08ae'a356U: // x = 0x1.5d46acp-110f + return fputil::round_result_slightly_up(-0x1.07d3b4p+5f); + case 0x1c7d'a337U: // x = 0x1.fb466ep-71f + return fputil::round_result_slightly_up(-0x1.5137dp+4f); + case 0x69c8'c583U: // x = 0x1.918b06p+84f + return fputil::round_result_slightly_down(0x1.97b652p+4f); +#endif // LIBC_TARGET_HAS_FMA } - if (xbits.uintval() < FPBits::MIN_NORMAL || - xbits.uintval() > FPBits::MAX_NORMAL) { + int m = -FPBits::EXPONENT_BIAS; + + if (LIBC_UNLIKELY(x_u < FPBits::MIN_NORMAL || x_u > FPBits::MAX_NORMAL)) { if (xbits.is_zero()) { + // Return -inf and raise FE_DIVBYZERO + fputil::raise_except(FE_DIVBYZERO); return static_cast(FPBits::neg_inf()); } if (xbits.get_sign() && !xbits.is_nan()) { - return FPBits::build_nan(1 << (fputil::MantissaWidth::VALUE - 1)); + // Return NaN and raise FE_INVALID + fputil::raise_except(FE_INVALID); + return FPBits::build_quiet_nan(0); } if (xbits.is_inf_or_nan()) { return x; } // Normalize denormal inputs. xbits.set_val(xbits.get_val() * 0x1.0p23f); - m -= 23.0; + m -= 23; } - m += static_cast(xbits.get_exponent()); + m += xbits.get_unbiased_exponent(); + int f_index = xbits.get_mantissa() >> 16; // Set bits to 1.m xbits.set_unbiased_exponent(0x7F); - int f_index = xbits.get_mantissa() >> 16; FPBits f = xbits; f.bits &= ~0x0000'FFFF; @@ -169,7 +178,8 @@ double d = static_cast(xbits) - static_cast(f); d *= ONE_OVER_F[f_index]; - double extra_factor = fputil::multiply_add(m, LOG10_2, LOG10_F[f_index]); + double extra_factor = + fputil::multiply_add(static_cast(m), LOG10_2, LOG10_F[f_index]); double r = fputil::polyeval(d, extra_factor, 0x1.bcb7b1526e4c5p-2, -0x1.bcb7b1518a5e9p-3, 0x1.287a72a6f716p-3, diff --git a/libc/src/math/generic/log1pf.cpp b/libc/src/math/generic/log1pf.cpp --- a/libc/src/math/generic/log1pf.cpp +++ b/libc/src/math/generic/log1pf.cpp @@ -8,58 +8,59 @@ #include "src/math/log1pf.h" #include "common_constants.h" // Lookup table for (1/f) and log(f) -#include "src/__support/FPUtil/BasicOperations.h" #include "src/__support/FPUtil/FEnvImpl.h" #include "src/__support/FPUtil/FMA.h" #include "src/__support/FPUtil/FPBits.h" #include "src/__support/FPUtil/PolyEval.h" +#include "src/__support/FPUtil/except_value_utils.h" +#include "src/__support/FPUtil/multiply_add.h" #include "src/__support/common.h" +#include "src/__support/macros/attributes.h" // LIBC_UNLIKELY // This is an algorithm for log10(x) in single precision which is // correctly rounded for all rounding modes. // - An exhaustive test show that when x >= 2^45, log1pf(x) == logf(x) // for all rounding modes. -// - When 2^(-8) <= |x| < 2^45, the sum (double(x) + 1.0) is exact, +// - When 2^(-6) <= |x| < 2^45, the sum (double(x) + 1.0) is exact, // so we can adapt the correctly rounded algorithm of logf to compute // log(double(x) + 1.0) correctly. For more information about the logf // algorithm, see `libc/src/math/generic/logf.cpp`. -// - When |x| < 2^(-8), we use a degree-6 polynomial in double precision +// - When |x| < 2^(-6), we use a degree-8 polynomial in double precision // generated with Sollya using the following command: -// fpminimax(log(1 + x)/x, 5, [|D...|], [-2^-8; 2^-8]); +// fpminimax(log(1 + x)/x, 7, [|D...|], [-2^-6; 2^-6]); namespace __llvm_libc { namespace internal { -// We don't need to treat denormal -static inline float log(double x) { +// We don't need to treat denormal and 0 +LIBC_INLINE float log(double x) { constexpr double LOG_2 = 0x1.62e42fefa39efp-1; using FPBits = typename fputil::FPBits; FPBits xbits(x); - if (xbits.is_zero()) { - return static_cast(fputil::FPBits::neg_inf()); - } + uint64_t x_u = xbits.uintval(); - if (xbits.uintval() > FPBits::MAX_NORMAL) { + if (LIBC_UNLIKELY(x_u > FPBits::MAX_NORMAL)) { if (xbits.get_sign() && !xbits.is_nan()) { - return fputil::FPBits::build_nan( - 1 << (fputil::MantissaWidth::VALUE - 1)); + fputil::raise_except(FE_INVALID); + return fputil::FPBits::build_quiet_nan(0); } return static_cast(x); } double m = static_cast(xbits.get_exponent()); - // Set bits to 1.m - xbits.set_unbiased_exponent(0x3FF); // Get the 8 highest bits, use 7 bits (excluding the implicit hidden bit) for // lookup tables. int f_index = xbits.get_mantissa() >> 45; // fputil::MantissaWidth::VALUE - 7 + // Set bits to 1.m + xbits.set_unbiased_exponent(0x3FF); FPBits f = xbits; + // Clear the lowest 45 bits. f.bits &= ~0x0000'1FFF'FFFF'FFFFULL; @@ -80,87 +81,73 @@ LLVM_LIBC_FUNCTION(float, log1pf, (float x)) { using FPBits = typename fputil::FPBits; FPBits xbits(x); + uint32_t x_u = xbits.uintval(); + uint32_t x_a = x_u & 0x7fff'ffffU; double xd = static_cast(x); - if (xbits.get_exponent() >= -8) { + // Use log1p(x) = log(1 + x) for |x| > 2^-6; + if (x_a > 0x3c80'0000U) { // Hard-to-round cases. - switch (xbits.uintval()) { - case 0x3b9315c8U: // x = 0x1.262b9p-8f - if (fputil::get_round() != FE_UPWARD) - return 0x1.25830cp-8f; - break; - case 0x3c6eb7afU: // x = 0x1.dd6f5ep-7f - if (fputil::get_round() == FE_UPWARD) - return 0x1.d9fd86p-7f; - return 0x1.d9fd84p-7f; - case 0x41078febU: // x = 0x1.0f1fd6p+3f - if (fputil::get_round() != FE_UPWARD) - return 0x1.1fcbcep+1f; - break; + switch (x_u) { + case 0x41078febU: // x = 0x1.0f1fd6p3 + return fputil::round_result_slightly_up(0x1.1fcbcep1f); case 0x5cd69e88U: // x = 0x1.ad3d1p+58f - if (fputil::get_round() != FE_UPWARD) - return 0x1.45c146p+5f; - break; + return fputil::round_result_slightly_up(0x1.45c146p+5f); case 0x65d890d3U: // x = 0x1.b121a6p+76f - if (fputil::get_round() == FE_TONEAREST) - return 0x1.a9a3f2p+5f; - break; + return fputil::round_result_slightly_down(0x1.a9a3f2p+5f); case 0x6f31a8ecU: // x = 0x1.6351d8p+95f - if (fputil::get_round() == FE_TONEAREST) - return 0x1.08b512p+6f; - break; + return fputil::round_result_slightly_down(0x1.08b512p+6f); case 0x7a17f30aU: // x = 0x1.2fe614p+117f - if (fputil::get_round() != FE_UPWARD) - return 0x1.451436p+6f; - break; - case 0xbc4d092cU: // x = -0x1.9a1258p-7f - if (fputil::get_round() == FE_TONEAREST) - return -0x1.9ca8bep-7f; - break; - case 0xbc657728U: // x = -0x1.caee5p-7f - if (fputil::get_round() != FE_DOWNWARD) - return -0x1.ce2cccp-7f; - break; + return fputil::round_result_slightly_up(0x1.451436p+6f); case 0xbd1d20afU: // x = -0x1.3a415ep-5f - int round_mode = fputil::get_round(); - if (round_mode == FE_UPWARD || round_mode == FE_TOWARDZERO) - return -0x1.40711p-5f; - return -0x1.407112p-5f; + return fputil::round_result_slightly_up(-0x1.407112p-5f); + case 0xbf800000U: // x = -1.0 + fputil::raise_except(FE_DIVBYZERO); + return static_cast(fputil::FPBits::neg_inf()); +#ifndef LIBC_TARGET_HAS_FMA + case 0x4cc1c80bU: // x = 0x1.839016p+26f + return fputil::round_result_slightly_down(0x1.26fc04p+4f); + case 0x5ee8984eU: // x = 0x1.d1309cp+62f + return fputil::round_result_slightly_up(0x1.5c9442p+5f); + case 0x665e7ca6U: // x = 0x1.bcf94cp+77f + return fputil::round_result_slightly_up(0x1.af66cp+5f); + case 0x79e7ec37U: // x = 0x1.cfd86ep+116f + return fputil::round_result_slightly_up(0x1.43ff6ep+6); +#endif // LIBC_TARGET_HAS_FMA } return internal::log(xd + 1.0); } + // |x| <= 2^-6. // Hard-to round cases. - switch (xbits.uintval()) { + switch (x_u) { case 0x35400003U: // x = 0x1.800006p-21f - if (fputil::get_round() == FE_TONEAREST) - return 0x1.7ffffep-21f; - break; + return fputil::round_result_slightly_down(0x1.7ffffep-21f); case 0x3710001bU: // x = 0x1.200036p-17f - if (fputil::get_round() == FE_TONEAREST) - return 0x1.1fffe6p-17f; - break; - case 0xb53ffffdU: // x = -0x1.7ffffap-21f - if (fputil::get_round() != FE_DOWNWARD) - return -0x1.800002p-21f; - break; - case 0xb70fffe5U: // x = -0x1.1fffcap-17f - if (fputil::get_round() != FE_DOWNWARD) - return -0x1.20001ap-17f; - break; - case 0xbb0ec8c4U: // x = -0x1.1d9188p-9f - if (fputil::get_round() == FE_TONEAREST) - return -0x1.1de14ap-9f; - break; + return fputil::round_result_slightly_down(0x1.1fffe6p-17f); + case 0xb53ffffdU: // x = -0x1.7ffffap-21 + return fputil::round_result_slightly_down(-0x1.800002p-21f); + case 0xb70fffe5U: // x = -0x1.1fffcap-17 + return fputil::round_result_slightly_down(-0x1.20001ap-17f); + case 0xbb0ec8c4U: // x = -0x1.1d9188p-9 + return fputil::round_result_slightly_up(-0x1.1de14ap-9f); } - double r; - // Polymial generated with Sollya: - // > fpminimax(log(1 + x)/x, 5, [|D...|], [-2^-8; 2^-8]); - r = fputil::polyeval(xd, -0x1p-1, 0x1.5555555515551p-2, -0x1.ffffffff82bdap-3, - 0x1.999b33348d3aep-3, -0x1.5556cae3adcc3p-3); - return static_cast(fputil::multiply_add(r, xd * xd, xd)); + // Polymial generated by Sollya with: + // > fpminimax(log(1 + x)/x, 7, [|D...|], [-2^-6; 2^-6]); + const double COEFFS[7] = {-0x1.0000000000000p-1, 0x1.5555555556aadp-2, + -0x1.000000000181ap-2, 0x1.999998998124ep-3, + -0x1.55555452e2a2bp-3, 0x1.24adb8cde4aa7p-3, + -0x1.0019db915ef6fp-3}; + + double xsq = xd * xd; + double c0 = fputil::multiply_add(xd, COEFFS[1], COEFFS[0]); + double c1 = fputil::multiply_add(xd, COEFFS[3], COEFFS[2]); + double c2 = fputil::multiply_add(xd, COEFFS[5], COEFFS[4]); + double r = fputil::polyeval(xsq, xd, c0, c1, c2, COEFFS[6]); + + return r; } } // namespace __llvm_libc diff --git a/libc/src/math/generic/log2f.cpp b/libc/src/math/generic/log2f.cpp --- a/libc/src/math/generic/log2f.cpp +++ b/libc/src/math/generic/log2f.cpp @@ -8,12 +8,13 @@ #include "src/math/log2f.h" #include "common_constants.h" // Lookup table for (1/f) -#include "src/__support/FPUtil/BasicOperations.h" #include "src/__support/FPUtil/FEnvImpl.h" -#include "src/__support/FPUtil/FMA.h" #include "src/__support/FPUtil/FPBits.h" #include "src/__support/FPUtil/PolyEval.h" +#include "src/__support/FPUtil/except_value_utils.h" +#include "src/__support/FPUtil/multiply_add.h" #include "src/__support/common.h" +#include "src/__support/macros/attributes.h" // LIBC_UNLIKELY // This is a correctly-rounded algorithm for log2(x) in single precision with // round-to-nearest, tie-to-even mode from the RLIBM project at: @@ -101,52 +102,46 @@ LLVM_LIBC_FUNCTION(float, log2f, (float x)) { using FPBits = typename fputil::FPBits; FPBits xbits(x); - int m = 0; + uint32_t x_u = xbits.uintval(); // Hard to round value(s). - switch (FPBits(x).uintval()) { - case 0x3f81d0b5U: { - int rounding_mode = fputil::get_round(); - if (rounding_mode == FE_DOWNWARD || rounding_mode == FE_TOWARDZERO) { - return 0x1.4cdc4cp-6f; - } - break; - } - case 0x3f7e3274U: - if (fputil::get_round() == FE_TONEAREST) { - return -0x1.4e1d16p-7f; - } - break; - case 0x3f7d57f5U: - if (fputil::get_round() == FE_TOWARDZERO) { - return -0x1.ed1c32p-7f; - } - break; + using fputil::round_result_slightly_up; + + switch (x_u) { + case 0x3f7d57f5U: // x = 0x1.faafeap-1 + return round_result_slightly_up(-0x1.ed1c34p-7f); + case 0x3f7e3274U: // x = 0x1.fc64e8p-1f + return round_result_slightly_up(-0x1.4e1d16p-7f); + case 0x3f81d0b5U: // x = 0x1.03a16ap0f + return round_result_slightly_up(0x1.4cdc4cp-6f); } + int m = -FPBits::EXPONENT_BIAS; + // Exceptional inputs. - if (xbits.uintval() < FPBits::MIN_NORMAL || - xbits.uintval() > FPBits::MAX_NORMAL) { + if (LIBC_UNLIKELY(x_u < FPBits::MIN_NORMAL || x_u > FPBits::MAX_NORMAL)) { if (xbits.is_zero()) { + fputil::raise_except(FE_DIVBYZERO); return static_cast(FPBits::neg_inf()); } if (xbits.get_sign() && !xbits.is_nan()) { - return FPBits::build_nan(1 << (fputil::MantissaWidth::VALUE - 1)); + fputil::raise_except(FE_INVALID); + return FPBits::build_quiet_nan(0); } if (xbits.is_inf_or_nan()) { return x; } // Normalize denormal inputs. xbits.set_val(xbits.get_val() * 0x1.0p23f); - m = -23; + m -= 23; } - m += xbits.get_exponent(); + m += xbits.get_unbiased_exponent(); + int f_index = xbits.get_mantissa() >> 16; // Set bits to 1.m xbits.set_unbiased_exponent(0x7F); // Get the 8 highest bits, use 7 bits (excluding the implicit hidden bit) for // lookup tables. - int f_index = xbits.get_mantissa() >> 16; FPBits f = xbits; // Clear the lowest 16 bits. @@ -156,9 +151,17 @@ d *= ONE_OVER_F[f_index]; double extra_factor = static_cast(m) + LOG2_F[f_index]; - double r = __llvm_libc::fputil::polyeval( - d, extra_factor, 0x1.71547652bd4fp+0, -0x1.7154769b978c7p-1, - 0x1.ec71a99e349c8p-2, -0x1.720d90e6aac6cp-2, 0x1.5132da3583dap-2); + + const double COEFFS[5] = {0x1.71547652bd4fp+0, -0x1.7154769b978c7p-1, + 0x1.ec71a99e349c8p-2, -0x1.720d90e6aac6cp-2, + 0x1.5132da3583dap-2}; + + double dsq = d * d; + double c0 = fputil::multiply_add(d, COEFFS[0], extra_factor); + double c1 = fputil::multiply_add(d, COEFFS[2], COEFFS[1]); + double c2 = fputil::multiply_add(d, COEFFS[4], COEFFS[3]); + + double r = fputil::polyeval(dsq, c0, c1, c2); return static_cast(r); } diff --git a/libc/src/math/generic/logf.cpp b/libc/src/math/generic/logf.cpp --- a/libc/src/math/generic/logf.cpp +++ b/libc/src/math/generic/logf.cpp @@ -8,12 +8,13 @@ #include "src/math/logf.h" #include "common_constants.h" // Lookup table for (1/f) and log(f) -#include "src/__support/FPUtil/BasicOperations.h" #include "src/__support/FPUtil/FEnvImpl.h" -#include "src/__support/FPUtil/FMA.h" #include "src/__support/FPUtil/FPBits.h" #include "src/__support/FPUtil/PolyEval.h" +#include "src/__support/FPUtil/except_value_utils.h" +#include "src/__support/FPUtil/multiply_add.h" #include "src/__support/common.h" +#include "src/__support/macros/attributes.h" // LIBC_UNLIKELY // This is an algorithm for log(x) in single precision which is correctly // rounded for all rounding modes, based on the implementation of log(x) from @@ -53,64 +54,67 @@ constexpr double LOG_2 = 0x1.62e42fefa39efp-1; using FPBits = typename fputil::FPBits; FPBits xbits(x); + uint32_t x_u = xbits.uintval(); - switch (FPBits(x).uintval()) { + using fputil::round_result_slightly_down; + using fputil::round_result_slightly_up; + + switch (x_u) { + case 0x3f800001U: // x = 0x1.000002p+0f + return round_result_slightly_up(0x1.fffffep-24f); case 0x41178febU: // x = 0x1.2f1fd6p+3f - if (fputil::get_round() == FE_TONEAREST) - return 0x1.1fcbcep+1f; - break; + return round_result_slightly_up(0x1.1fcbcep+1f); case 0x4c5d65a5U: // x = 0x1.bacb4ap+25f - if (fputil::get_round() == FE_TONEAREST) - return 0x1.1e0696p+4f; - break; + return round_result_slightly_down(0x1.1e0696p+4f); + case 0x500ffb03U: // x = 0x1.1ff606p+33f + return round_result_slightly_up(0x1.6fdd34p+4f); + case 0x5cd69e88U: // x = 0x1.ad3d1p+58f + return round_result_slightly_up(0x1.45c146p+5f); case 0x65d890d3U: // x = 0x1.b121a6p+76f - if (fputil::get_round() == FE_TONEAREST) - return 0x1.a9a3f2p+5f; - break; + return round_result_slightly_down(0x1.a9a3f2p+5f); case 0x6f31a8ecU: // x = 0x1.6351d8p+95f - if (fputil::get_round() == FE_TONEAREST) - return 0x1.08b512p+6f; - break; - case 0x3f800001U: // x = 0x1.000002p+0f - if (fputil::get_round() == FE_UPWARD) - return 0x1p-23f; - return 0x1.fffffep-24f; - case 0x500ffb03U: // x = 0x1.1ff606p+33f - if (fputil::get_round() != FE_UPWARD) - return 0x1.6fdd34p+4f; - break; + return round_result_slightly_down(0x1.08b512p+6f); case 0x7a17f30aU: // x = 0x1.2fe614p+117f - if (fputil::get_round() != FE_UPWARD) - return 0x1.451436p+6f; - break; - case 0x5cd69e88U: // x = 0x1.ad3d1p+58f - if (fputil::get_round() != FE_UPWARD) - return 0x1.45c146p+5f; - break; + return round_result_slightly_up(0x1.451436p+6f); +#ifndef LIBC_TARGET_HAS_FMA + case 0x1b7679ffU: // x = 0x1.ecf3fep-73f + return round_result_slightly_up(-0x1.8f8e5ap+5f); + case 0x1e88452dU: // x = 0x1.108a5ap-66f + return round_result_slightly_up(-0x1.6d7b18p+5f); + case 0x5ee8984eU: // x = 0x1.d1309cp+62f; + return round_result_slightly_up(0x1.5c9442p+5f); + case 0x665e7ca6U: // x = 0x1.bcf94cp+77f + return round_result_slightly_up(0x1.af66cp+5f); + case 0x79e7ec37U: // x = 0x1.cfd86ep+116f + return round_result_slightly_up(0x1.43ff6ep+6f); +#endif // LIBC_TARGET_HAS_FMA } - int m = 0; + int m = -FPBits::EXPONENT_BIAS; - if (xbits.uintval() < FPBits::MIN_NORMAL || - xbits.uintval() > FPBits::MAX_NORMAL) { + if (LIBC_UNLIKELY(x_u < FPBits::MIN_NORMAL || x_u > FPBits::MAX_NORMAL)) { if (xbits.is_zero()) { + // Return -inf and raise FE_DIVBYZERO + fputil::raise_except(FE_DIVBYZERO); return static_cast(FPBits::neg_inf()); } if (xbits.get_sign() && !xbits.is_nan()) { - return FPBits::build_nan(1 << (fputil::MantissaWidth::VALUE - 1)); + // Return NaN and raise FE_INVALID + fputil::raise_except(FE_INVALID); + return FPBits::build_quiet_nan(0); } if (xbits.is_inf_or_nan()) { return x; } // Normalize denormal inputs. xbits.set_val(xbits.get_val() * 0x1.0p23f); - m = -23; + m -= 23; } - m += xbits.get_exponent(); + m += xbits.get_unbiased_exponent(); + int f_index = xbits.get_mantissa() >> 16; // Set bits to 1.m xbits.set_unbiased_exponent(0x7F); - int f_index = xbits.get_mantissa() >> 16; FPBits f = xbits; f.bits &= ~0x0000'FFFF; diff --git a/libc/test/UnitTest/FPMatcher.h b/libc/test/UnitTest/FPMatcher.h --- a/libc/test/UnitTest/FPMatcher.h +++ b/libc/test/UnitTest/FPMatcher.h @@ -83,6 +83,8 @@ __llvm_libc::fputil::testing::getMatcher<__llvm_libc::testing::Cond_EQ>( \ expected)) +#define EXPECT_FP_IS_NAN(actual) EXPECT_TRUE((actual) != (actual)) + #define ASSERT_FP_EQ(expected, actual) \ ASSERT_THAT( \ actual, \ @@ -133,6 +135,28 @@ } \ } while (0) +#define EXPECT_FP_EQ_WITH_EXCEPTION(expected_val, actual_val, expected_except) \ + do { \ + __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT); \ + EXPECT_FP_EQ(expected_val, actual_val); \ + if (math_errhandling & MATH_ERREXCEPT) { \ + EXPECT_EQ(__llvm_libc::fputil::test_except(FE_ALL_EXCEPT), \ + expected_except); \ + __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT); \ + } \ + } while (0) + +#define EXPECT_FP_IS_NAN_WITH_EXCEPTION(actual_val, expected_except) \ + do { \ + __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT); \ + EXPECT_FP_IS_NAN(actual_val); \ + if (math_errhandling & MATH_ERREXCEPT) { \ + EXPECT_EQ(__llvm_libc::fputil::test_except(FE_ALL_EXCEPT), \ + expected_except); \ + __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT); \ + } \ + } while (0) + #define EXPECT_FP_EQ_ALL_ROUNDING(expected, actual) \ do { \ using namespace __llvm_libc::testutils; \ diff --git a/libc/test/src/math/exhaustive/log1pf_test.cpp b/libc/test/src/math/exhaustive/log1pf_test.cpp --- a/libc/test/src/math/exhaustive/log1pf_test.cpp +++ b/libc/test/src/math/exhaustive/log1pf_test.cpp @@ -16,7 +16,7 @@ namespace mpfr = __llvm_libc::testing::mpfr; -struct LlvmLibclog1pfExhaustiveTest : public LlvmLibcExhaustiveTest { +struct LlvmLibcLog1pfExhaustiveTest : public LlvmLibcExhaustiveTest { bool check(uint32_t start, uint32_t stop, mpfr::RoundingMode rounding) override { mpfr::ForceRoundingMode r(rounding); @@ -35,22 +35,39 @@ // Range: All non-negative; static constexpr uint32_t START = 0x0000'0000U; static constexpr uint32_t STOP = 0x7f80'0000U; -// Range: [-1, 0]; -// static constexpr uint32_t START = 0x8000'0000U; -// static constexpr uint32_t STOP = 0xbf80'0000U; -TEST_F(LlvmLibclog1pfExhaustiveTest, RoundNearestTieToEven) { +TEST_F(LlvmLibcLog1pfExhaustiveTest, RoundNearestTieToEven) { test_full_range(START, STOP, mpfr::RoundingMode::Nearest); } -TEST_F(LlvmLibclog1pfExhaustiveTest, RoundUp) { +TEST_F(LlvmLibcLog1pfExhaustiveTest, RoundUp) { test_full_range(START, STOP, mpfr::RoundingMode::Upward); } -TEST_F(LlvmLibclog1pfExhaustiveTest, RoundDown) { +TEST_F(LlvmLibcLog1pfExhaustiveTest, RoundDown) { test_full_range(START, STOP, mpfr::RoundingMode::Downward); } -TEST_F(LlvmLibclog1pfExhaustiveTest, RoundTowardZero) { +TEST_F(LlvmLibcLog1pfExhaustiveTest, RoundTowardZero) { test_full_range(START, STOP, mpfr::RoundingMode::TowardZero); } + +// Range: [-1, 0]; +static constexpr uint32_t NEG_START = 0x8000'0000U; +static constexpr uint32_t NEG_STOP = 0xbf7f'ffffU; + +TEST_F(LlvmLibcLog1pfExhaustiveTest, NegativeRoundNearestTieToEven) { + test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Nearest); +} + +TEST_F(LlvmLibcLog1pfExhaustiveTest, NegativeRoundUp) { + test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Upward); +} + +TEST_F(LlvmLibcLog1pfExhaustiveTest, NegativeRoundDown) { + test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Downward); +} + +TEST_F(LlvmLibcLog1pfExhaustiveTest, NegativeRoundTowardZero) { + test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::TowardZero); +} diff --git a/libc/test/src/math/log10f_test.cpp b/libc/test/src/math/log10f_test.cpp --- a/libc/test/src/math/log10f_test.cpp +++ b/libc/test/src/math/log10f_test.cpp @@ -23,15 +23,16 @@ TEST(LlvmLibcLog10fTest, SpecialNumbers) { EXPECT_FP_EQ(aNaN, __llvm_libc::log10f(aNaN)); EXPECT_FP_EQ(inf, __llvm_libc::log10f(inf)); - EXPECT_TRUE(FPBits(__llvm_libc::log10f(neg_inf)).is_nan()); - EXPECT_FP_EQ(neg_inf, __llvm_libc::log10f(0.0f)); - EXPECT_FP_EQ(neg_inf, __llvm_libc::log10f(-0.0f)); - EXPECT_TRUE(FPBits(__llvm_libc::log10f(-1.0f)).is_nan()); + EXPECT_FP_IS_NAN_WITH_EXCEPTION(__llvm_libc::log10f(neg_inf), FE_INVALID); + EXPECT_FP_EQ_WITH_EXCEPTION(neg_inf, __llvm_libc::log10f(0.0f), FE_DIVBYZERO); + EXPECT_FP_EQ_WITH_EXCEPTION(neg_inf, __llvm_libc::log10f(-0.0f), + FE_DIVBYZERO); + EXPECT_FP_IS_NAN_WITH_EXCEPTION(__llvm_libc::log10f(-1.0f), FE_INVALID); EXPECT_FP_EQ(zero, __llvm_libc::log10f(1.0f)); } TEST(LlvmLibcLog10fTest, TrickyInputs) { - constexpr int N = 12; + constexpr int N = 15; constexpr uint32_t INPUTS[N] = { 0x41200000U /*10.0f*/, 0x42c80000U /*100.0f*/, @@ -44,7 +45,11 @@ 0x4e6e6b28U /*1,000,000,000.0f*/, 0x501502f9U /*10,000,000,000.0f*/, 0x4f134f83U /*2471461632.0f*/, - 0x7956ba5eU /*69683218960000541503257137270226944.0f*/}; + 0x7956ba5eU /*69683218960000541503257137270226944.0f*/, + 0x08ae'a356U /*0x1.5d46acp-110f*/, + 0x1c7d'a337U /*0x1.fb466ep-71f*/, + 0x69c8'c583U /*0x1.918b06p+84f*/, + }; for (int i = 0; i < N; ++i) { float x = float(FPBits(INPUTS[i])); @@ -58,16 +63,6 @@ constexpr uint32_t STEP = UINT32_MAX / COUNT; for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) { float x = float(FPBits(v)); - if (isnan(x) || isinf(x)) - continue; - errno = 0; - float result = __llvm_libc::log10f(x); - // If the computation resulted in an error or did not produce valid result - // in the single-precision floating point range, then ignore comparing with - // MPFR result as MPFR can still produce valid results because of its - // wider precision. - if (isnan(result) || isinf(result) || errno != 0) - continue; ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Log10, x, __llvm_libc::log10f(x), 0.5); } diff --git a/libc/test/src/math/log1pf_test.cpp b/libc/test/src/math/log1pf_test.cpp --- a/libc/test/src/math/log1pf_test.cpp +++ b/libc/test/src/math/log1pf_test.cpp @@ -23,14 +23,15 @@ TEST(LlvmLibclog1pfTest, SpecialNumbers) { EXPECT_FP_EQ(aNaN, __llvm_libc::log1pf(aNaN)); EXPECT_FP_EQ(inf, __llvm_libc::log1pf(inf)); - EXPECT_TRUE(FPBits((__llvm_libc::log1pf(neg_inf))).is_nan()); + EXPECT_FP_IS_NAN_WITH_EXCEPTION(__llvm_libc::log1pf(neg_inf), FE_INVALID); EXPECT_FP_EQ(zero, __llvm_libc::log1pf(0.0f)); EXPECT_FP_EQ(neg_zero, __llvm_libc::log1pf(-0.0f)); - EXPECT_FP_EQ(neg_inf, __llvm_libc::log1pf(-1.0f)); + EXPECT_FP_EQ_WITH_EXCEPTION(neg_inf, __llvm_libc::log1pf(-1.0f), + FE_DIVBYZERO); } TEST(LlvmLibclog1pfTest, TrickyInputs) { - constexpr int N = 20; + constexpr int N = 27; constexpr uint32_t INPUTS[N] = { 0x35c00006U, /*0x1.80000cp-20f*/ 0x35400003U, /*0x1.800006p-21f*/ @@ -41,10 +42,16 @@ 0x3770004bU, /*0x1.e00096p-17f*/ 0x3b9315c8U, /*0x1.262b9p-8f*/ 0x3c6eb7afU, /*0x1.dd6f5ep-7f*/ + 0x3ddbfec3U, /*0x1.b7fd86p-4f*/ + 0x3efd81adU, /*0x1.fb035ap-2f*/ 0x41078febU, /*0x1.0f1fd6p+3f*/ + 0x4cc1c80bU, /*0x1.839016p+26f*/ 0x5cd69e88U, /*0x1.ad3d1p+58f*/ + 0x5ee8984eU, /*0x1.d1309cp+62f*/ 0x65d890d3U, /*0x1.b121a6p+76f*/ + 0x665e7ca6U, /*0x1.bcf94cp+77f*/ 0x6f31a8ecU, /*0x1.6351d8p+95f*/ + 0x79e7ec37U, /*0x1.cfd86ep+116f*/ 0x7a17f30aU, /*0x1.2fe614p+117f*/ 0xb53ffffdU, /*-0x1.7ffffap-21f*/ 0xb70fffe5U, /*-0x1.1fffcap-17f*/ @@ -52,6 +59,7 @@ 0xbc4d092cU, /*-0x1.9a1258p-7f*/ 0xbc657728U, /*-0x1.caee5p-7f*/ 0xbd1d20afU, /*-0x1.3a415ep-5f*/ + 0xbf800000U, /*-1.0f*/ }; for (int i = 0; i < N; ++i) { float x = float(FPBits(INPUTS[i])); @@ -68,13 +76,6 @@ if (isnan(x) || isinf(x)) continue; errno = 0; - float result = __llvm_libc::log1pf(x); - // If the computation resulted in an error or did not produce valid result - // in the single-precision floating point range, then ignore comparing with - // MPFR result as MPFR can still produce valid results because of its - // wider precision. - if (isnan(result) || isinf(result) || errno != 0) - continue; ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Log1p, x, __llvm_libc::log1pf(x), 0.5); } diff --git a/libc/test/src/math/log2f_test.cpp b/libc/test/src/math/log2f_test.cpp --- a/libc/test/src/math/log2f_test.cpp +++ b/libc/test/src/math/log2f_test.cpp @@ -23,10 +23,10 @@ TEST(LlvmLibcLog2fTest, SpecialNumbers) { EXPECT_FP_EQ(aNaN, __llvm_libc::log2f(aNaN)); EXPECT_FP_EQ(inf, __llvm_libc::log2f(inf)); - EXPECT_TRUE(FPBits(__llvm_libc::log2f(neg_inf)).is_nan()); - EXPECT_FP_EQ(neg_inf, __llvm_libc::log2f(0.0f)); - EXPECT_FP_EQ(neg_inf, __llvm_libc::log2f(-0.0f)); - EXPECT_TRUE(FPBits(__llvm_libc::log2f(-1.0f)).is_nan()); + EXPECT_FP_IS_NAN_WITH_EXCEPTION(__llvm_libc::log2f(neg_inf), FE_INVALID); + EXPECT_FP_EQ_WITH_EXCEPTION(neg_inf, __llvm_libc::log2f(0.0f), FE_DIVBYZERO); + EXPECT_FP_EQ_WITH_EXCEPTION(neg_inf, __llvm_libc::log2f(-0.0f), FE_DIVBYZERO); + EXPECT_FP_IS_NAN_WITH_EXCEPTION(__llvm_libc::log2f(-1.0f), FE_INVALID); EXPECT_FP_EQ(zero, __llvm_libc::log2f(1.0f)); } diff --git a/libc/test/src/math/logf_test.cpp b/libc/test/src/math/logf_test.cpp --- a/libc/test/src/math/logf_test.cpp +++ b/libc/test/src/math/logf_test.cpp @@ -23,18 +23,22 @@ TEST(LlvmLibcLogfTest, SpecialNumbers) { EXPECT_FP_EQ(aNaN, __llvm_libc::logf(aNaN)); EXPECT_FP_EQ(inf, __llvm_libc::logf(inf)); - EXPECT_TRUE(FPBits((__llvm_libc::logf(neg_inf))).is_nan()); - EXPECT_FP_EQ(neg_inf, __llvm_libc::logf(0.0f)); - EXPECT_FP_EQ(neg_inf, __llvm_libc::logf(-0.0f)); - EXPECT_TRUE(FPBits(__llvm_libc::logf(-1.0f)).is_nan()); + EXPECT_FP_IS_NAN_WITH_EXCEPTION(__llvm_libc::logf(neg_inf), FE_INVALID); + EXPECT_FP_EQ_WITH_EXCEPTION(neg_inf, __llvm_libc::logf(0.0f), FE_DIVBYZERO); + EXPECT_FP_EQ_WITH_EXCEPTION(neg_inf, __llvm_libc::logf(-0.0f), FE_DIVBYZERO); + EXPECT_FP_IS_NAN_WITH_EXCEPTION(__llvm_libc::logf(-1.0f), FE_INVALID); EXPECT_FP_EQ(zero, __llvm_libc::logf(1.0f)); } TEST(LlvmLibcLogfTest, TrickyInputs) { - constexpr int N = 28; + constexpr int N = 35; constexpr uint32_t INPUTS[N] = { + 0x1b7679ffU, /*0x1.ecf3fep-73f*/ + 0x1e88452dU, /*0x1.108a5ap-66f*/ + 0x3f800001U, /*0x1.000002p+0f*/ 0x3509dcf6U, /*0x1.13b9ecp-21f*/ 0x3bf86ef0U, /*0x1.f0ddep-8f*/ + 0x3c413d3aU, /*0x1.827a74p-7f*/ 0x3ca1c99fU, /*0x1.43933ep-6f*/ 0x3d13e105U, /*0x1.27c20ap-5f*/ 0x3f7ff1f2U, /*0x1.ffe3e4p-1f*/ @@ -58,13 +62,17 @@ 0x4e85f412U, /*0x1.0be824p+30f*/ 0x500ffb03U, /*0x1.1ff606p+33f*/ 0x5cd69e88U, /*0x1.ad3d1p+58f*/ + 0x5ee8984eU, /*0x1.d1309cp+62f*/ 0x65d890d3U, /*0x1.b121a6p+76f*/ + 0x665e7ca6U, /*0x1.bcf94cp+77f*/ 0x6f31a8ecU, /*0x1.6351d8p+95f*/ + 0x79e7ec37U, /*0x1.cfd86ep+116f*/ 0x7a17f30aU, /*0x1.2fe614p+117f*/ }; for (int i = 0; i < N; ++i) { float x = float(FPBits(INPUTS[i])); - EXPECT_MPFR_MATCH(mpfr::Operation::Log, x, __llvm_libc::logf(x), 0.5); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Log, x, + __llvm_libc::logf(x), 0.5); } } @@ -75,14 +83,7 @@ float x = float(FPBits(v)); if (isnan(x) || isinf(x)) continue; - errno = 0; - float result = __llvm_libc::logf(x); - // If the computation resulted in an error or did not produce valid result - // in the single-precision floating point range, then ignore comparing with - // MPFR result as MPFR can still produce valid results because of its - // wider precision. - if (isnan(result) || isinf(result) || errno != 0) - continue; - ASSERT_MPFR_MATCH(mpfr::Operation::Log, x, __llvm_libc::logf(x), 0.5); + ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Log, x, + __llvm_libc::logf(x), 0.5); } } diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel --- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel +++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel @@ -845,7 +845,6 @@ name = "logf", additional_deps = [ ":__support_fputil_fma", - ":__support_fputil_nearest_integer", ":__support_fputil_multiply_add", ":__support_fputil_polyeval", ":common_constants", @@ -856,7 +855,6 @@ name = "log2f", additional_deps = [ ":__support_fputil_fma", - ":__support_fputil_nearest_integer", ":__support_fputil_multiply_add", ":__support_fputil_polyeval", ":common_constants", @@ -867,7 +865,6 @@ name = "log10f", additional_deps = [ ":__support_fputil_fma", - ":__support_fputil_nearest_integer", ":__support_fputil_multiply_add", ":__support_fputil_polyeval", ":common_constants", @@ -878,7 +875,6 @@ name = "log1pf", additional_deps = [ ":__support_fputil_fma", - ":__support_fputil_nearest_integer", ":__support_fputil_multiply_add", ":__support_fputil_polyeval", ":common_constants",