[libc-commits] [libc] [libc][math] Make pow function correctly rounded for all rounding modes. (PR #222827)
via libc-commits
libc-commits at lists.llvm.org
Thu Sep 17 10:56:33 PDT 2026
================
@@ -451,96 +275,230 @@ LIBC_INLINE double pow(double x, double y) {
// and 2^lo ~ 1 + lo * P(lo).
// Thus, we have:
// hi + mid = 2^-6 * round( 2^6 * y * log2(x) )
- // If we restrict the output such that |hi| < 150, (hi + mid) uses (8 + 6)
+ // If we restrict the output such that |hi| < 512, (hi + mid) uses (9 + 6)
// bits, hence, if we use double precision to perform
// round( 2^6 * y * log2(x))
- // the lo part is bounded by 2^-7 + 2^(-(52 - 14)) = 2^-7 + 2^-38
+ // the lo part is bounded by 2^-7 + 2^(-(52 - 15)) = 2^-7 + 2^-37
// In the following computations:
// y6 = 2^6 * y
// hm = 2^6 * (hi + mid) = round(2^6 * y * log2(x)) ~ round(y6 * s)
// lo6 = 2^6 * lo = 2^6 * (y - (hi + mid)) = y6 * log2(x) - hm.
- double y6 = y * 0x1.0p6; // Exact.
+ constexpr double SCALE = 0x1.0p6;
+ double y6 = y * SCALE; // Exact.
DoubleDouble y6_log2_x = fputil::exact_mult(y6, log2_x.hi);
y6_log2_x.lo = fputil::multiply_add(y6, log2_x.lo, y6_log2_x.lo);
// Check overflow/underflow.
double scale = 1.0;
+ bool is_denorm = false;
// |2^(hi + mid) - exp2_hi_mid| <= ulp(exp2_hi_mid) / 2
- // Clamp the exponent part into smaller range that fits double precision.
- // For those exponents that are out of range, the final conversion will round
- // them correctly to inf/max float or 0/min float accordingly.
- constexpr double UPPER_EXP_BOUND = 512.0 * 0x1.0p6;
+
+ // The fast computation for 2^hi below requires that:
+ // |hi| < 512, or equivalently, |hm| < 512 * 2^6.
+ // This guarantees that the biased exponent:
+ // exp_biased = (hm_i >> 6) + EXP_BIAS = hi + 1023
+ // is strictly within the normal double range:
+ // 1023 - 511 <= exp_biased <= 1023 + 511, or 512 <= exp_biased <= 1534.
+ // Hence, 2^hi is always a normal, non-zero, finite power of 2, and the
+ // multiplication upper * exp2_hi is exact and never underflows or overflows.
+ //
+ // From the edge case checks above:
+ // 2^(-54) / 1074 <= |y| <= 1075 * 2^53, and
+ // |log_2(1 - 2^(-53))| <= |log2(x)| <= 1074.
+ // So their product is bounded by:
+ // 2^-117 < |y * log2(x)| < 2^74.
+ //
+ // The meaningful range for y * log2(x) in double precision is:
+ // -1075 <= y * log2(x) <= 1024.
+ // Any value > 1024 overflows, and any value < -1075 underflows.
+ //
+ // When |y * log2(x)| >= 511, we shift the exponent by an offset S:
+ // y * log2(x) = (y * log2(x) - S) + S
+ // and multiply by scale = 2^S at the end.
+ // To ensure the shifted exponent (y * log2(x) - S) stays within [-511, 511]:
+ // - For positive range [511, 1024]:
+ // 511 - S > -511 ==> S < 1022
+ // 1024 - S < 511 ==> S > 513
+ // - For negative range [-1076, -511]:
+ // -511 + S < 511 ==> S <= 1022
+ // -1076 + S > -511 ==> S >= 565
+ // Combined, the shift must satisfy: 565 <= S <= 1022.
+ // We choose S = 600, which yields:
+ // y * log2(x) - 600 in [-89, 424] for the positive range, and
+ // y * log2(x) + 600 in [-476, 89] for the negative range,
+ // both fitting well within [-511, 511].
+ //
+ // For exponents that are completely out of range:
+ // y * log2(x) > 1025 or y * log2(x) < -1076,
+ // the product can be as large as 2^74, so subtracting or adding 600 * 64
+ // would still overflow a 32-bit int when computing hm_i.
+ // We return early with overflow or underflow in these cases.
+ //
+ // Alternatively, for less branching or in SIMD / vector implementations, one
+ // could clamp y6_log2_x.hi to:
+ // - UPPER_EXP_BOUND (511 * 64) for overflow, which when multiplied by
+ // scale = 2^600 yields 2^1111 and correctly overflows.
+ // - -500 * 64 for underflow, which when multiplied by scale = 2^-600 yields
+ // 2^-1100 and correctly underflows.
+ // However, in this correctly rounded version, such clamping can cause the
+ // Ziv accuracy test to fail on the clamped exponent, unnecessarily
+ // triggering the slow accurate paths for overflow or underflow cases.
+
+ constexpr double UPPER_EXP_BOUND = 511.0 * SCALE;
if (LIBC_UNLIKELY(FPBits(y6_log2_x.hi).abs().get_val() >= UPPER_EXP_BOUND)) {
if (FPBits(y6_log2_x.hi).sign() == Sign::POS) {
- scale = 0x1.0p512;
- y6_log2_x.hi -= 512.0 * 64.0;
- if (y6_log2_x.hi > 513.0 * 64.0)
- y6_log2_x.hi = 513.0 * 64.0;
+ if (y6_log2_x.hi > 1025.0 * SCALE)
+ return set_overflow(is_neg);
+ scale = 0x1.0p600;
+ y6_log2_x.hi -= 600.0 * SCALE;
} else {
- scale = 0x1.0p-512;
- y6_log2_x.hi += 512.0 * 64.0;
- if (y6_log2_x.hi < (-1076.0 + 512.0) * 64.0)
- y6_log2_x.hi = -564.0 * 64.0;
+ if (y6_log2_x.hi <= -1021.0 * SCALE) {
+ if (y6_log2_x.hi < -1076.0 * SCALE)
+ return set_underflow(is_neg);
+ is_denorm = true;
+ } else {
+ scale = 0x1.0p-600;
+ y6_log2_x.hi += 600.0 * SCALE;
+ }
}
}
double hm = fputil::nearest_integer(y6_log2_x.hi);
// lo6 = 2^6 * lo.
- double lo6_hi = y6_log2_x.hi - hm;
- double lo6 = lo6_hi + y6_log2_x.lo;
+ DoubleDouble lo6 = fputil::exact_add(y6_log2_x.hi - hm, y6_log2_x.lo);
int hm_i = static_cast<int>(hm);
unsigned idx_y = static_cast<unsigned>(hm_i) & 0x3f;
// 2^hi
- int64_t exp2_hi_i = static_cast<int64_t>(
- static_cast<uint64_t>(static_cast<int64_t>(hm_i >> 6))
- << FPBits::FRACTION_LEN);
+ int hi = hm_i >> 6;
+ double exp2_hi = 1.0;
+ if (LIBC_LIKELY(!is_denorm)) {
+ int64_t exp2_hi_i = static_cast<int64_t>(
+ static_cast<uint64_t>(hi + FPBits::EXP_BIAS) << FPBits::FRACTION_LEN);
+ exp2_hi = FPBits(static_cast<uint64_t>(exp2_hi_i)).get_val();
+ }
+
// 2^mid
- int64_t exp2_mid_hi_i =
- static_cast<int64_t>(FPBits(EXP2_MID1[idx_y].hi).uintval());
- int64_t exp2_mid_lo_i =
- static_cast<int64_t>(FPBits(EXP2_MID1[idx_y].mid).uintval());
- // (-1)^sign * 2^hi * 2^mid
- // Error <= 2^hi * 2^-53
- uint64_t exp2_hm_hi_i =
- static_cast<uint64_t>(exp2_hi_i + exp2_mid_hi_i) + sign;
- // The low part could be 0.
- uint64_t exp2_hm_lo_i =
- idx_y != 0 ? static_cast<uint64_t>(exp2_hi_i + exp2_mid_lo_i) + sign
- : sign;
- double exp2_hm_hi = FPBits(exp2_hm_hi_i).get_val();
- double exp2_hm_lo = FPBits(exp2_hm_lo_i).get_val();
-
- // Degree-5 polynomial approximation P(lo6) ~ 2^(lo6 / 2^6) = 2^(lo).
- // Generated by Sollya with:
- // > P = fpminimax(2^(x/64), 5, [|1, D...|], [-2^-1, 2^-1]);
- // > dirtyinfnorm(2^(x/64) - P, [-0.5, 0.5]);
- // 0x1.a2b77e618f5c4c176fd11b7659016cde5de83cb72p-60
- constexpr double EXP2_COEFFS[] = {0x1p0,
- 0x1.62e42fefa39efp-7,
- 0x1.ebfbdff82a23ap-15,
- 0x1.c6b08d7076268p-23,
- 0x1.3b2ad33f8b48bp-31,
- 0x1.5d870c4d84445p-40};
-
- double lo6_sqr = lo6 * lo6;
-
- double d0 = fputil::multiply_add(lo6, EXP2_COEFFS[2], EXP2_COEFFS[1]);
- double d1 = fputil::multiply_add(lo6, EXP2_COEFFS[4], EXP2_COEFFS[3]);
- double pp = fputil::polyeval(lo6_sqr, d0, d1, EXP2_COEFFS[5]);
-
- double r = fputil::multiply_add(exp2_hm_hi * lo6, pp, exp2_hm_lo);
- r += exp2_hm_hi;
-
- return r * scale;
+ DoubleDouble exp2_mid{EXP2_MID1[idx_y].mid * sign_d,
+ EXP2_MID1[idx_y].hi * sign_d};
+
+ // Polynomial expansion for 2^(lo6/64):
+ // 2^(lo6/64) ~ 1 + lo6 * (log(2)/64) + lo6^2 * P(lo6)
+ // The linear term is computed in DoubleDouble, and the degree-3 polynomial
+ // P(lo6) is evaluated with standard double precision.
+ //
+ // hi and lo parts of log(2)/64, generated by Sollya with:
+ // > a = D(log(2)/64);
+ // > b = D(log(2)/64 - a);
+ constexpr DoubleDouble LOG_2_OVER_64 = {0x1.abc9e3b39803fp-62,
+ 0x1.62e42fefa39efp-7};
+
+ // Degree-3 polynomial approximation for (2^(lo6/64) - 1 - lo6*log(2)/64) /
+ // lo6^2: Generated by Sollya with:
+ // > f = (2^(x/64) - 1 - x*log(2)/64) / x^2;
+ // > P = fpminimax(f, 5, [|D...|], [-0.5, 0.5]);
+ // > dirtyinfnorm((f - P) * x^2, [-0.5, 0.5]);
+ // 0x1.1778...p-73
+ constexpr double EXP2_COEFFS[] = {
+ 0x1.ebfbdff82c58fp-15, 0x1.c6b08d704a0cp-23, 0x1.3b2ab6fb4d08fp-31,
+ 0x1.5d87fe77a735dp-40, 0x1.430d835610044p-49, 0x1.ffe67c38112c3p-59};
+
+ DoubleDouble lo_log_2 = fputil::quick_mult(lo6, LOG_2_OVER_64);
+ DoubleDouble lo_log_2_p1 = fputil::exact_add(1.0, lo_log_2.hi);
+ lo_log_2_p1.lo += lo_log_2.lo;
+
+ double lo6_sq = lo6.hi * lo6.hi;
+ double e0 = fputil::multiply_add(lo6.hi, EXP2_COEFFS[1], EXP2_COEFFS[0]);
+ double e1 = fputil::multiply_add(lo6.hi, EXP2_COEFFS[3], EXP2_COEFFS[2]);
+ double e2 = fputil::multiply_add(lo6.hi, EXP2_COEFFS[5], EXP2_COEFFS[4]);
+
+ double lo6_4 = lo6_sq * lo6_sq;
+ double f0 = fputil::multiply_add(lo6_sq, e0, lo_log_2_p1.lo);
+ double f1 = fputil::multiply_add(lo6_sq, e2, e1);
+
+ lo_log_2_p1.lo = fputil::multiply_add(lo6_4, f1, f0);
+ DoubleDouble r = fputil::quick_mult(exp2_mid, lo_log_2_p1);
+
+ // Absolute error bound for r:
+ // - Base error from exp2 stage is bounded by 0x1.0p-64.
+ // - Propagated error from log2(x):
+ // |y| * err(log2(x)) * (2^mid * log(2)) <= |y| * 2^-73.5
+ // Dynamic error bound adapting to exponent scale |y|:
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ double err_r =
+ fputil::multiply_add(FPBits(y).abs().get_val(), 0x1.8p-73, 0x1.0p-64);
+#else
+ // Without FMA, intermediate roundings increase log2(x) and exp2 errors.
+ double err_r =
+ fputil::multiply_add(FPBits(y).abs().get_val(), 0x1.cp-72, 0x1.0p-63);
+#endif
----------------
lntue wrote:
Done.
https://github.com/llvm/llvm-project/pull/222827
More information about the libc-commits
mailing list