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 @@ -561,6 +561,7 @@ ../exp2f.h DEPENDS .explogxf + libc.src.__support.FPUtil.except_value_utils libc.src.__support.FPUtil.fenv_impl libc.src.__support.FPUtil.fp_bits libc.src.__support.FPUtil.multiply_add diff --git a/libc/src/math/generic/exp2f.cpp b/libc/src/math/generic/exp2f.cpp --- a/libc/src/math/generic/exp2f.cpp +++ b/libc/src/math/generic/exp2f.cpp @@ -10,6 +10,7 @@ #include "src/__support/FPUtil/FEnvImpl.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/FPUtil/nearest_integer.h" #include "src/__support/FPUtil/rounding_mode.h" @@ -33,13 +34,38 @@ uint32_t x_u = xbits.uintval(); uint32_t x_abs = x_u & 0x7fff'ffffU; - // |x| < 2^-25 - if (LIBC_UNLIKELY(x_abs <= 0x3280'0000U)) { - return 1.0f + x; - } + // When |x| >= 128, or x is nan, or |x| <= 2^-5 + if (LIBC_UNLIKELY(x_abs >= 0x4300'0000U || x_abs <= 0x3d00'0000U)) { + // |x| <= 2^-5 + if (x_abs <= 0x3d00'0000) { + // |x| < 2^-25 + if (LIBC_UNLIKELY(x_abs <= 0x3280'0000U)) { + return 1.0f + x; + } + + // Check exceptional values. + if (LIBC_UNLIKELY((x_u & EXVAL_MASK) == EXVAL_MASK)) { + if (LIBC_UNLIKELY(x_u == EXVAL1)) { // x = 0x1.853a6ep-9f + return fputil::round_result_slightly_down(0x1.00870ap+0f); + } else if (LIBC_UNLIKELY(x_u == EXVAL2)) { // x = -0x1.e7526ep-6f + return fputil::round_result_slightly_down(0x1.f58d62p-1f); + } + } - // // When |x| >= 128, or x is nan - if (LIBC_UNLIKELY(x_abs >= 0x4300'0000U)) { + // Minimax polynomial generated by Sollya with: + // > P = fpminimax((2^x - 1)/x, 5, [|D...|], [-2^-5, 2^-5]); + constexpr double COEFFS[] = { + 0x1.62e42fefa39f3p-1, 0x1.ebfbdff82c57bp-3, 0x1.c6b08d6f2d7aap-5, + 0x1.3b2ab6fc92f5dp-7, 0x1.5d897cfe27125p-10, 0x1.43090e61e6af1p-13}; + double xd = static_cast(x); + 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 p = fputil::polyeval(xsq, c0, c1, c2); + double r = fputil::multiply_add(p, xd, 1.0); + return r; + } // x >= 128 if (!xbits.get_sign()) { @@ -73,17 +99,6 @@ } } - // Check exceptional values. - if (LIBC_UNLIKELY((x_u & EXVAL_MASK) == EXVAL_MASK)) { - if (LIBC_UNLIKELY(x_u == EXVAL1)) { // x = 0x1.853a6ep-9f - if (fputil::fenv_is_round_to_nearest()) - return 0x1.00870ap+0f; - } else if (LIBC_UNLIKELY(x_u == EXVAL2)) { // x = -0x1.e7526ep-6f - if (fputil::fenv_is_round_to_nearest()) - return 0x1.f58d62p-1f; - } - } - // 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