diff options
author | Tue Ly <lntue@google.com> | 2022-03-11 18:08:47 +0300 |
---|---|---|
committer | Tue Ly <lntue@google.com> | 2022-03-14 16:42:37 +0300 |
commit | 64721a3312c6e1b7189ffdcf05f9c6ee5b7e6777 (patch) | |
tree | 20d4192758f4e079291d2ae4ee3700ae089d9ce7 /libc | |
parent | 8cba72177dcd8de5d37177dbaf2347e5c1f0f1e8 (diff) |
[libc] Implement exp2f function that is correctly rounded for all rounding modes.
Implement exp2f function that is correctly rounded for all rounding modes.
Reviewed By: sivachandra, zimmermann6
Differential Revision: https://reviews.llvm.org/D121463
Diffstat (limited to 'libc')
-rw-r--r-- | libc/src/math/generic/CMakeLists.txt | 5 | ||||
-rw-r--r-- | libc/src/math/generic/exp2f.cpp | 183 | ||||
-rw-r--r-- | libc/test/src/math/exhaustive/CMakeLists.txt | 17 | ||||
-rw-r--r-- | libc/test/src/math/exhaustive/exp2f_test.cpp | 79 | ||||
-rw-r--r-- | libc/test/src/math/exp2f_test.cpp | 50 |
5 files changed, 278 insertions, 56 deletions
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index e7449bdf1753..e925e012f18e 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -489,9 +489,10 @@ add_entrypoint_object( HDRS ../exp2f.h DEPENDS - .exp_utils - .math_utils + libc.src.__support.FPUtil.fputil libc.include.math + COMPILE_OPTIONS + -O3 ) add_entrypoint_object( diff --git a/libc/src/math/generic/exp2f.cpp b/libc/src/math/generic/exp2f.cpp index f9cf8d81df20..b19cf8ebbd06 100644 --- a/libc/src/math/generic/exp2f.cpp +++ b/libc/src/math/generic/exp2f.cpp @@ -7,58 +7,159 @@ //===----------------------------------------------------------------------===// #include "src/math/exp2f.h" -#include "exp_utils.h" -#include "math_utils.h" - +#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/common.h" -#include <math.h> - -#include <stdint.h> -#define T exp2f_data.tab -#define C exp2f_data.poly -#define SHIFT exp2f_data.shift_scaled +#include <errno.h> namespace __llvm_libc { +// Lookup table for 2^(m * 2^(-7)) with m = 0, ..., 127. +// Table is generated with Sollya as follow: +// > display = hexadecimal; +// > for i from 0 to 127 do { D(2^(i / 128)); }; +static constexpr double EXP_M[128] = { + 0x1.0000000000000p0, 0x1.0163da9fb3335p0, 0x1.02c9a3e778061p0, + 0x1.04315e86e7f85p0, 0x1.059b0d3158574p0, 0x1.0706b29ddf6dep0, + 0x1.0874518759bc8p0, 0x1.09e3ecac6f383p0, 0x1.0b5586cf9890fp0, + 0x1.0cc922b7247f7p0, 0x1.0e3ec32d3d1a2p0, 0x1.0fb66affed31bp0, + 0x1.11301d0125b51p0, 0x1.12abdc06c31ccp0, 0x1.1429aaea92de0p0, + 0x1.15a98c8a58e51p0, 0x1.172b83c7d517bp0, 0x1.18af9388c8deap0, + 0x1.1a35beb6fcb75p0, 0x1.1bbe084045cd4p0, 0x1.1d4873168b9aap0, + 0x1.1ed5022fcd91dp0, 0x1.2063b88628cd6p0, 0x1.21f49917ddc96p0, + 0x1.2387a6e756238p0, 0x1.251ce4fb2a63fp0, 0x1.26b4565e27cddp0, + 0x1.284dfe1f56381p0, 0x1.29e9df51fdee1p0, 0x1.2b87fd0dad990p0, + 0x1.2d285a6e4030bp0, 0x1.2ecafa93e2f56p0, 0x1.306fe0a31b715p0, + 0x1.32170fc4cd831p0, 0x1.33c08b26416ffp0, 0x1.356c55f929ff1p0, + 0x1.371a7373aa9cbp0, 0x1.38cae6d05d866p0, 0x1.3a7db34e59ff7p0, + 0x1.3c32dc313a8e5p0, 0x1.3dea64c123422p0, 0x1.3fa4504ac801cp0, + 0x1.4160a21f72e2ap0, 0x1.431f5d950a897p0, 0x1.44e086061892dp0, + 0x1.46a41ed1d0057p0, 0x1.486a2b5c13cd0p0, 0x1.4a32af0d7d3dep0, + 0x1.4bfdad5362a27p0, 0x1.4dcb299fddd0dp0, 0x1.4f9b2769d2ca7p0, + 0x1.516daa2cf6642p0, 0x1.5342b569d4f82p0, 0x1.551a4ca5d920fp0, + 0x1.56f4736b527dap0, 0x1.58d12d497c7fdp0, 0x1.5ab07dd485429p0, + 0x1.5c9268a5946b7p0, 0x1.5e76f15ad2148p0, 0x1.605e1b976dc09p0, + 0x1.6247eb03a5585p0, 0x1.6434634ccc320p0, 0x1.6623882552225p0, + 0x1.68155d44ca973p0, 0x1.6a09e667f3bcdp0, 0x1.6c012750bdabfp0, + 0x1.6dfb23c651a2fp0, 0x1.6ff7df9519484p0, 0x1.71f75e8ec5f74p0, + 0x1.73f9a48a58174p0, 0x1.75feb564267c9p0, 0x1.780694fde5d3fp0, + 0x1.7a11473eb0187p0, 0x1.7c1ed0130c132p0, 0x1.7e2f336cf4e62p0, + 0x1.80427543e1a12p0, 0x1.82589994cce13p0, 0x1.8471a4623c7adp0, + 0x1.868d99b4492edp0, 0x1.88ac7d98a6699p0, 0x1.8ace5422aa0dbp0, + 0x1.8cf3216b5448cp0, 0x1.8f1ae99157736p0, 0x1.9145b0b91ffc6p0, + 0x1.93737b0cdc5e5p0, 0x1.95a44cbc8520fp0, 0x1.97d829fde4e50p0, + 0x1.9a0f170ca07bap0, 0x1.9c49182a3f090p0, 0x1.9e86319e32323p0, + 0x1.a0c667b5de565p0, 0x1.a309bec4a2d33p0, 0x1.a5503b23e255dp0, + 0x1.a799e1330b358p0, 0x1.a9e6b5579fdbfp0, 0x1.ac36bbfd3f37ap0, + 0x1.ae89f995ad3adp0, 0x1.b0e07298db666p0, 0x1.b33a2b84f15fbp0, + 0x1.b59728de5593ap0, 0x1.b7f76f2fb5e47p0, 0x1.ba5b030a1064ap0, + 0x1.bcc1e904bc1d2p0, 0x1.bf2c25bd71e09p0, 0x1.c199bdd85529cp0, + 0x1.c40ab5fffd07ap0, 0x1.c67f12e57d14bp0, 0x1.c8f6d9406e7b5p0, + 0x1.cb720dcef9069p0, 0x1.cdf0b555dc3fap0, 0x1.d072d4a07897cp0, + 0x1.d2f87080d89f2p0, 0x1.d5818dcfba487p0, 0x1.d80e316c98398p0, + 0x1.da9e603db3285p0, 0x1.dd321f301b460p0, 0x1.dfc97337b9b5fp0, + 0x1.e264614f5a129p0, 0x1.e502ee78b3ff6p0, 0x1.e7a51fbc74c83p0, + 0x1.ea4afa2a490dap0, 0x1.ecf482d8e67f1p0, 0x1.efa1bee615a27p0, + 0x1.f252b376bba97p0, 0x1.f50765b6e4540p0, 0x1.f7bfdad9cbe14p0, + 0x1.fa7c1819e90d8p0, 0x1.fd3c22b8f71f1p0, +}; + +INLINE_FMA LLVM_LIBC_FUNCTION(float, exp2f, (float x)) { - uint32_t abstop; - uint64_t ki, t; - // double_t for better performance on targets with FLT_EVAL_METHOD==2. - double_t kd, xd, z, r, r2, y, s; + using FPBits = typename fputil::FPBits<float>; + FPBits xbits(x); - xd = static_cast<double_t>(x); - abstop = top12_bits(x) & 0x7ff; - if (unlikely(abstop >= top12_bits(128.0f))) { - // |x| >= 128 or x is nan. - if (as_uint32_bits(x) == as_uint32_bits(-INFINITY)) + // When x < -150 or nan + if (unlikely(xbits.uintval() > 0xc316'0000U)) { + // exp(-Inf) = 0 + if (xbits.is_inf()) return 0.0f; - if (abstop >= top12_bits(INFINITY)) - return x + x; - if (x > 0.0f) - return overflow<float>(0); - if (x <= -150.0f) - return underflow<float>(0); - if (x < -149.0f) - return may_underflow<float>(0); + // exp(nan) = nan + if (xbits.is_nan()) + return x; + if (fputil::get_round() == FE_UPWARD) + return static_cast<float>(FPBits(FPBits::MIN_SUBNORMAL)); + if (x != 0.0f) + errno = ERANGE; + return 0.0f; } + // x >= 128 or nan + if (unlikely(!xbits.get_sign() && (xbits.uintval() >= 0x4300'0000U))) { + if (xbits.uintval() < 0x7f80'0000U) { + int rounding = fputil::get_round(); + if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO) + return static_cast<float>(FPBits(FPBits::MAX_NORMAL)); - // x = k/N + r with r in [-1/(2N), 1/(2N)] and int k. - kd = static_cast<double>(xd + SHIFT); - ki = as_uint64_bits(kd); - kd -= SHIFT; // k/N for int k. - r = xd - kd; + errno = ERANGE; + } + return x + static_cast<float>(FPBits::inf()); + } + // |x| < 2^-25 + if (unlikely(xbits.get_unbiased_exponent() <= 101)) { + return 1.0f + x; + } + // Exceptional values. + switch (xbits.uintval()) { + case 0x3b42'9d37U: // x = 0x1.853a6ep-9f + if (fputil::get_round() == FE_TONEAREST) + return 0x1.00870ap+0f; + break; + case 0xbcf3'a937U: // x = -0x1.e7526ep-6f + if (fputil::get_round() == FE_TONEAREST) + return 0x1.f58d62p-1f; + break; + } + + // For -150 <= x < 128, to compute 2^x, we perform the following range + // reduction: find hi, mid, lo such that: + // x = hi + mid + lo, in which + // hi is an integer, + // mid * 2^7 is an integer + // -2^(-8) <= lo < 2^-8. + // In particular, + // hi + mid = round(x * 2^7) * 2^(-7). + // Then, + // 2^(x) = 2^(hi + mid + lo) = 2^hi * 2^mid * 2^lo. + // Multiply by 2^hi is simply adding hi to the exponent field. We store + // exp(mid) in the lookup tables EXP_M. exp(lo) is computed using a degree-7 + // minimax polynomial generated by Sollya. - // exp2(x) = 2^(k/N) * 2^r ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) - t = T[ki % N]; - t += ki << (52 - EXP2F_TABLE_BITS); - s = as_double(t); - z = C[0] * r + C[1]; - r2 = r * r; - y = C[2] * r + 1; - y = z * r2 + y; - y = y * s; - return static_cast<float>(y); + // x_hi = hi + mid. + int x_hi = static_cast<int>(x * 0x1.0p7f); + // Subtract (hi + mid) from x to get lo. + x -= static_cast<float>(x_hi) * 0x1.0p-7f; + double xd = static_cast<double>(x); + // Make sure that -2^(-8) <= lo < 2^-8. + if (x >= 0x1.0p-8f) { + ++x_hi; + xd -= 0x1.0p-7; + } + if (x < -0x1.0p-8f) { + --x_hi; + xd += 0x1.0p-7; + } + // For 2-complement integers, arithmetic right shift is the same as dividing + // by a power of 2 and then round down (toward negative infinity). + int hi = x_hi >> 7; + // mid = x_hi & 0x0000'007fU; + double exp_mid = EXP_M[x_hi & 0x7f]; + // Degree-6 minimax polynomial generated by Sollya with the following + // commands: + // > display = hexadecimal; + // > Q = fpminimax((2^x - 1)/x, 5, [|D...|], [-2^-8, 2^-8]); + // > Q; + double exp_lo = + fputil::polyeval(xd, 0x1p0, 0x1.62e42fefa39efp-1, 0x1.ebfbdff82c58ep-3, + 0x1.c6b08d711fe2fp-5, 0x1.3b2ab6fe3deb5p-7, + 0x1.5d72a05f45c04p-10, 0x1.4284d40c33326p-13); + fputil::FPBits<double> result(exp_mid * exp_lo); + result.set_unbiased_exponent(static_cast<uint16_t>( + static_cast<int>(result.get_unbiased_exponent()) + hi)); + return static_cast<float>(static_cast<double>(result)); } } // namespace __llvm_libc diff --git a/libc/test/src/math/exhaustive/CMakeLists.txt b/libc/test/src/math/exhaustive/CMakeLists.txt index 01ffc79dde21..c4ee2bb12315 100644 --- a/libc/test/src/math/exhaustive/CMakeLists.txt +++ b/libc/test/src/math/exhaustive/CMakeLists.txt @@ -65,6 +65,23 @@ add_fp_unittest( ) add_fp_unittest( + exp2f_test + NO_RUN_POSTBUILD + NEED_MPFR + SUITE + libc_math_exhaustive_tests + SRCS + exp2f_test.cpp + DEPENDS + .exhaustive_test + libc.include.math + libc.src.math.exp2f + libc.src.__support.FPUtil.fputil + LINK_OPTIONS + -lpthread +) + +add_fp_unittest( expm1f_test NEED_MPFR SUITE diff --git a/libc/test/src/math/exhaustive/exp2f_test.cpp b/libc/test/src/math/exhaustive/exp2f_test.cpp new file mode 100644 index 000000000000..48fa434ff332 --- /dev/null +++ b/libc/test/src/math/exhaustive/exp2f_test.cpp @@ -0,0 +1,79 @@ +//===-- Exhaustive test for exp2f -----------------------------------------===// +// +// 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 "exhaustive_test.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/math/exp2f.h" +#include "utils/MPFRWrapper/MPFRUtils.h" +#include "utils/UnitTest/FPMatcher.h" + +using FPBits = __llvm_libc::fputil::FPBits<float>; + +namespace mpfr = __llvm_libc::testing::mpfr; + +struct LlvmLibcExp2fExhaustiveTest : public LlvmLibcExhaustiveTest<uint32_t> { + void check(uint32_t start, uint32_t stop, + mpfr::RoundingMode rounding) override { + mpfr::ForceRoundingMode r(rounding); + uint32_t bits = start; + do { + FPBits xbits(bits); + float x = float(xbits); + EXPECT_MPFR_MATCH(mpfr::Operation::Exp2, x, __llvm_libc::exp2f(x), 0.5, + rounding); + } while (bits++ < stop); + } +}; + +static constexpr int NUM_THREADS = 16; + +// Range: [0, 128]; +static constexpr uint32_t POS_START = 0x0000'0000U; +static constexpr uint32_t POS_STOP = 0x4300'0000U; + +TEST_F(LlvmLibcExp2fExhaustiveTest, PostiveRangeRoundNearestTieToEven) { + test_full_range(POS_START, POS_STOP, NUM_THREADS, + mpfr::RoundingMode::Nearest); +} + +TEST_F(LlvmLibcExp2fExhaustiveTest, PostiveRangeRoundUp) { + test_full_range(POS_START, POS_STOP, NUM_THREADS, mpfr::RoundingMode::Upward); +} + +TEST_F(LlvmLibcExp2fExhaustiveTest, PostiveRangeRoundDown) { + test_full_range(POS_START, POS_STOP, NUM_THREADS, + mpfr::RoundingMode::Downward); +} + +TEST_F(LlvmLibcExp2fExhaustiveTest, PostiveRangeRoundTowardZero) { + test_full_range(POS_START, POS_STOP, NUM_THREADS, + mpfr::RoundingMode::TowardZero); +} + +// Range: [-150, 0]; +static constexpr uint32_t NEG_START = 0x8000'0000U; +static constexpr uint32_t NEG_STOP = 0xc316'0000U; + +TEST_F(LlvmLibcExp2fExhaustiveTest, NegativeRangeRoundNearestTieToEven) { + test_full_range(NEG_START, NEG_STOP, NUM_THREADS, + mpfr::RoundingMode::Nearest); +} + +TEST_F(LlvmLibcExp2fExhaustiveTest, NegativeRangeRoundUp) { + test_full_range(NEG_START, NEG_STOP, NUM_THREADS, mpfr::RoundingMode::Upward); +} + +TEST_F(LlvmLibcExp2fExhaustiveTest, NegativeRangeRoundDown) { + test_full_range(NEG_START, NEG_STOP, NUM_THREADS, + mpfr::RoundingMode::Downward); +} + +TEST_F(LlvmLibcExp2fExhaustiveTest, NegativeRangeRoundTowardZero) { + test_full_range(NEG_START, NEG_STOP, NUM_THREADS, + mpfr::RoundingMode::TowardZero); +} diff --git a/libc/test/src/math/exp2f_test.cpp b/libc/test/src/math/exp2f_test.cpp index 4cc73830f7d5..20c22c6bbb9e 100644 --- a/libc/test/src/math/exp2f_test.cpp +++ b/libc/test/src/math/exp2f_test.cpp @@ -39,7 +39,7 @@ TEST(LlvmLibcExp2fTest, SpecialNumbers) { EXPECT_MATH_ERRNO(0); } -TEST(LlvmLibcExpfTest, Overflow) { +TEST(LlvmLibcExp2fTest, Overflow) { errno = 0; EXPECT_FP_EQ(inf, __llvm_libc::exp2f(float(FPBits(0x7f7fffffU)))); EXPECT_MATH_ERRNO(ERANGE); @@ -53,50 +53,73 @@ TEST(LlvmLibcExpfTest, Overflow) { // Test with inputs which are the borders of underflow/overflow but still // produce valid results without setting errno. -TEST(LlvmLibcExpfTest, Borderline) { +TEST(LlvmLibcExp2fTest, Borderline) { float x; errno = 0; x = float(FPBits(0x42fa0001U)); - EXPECT_MPFR_MATCH(mpfr::Operation::Exp2, x, __llvm_libc::exp2f(x), 1.0); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x, + __llvm_libc::exp2f(x), 0.5); EXPECT_MATH_ERRNO(0); x = float(FPBits(0x42ffffffU)); - EXPECT_MPFR_MATCH(mpfr::Operation::Exp2, x, __llvm_libc::exp2f(x), 1.0); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x, + __llvm_libc::exp2f(x), 0.5); EXPECT_MATH_ERRNO(0); x = float(FPBits(0xc2fa0001U)); - EXPECT_MPFR_MATCH(mpfr::Operation::Exp2, x, __llvm_libc::exp2f(x), 1.0); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x, + __llvm_libc::exp2f(x), 0.5); EXPECT_MATH_ERRNO(0); x = float(FPBits(0xc2fc0000U)); - EXPECT_MPFR_MATCH(mpfr::Operation::Exp2, x, __llvm_libc::exp2f(x), 1.0); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x, + __llvm_libc::exp2f(x), 0.5); EXPECT_MATH_ERRNO(0); x = float(FPBits(0xc2fc0001U)); - EXPECT_MPFR_MATCH(mpfr::Operation::Exp2, x, __llvm_libc::exp2f(x), 1.0); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x, + __llvm_libc::exp2f(x), 0.5); EXPECT_MATH_ERRNO(0); x = float(FPBits(0xc3150000U)); - EXPECT_MPFR_MATCH(mpfr::Operation::Exp2, x, __llvm_libc::exp2f(x), 1.0); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x, + __llvm_libc::exp2f(x), 0.5); + EXPECT_MATH_ERRNO(0); + + x = float(FPBits(0x3b42'9d37U)); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x, + __llvm_libc::exp2f(x), 0.5); + EXPECT_MATH_ERRNO(0); + + x = float(FPBits(0xbcf3'a937U)); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x, + __llvm_libc::exp2f(x), 0.5); EXPECT_MATH_ERRNO(0); } -TEST(LlvmLibcExpfTest, Underflow) { +TEST(LlvmLibcExp2fTest, Underflow) { errno = 0; EXPECT_FP_EQ(0.0f, __llvm_libc::exp2f(float(FPBits(0xff7fffffU)))); EXPECT_MATH_ERRNO(ERANGE); float x = float(FPBits(0xc3158000U)); - EXPECT_MPFR_MATCH(mpfr::Operation::Exp2, x, __llvm_libc::exp2f(x), 1.0); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x, + __llvm_libc::exp2f(x), 0.5); + EXPECT_MATH_ERRNO(0); + + x = float(FPBits(0xc3160000U)); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x, + __llvm_libc::exp2f(x), 0.5); EXPECT_MATH_ERRNO(ERANGE); x = float(FPBits(0xc3165432U)); - EXPECT_MPFR_MATCH(mpfr::Operation::Exp2, x, __llvm_libc::exp2f(x), 1.0); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x, + __llvm_libc::exp2f(x), 0.5); EXPECT_MATH_ERRNO(ERANGE); } -TEST(LlvmLibcexp2fTest, InFloatRange) { +TEST(LlvmLibcExp2fTest, InFloatRange) { constexpr uint32_t COUNT = 1000000; constexpr uint32_t STEP = UINT32_MAX / COUNT; for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) { @@ -112,6 +135,7 @@ TEST(LlvmLibcexp2fTest, InFloatRange) { // wider precision. if (isnan(result) || isinf(result) || errno != 0) continue; - ASSERT_MPFR_MATCH(mpfr::Operation::Exp2, x, __llvm_libc::exp2f(x), 1.0); + ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x, + __llvm_libc::exp2f(x), 0.5); } } |