diff --git a/libc/src/math/generic/expf.cpp b/libc/src/math/generic/expf.cpp --- a/libc/src/math/generic/expf.cpp +++ b/libc/src/math/generic/expf.cpp @@ -8,6 +8,7 @@ #include "src/math/expf.h" #include "common_constants.h" // Lookup tables EXP_M1 and EXP_M2. +#include "expxf.h" #include "src/__support/FPUtil/BasicOperations.h" #include "src/__support/FPUtil/FEnvImpl.h" #include "src/__support/FPUtil/FPBits.h" @@ -27,11 +28,6 @@ uint32_t x_u = xbits.uintval(); uint32_t x_abs = x_u & 0x7fff'ffffU; - // Exceptional values - if (unlikely(x_u == 0xc236'bd8cU)) { // x = -0x1.6d7b18p+5f - return 0x1.108a58p-66f - x * 0x1.0p-95f; - } - // When |x| >= 89, |x| < 2^-25, or x is nan if (unlikely(x_abs >= 0x42b2'0000U || x_abs <= 0x3280'0000U)) { // |x| < 2^-25 @@ -42,13 +38,15 @@ // When x < log(2^-150) or nan if (xbits.uintval() >= 0xc2cf'f1b5U) { // exp(-Inf) = 0 - if (xbits.is_inf()) + if (xbits.is_inf()) { return 0.0f; + } // exp(nan) = nan if (xbits.is_nan()) return x; - if (fputil::get_round() == FE_UPWARD) - return static_cast(FPBits(FPBits::MIN_SUBNORMAL)); + + if (unlikely(fputil::get_round() == FE_UPWARD)) + return FPBits(FPBits::MIN_SUBNORMAL).get_val(); errno = ERANGE; return 0.0f; } @@ -57,48 +55,18 @@ // x is finite if (xbits.uintval() < 0x7f80'0000U) { int rounding = fputil::get_round(); - if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO) - return static_cast(FPBits(FPBits::MAX_NORMAL)); + if (unlikely(rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)) + return FPBits(FPBits::MAX_NORMAL).get_val(); errno = ERANGE; } // x is +inf or nan - return x + static_cast(FPBits::inf()); + return x + FPBits::inf().get_val(); } } - // For -104 < x < 89, to compute exp(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, - // exp(x) = exp(hi + mid + lo) = exp(hi) * exp(mid) * exp(lo). - // We store exp(hi) and exp(mid) in the lookup tables EXP_M1 and EXP_M2 - // respectively. exp(lo) is computed using a degree-4 minimax polynomial - // generated by Sollya. - // x_hi = (hi + mid) * 2^7 = round(x * 2^7). - float kf = fputil::nearest_integer(x * 0x1.0p7f); - // Subtract (hi + mid) from x to get lo. - double xd = static_cast(fputil::multiply_add(kf, -0x1.0p-7f, x)); - int x_hi = static_cast(kf); - x_hi += 104 << 7; - // hi = x_hi >> 7 - double exp_hi = EXP_M1[x_hi >> 7]; - // mid * 2^7 = x_hi & 0x0000'007fU; - double exp_mid = EXP_M2[x_hi & 0x7f]; - // Degree-4 minimax polynomial generated by Sollya with the following - // commands: - // > display = hexadecimal; - // > Q = fpminimax(expm1(x)/x, 3, [|D...|], [-2^-8, 2^-8]); - // > Q; - double exp_lo = - fputil::polyeval(xd, 0x1p0, 0x1.ffffffffff777p-1, 0x1.000000000071cp-1, - 0x1.555566668e5e7p-3, 0x1.55555555ef243p-5); - return static_cast(exp_hi * exp_mid * exp_lo); + auto ep = __llvm_libc::exp_eval(x); + return fputil::multiply_add(ep.mult_exp, ep.r, ep.mult_exp); } } // namespace __llvm_libc diff --git a/libc/src/math/generic/expm1f.cpp b/libc/src/math/generic/expm1f.cpp --- a/libc/src/math/generic/expm1f.cpp +++ b/libc/src/math/generic/expm1f.cpp @@ -8,6 +8,7 @@ #include "src/math/expm1f.h" #include "common_constants.h" // Lookup tables EXP_M1 and EXP_M2. +#include "expxf.h" #include "src/__support/FPUtil/BasicOperations.h" #include "src/__support/FPUtil/FEnvImpl.h" #include "src/__support/FPUtil/FMA.h" @@ -28,132 +29,45 @@ uint32_t x_u = xbits.uintval(); uint32_t x_abs = x_u & 0x7fff'ffffU; - // Exceptional value - if (unlikely(x_u == 0x3e35'bec5U)) { // x = 0x1.6b7d8ap-3f - int round_mode = fputil::get_round(); - if (round_mode == FE_TONEAREST || round_mode == FE_UPWARD) - return 0x1.8dbe64p-3f; - return 0x1.8dbe62p-3f; - } - -#if !defined(LIBC_TARGET_HAS_FMA) - if (unlikely(x_u == 0xbdc1'c6cbU)) { // x = -0x1.838d96p-4f - int round_mode = fputil::get_round(); - if (round_mode == FE_TONEAREST || round_mode == FE_DOWNWARD) - return -0x1.71c884p-4f; - return -0x1.71c882p-4f; - } -#endif // LIBC_TARGET_HAS_FMA + // When |x| >= 89, |x| < 2^-25, or x is nan or x < -18 + if (unlikely(x_abs >= 0x42b2'0000U || x_abs <= 0x3280'0000U || + x_u >= 0xc190'0000U)) { + // |x| < 2^-25 + if (xbits.get_unbiased_exponent() <= 101) { + return unlikely(x_abs == 0) ? x : (x + 0.5 * x * x); + } - // When |x| > 25*log(2), or nan - if (unlikely(x_abs >= 0x418a'a123U)) { - // x < log(2^-25) - if (xbits.get_sign()) { + // When x < log(2^-150) or nan + if (xbits.uintval() >= 0xc180'0000U) { // exp(-Inf) = 0 - if (xbits.is_inf()) + if (xbits.is_inf()) { return -1.0f; + } + // exp(nan) = nan if (xbits.is_nan()) return x; - int round_mode = fputil::get_round(); - if (round_mode == FE_UPWARD || round_mode == FE_TOWARDZERO) - return -0x1.ffff'fep-1f; // -1.0f + 0x1.0p-24f - return -1.0f; - } else { - // x >= 89 or nan - if (xbits.uintval() >= 0x42b2'0000) { - if (xbits.uintval() < 0x7f80'0000U) { - int rounding = fputil::get_round(); - if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO) - return static_cast(FPBits(FPBits::MAX_NORMAL)); - errno = ERANGE; - } - return x + static_cast(FPBits::inf()); - } + return -1.0f + opt_barrier(FPBits(FPBits::MIN_NORMAL).get_val()); } - } - // |x| < 2^-4 - if (x_abs < 0x3d80'0000U) { - // |x| < 2^-25 - if (x_abs < 0x3300'0000U) { - // x = -0.0f - if (unlikely(xbits.uintval() == 0x8000'0000U)) - return x; - // When |x| < 2^-25, the relative error of the approximation e^x - 1 ~ x - // is: - // |(e^x - 1) - x| / |e^x - 1| < |x^2| / |x| - // = |x| - // < 2^-25 - // < epsilon(1)/2. - // So the correctly rounded values of expm1(x) are: - // = x + eps(x) if rounding mode = FE_UPWARD, - // or (rounding mode = FE_TOWARDZERO and x is - // negative), - // = x otherwise. - // To simplify the rounding decision and make it more efficient, we use - // fma(x, x, x) ~ x + x^2 instead. - // Note: to use the formula x + x^2 to decide the correct rounding, we - // do need fma(x, x, x) to prevent underflow caused by x*x when |x| < - // 2^-76. For targets without FMA instructions, we simply use double for - // intermediate results as it is more efficient than using an emulated - // version of FMA. -#if defined(LIBC_TARGET_HAS_FMA) - return fputil::fma(x, x, x); -#else - double xd = x; - return static_cast(fputil::multiply_add(xd, xd, xd)); -#endif // LIBC_TARGET_HAS_FMA - } + // x >= 89 or nan + if (!xbits.get_sign() && (xbits.uintval() >= 0x42b2'0000)) { + // x is finite + if (xbits.uintval() < 0x7f80'0000U) { + int rounding = fputil::get_round(); + if (unlikely(rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)) + return FPBits(FPBits::MAX_NORMAL).get_val(); - // 2^-25 <= |x| < 2^-4 - double xd = static_cast(x); - double xsq = xd * xd; - // Degree-8 minimax polynomial generated by Sollya with: - // > display = hexadecimal; - // > P = fpminimax((expm1(x) - x)/x^2, 6, [|D...|], [-2^-4, 2^-4]); - double r = - fputil::polyeval(xd, 0x1p-1, 0x1.55555555557ddp-3, 0x1.55555555552fap-5, - 0x1.111110fcd58b7p-7, 0x1.6c16c1717660bp-10, - 0x1.a0241f0006d62p-13, 0x1.a01e3f8d3c06p-16); - return static_cast(fputil::multiply_add(r, xsq, xd)); + errno = ERANGE; + } + // x is +inf or nan + return x + FPBits::inf().get_val(); + } } - // For -18 < x < 89, to compute expm1(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, - // expm1(x) = exp(hi + mid + lo) - 1 = exp(hi) * exp(mid) * exp(lo) - 1. - // We store exp(hi) and exp(mid) in the lookup tables EXP_M1 and EXP_M2 - // respectively. exp(lo) is computed using a degree-4 minimax polynomial - // generated by Sollya. - - // x_hi = hi + mid. - float kf = fputil::nearest_integer(x * 0x1.0p7f); - int x_hi = static_cast(kf); - // Subtract (hi + mid) from x to get lo. - double xd = static_cast(fputil::multiply_add(kf, -0x1.0p-7f, x)); - x_hi += 104 << 7; - // hi = x_hi >> 7 - double exp_hi = EXP_M1[x_hi >> 7]; - // lo = x_hi & 0x0000'007fU; - double exp_mid = EXP_M2[x_hi & 0x7f]; - double exp_hi_mid = exp_hi * exp_mid; - // Degree-4 minimax polynomial generated by Sollya with the following - // commands: - // > display = hexadecimal; - // > Q = fpminimax(expm1(x)/x, 3, [|D...|], [-2^-8, 2^-8]); - // > Q; - double exp_lo = - fputil::polyeval(xd, 0x1.0p0, 0x1.ffffffffff777p-1, 0x1.000000000071cp-1, - 0x1.555566668e5e7p-3, 0x1.55555555ef243p-5); - return static_cast(fputil::multiply_add(exp_hi_mid, exp_lo, -1.0)); + auto ep = __llvm_libc::exp_eval(x); + return fputil::multiply_add(ep.mult_exp, ep.r, ep.mult_exp - 1.0); } } // namespace __llvm_libc