From 7963aecead968d126545d3730da5d9942c3f9518 Mon Sep 17 00:00:00 2001 From: Pieter Wuille Date: Fri, 26 Jul 2024 15:04:05 -0400 Subject: [PATCH] feefrac: add helper functions for 96-bit division These functions are needed to implement FeeFrac evaluation later: given a FeeFrac{fee, size}, its fee at at_size is (fee * at_size / size). --- src/test/fuzz/feefrac.cpp | 108 ++++++++++++++++++++++++++++++++++++++ src/util/feefrac.h | 40 ++++++++++++++ 2 files changed, 148 insertions(+) diff --git a/src/test/fuzz/feefrac.cpp b/src/test/fuzz/feefrac.cpp index 6363eb530eb..28ecfb14a55 100644 --- a/src/test/fuzz/feefrac.cpp +++ b/src/test/fuzz/feefrac.cpp @@ -9,6 +9,7 @@ #include #include +#include #include #include @@ -32,6 +33,18 @@ arith_uint256 Abs256(int64_t x) } } +/** Construct an arith_uint256 whose value equals abs(x), for 96-bit x. */ +arith_uint256 Abs256(std::pair x) +{ + if (x.first >= 0) { + // x.first and x.second are both non-negative; sum their absolute values. + return (Abs256(x.first) << 32) + Abs256(x.second); + } else { + // x.first is negative and x.second is non-negative; subtract the absolute values. + return (Abs256(x.first) << 32) - Abs256(x.second); + } +} + std::strong_ordering MulCompare(int64_t a1, int64_t a2, int64_t b1, int64_t b2) { // Compute and compare signs. @@ -92,3 +105,98 @@ FUZZ_TARGET(feefrac) assert((fr1 == fr2) == std::is_eq(cmp_total)); assert((fr1 != fr2) == std::is_neq(cmp_total)); } + +FUZZ_TARGET(feefrac_div_fallback) +{ + // Verify the behavior of FeeFrac::DivFallback over all possible inputs. + + // Construct a 96-bit signed value num, and positive 31-bit value den. + FuzzedDataProvider provider(buffer.data(), buffer.size()); + auto num_high = provider.ConsumeIntegral(); + auto num_low = provider.ConsumeIntegral(); + std::pair num{num_high, num_low}; + auto den = provider.ConsumeIntegralInRange(1, std::numeric_limits::max()); + + // Predict the sign of the actual result. + bool is_negative = num_high < 0; + // Evaluate absolute value using arith_uint256. If the actual result is negative, the absolute + // value of the quotient is the rounded-up quotient of the absolute values. + auto num_abs = Abs256(num); + auto den_abs = Abs256(den); + auto quot_abs = is_negative ? (num_abs + den_abs - 1) / den_abs : num_abs / den_abs; + + // If the result is not representable by an int64_t, bail out. + if ((is_negative && quot_abs > MAX_ABS_INT64) || (!is_negative && quot_abs >= MAX_ABS_INT64)) { + return; + } + + // Verify the behavior of FeeFrac::DivFallback. + auto res = FeeFrac::DivFallback(num, den); + assert((res < 0) == is_negative); + assert(Abs256(res) == quot_abs); + + // Compare approximately with floating-point. + long double expect = std::floorl(num_high * 4294967296.0L + num_low) / den; + // Expect to be accurate within 50 bits of precision, +- 1 sat. + if (expect == 0.0L) { + assert(res >= -1 && res <= 1); + } else if (expect > 0.0L) { + assert(res >= expect * 0.999999999999999L - 1.0L); + assert(res <= expect * 1.000000000000001L + 1.0L); + } else { + assert(res >= expect * 1.000000000000001L - 1.0L); + assert(res <= expect * 0.999999999999999L + 1.0L); + } +} + +FUZZ_TARGET(feefrac_mul_div) +{ + // Verify the behavior of: + // - The combination of FeeFrac::Mul + FeeFrac::Div. + // - The combination of FeeFrac::MulFallback + FeeFrac::DivFallback. + // - FeeFrac::Evaluate. + + // Construct a 32-bit signed multiplicand, a 64-bit signed multiplicand, and a positive 31-bit + // divisor. + FuzzedDataProvider provider(buffer.data(), buffer.size()); + auto mul32 = provider.ConsumeIntegral(); + auto mul64 = provider.ConsumeIntegral(); + auto div = provider.ConsumeIntegralInRange(1, std::numeric_limits::max()); + + // Predict the sign of the overall result. + bool is_negative = ((mul32 < 0) && (mul64 > 0)) || ((mul32 > 0) && (mul64 < 0)); + // Evaluate absolute value using arith_uint256. If the actual result is negative, the absolute + // value of the quotient is the rounded-up quotient of the absolute values. + auto prod_abs = Abs256(mul32) * Abs256(mul64); + auto div_abs = Abs256(div); + auto quot_abs = is_negative ? (prod_abs + div_abs - 1) / div_abs : prod_abs / div_abs; + + // If the result is not representable by an int64_t, bail out. + if ((is_negative && quot_abs > MAX_ABS_INT64) || (!is_negative && quot_abs >= MAX_ABS_INT64)) { + // If 0 <= mul32 <= div, then the result is guaranteed to be representable. + assert(mul32 < 0 || mul32 > div); + return; + } + + // Verify the behavior of FeeFrac::Mul + FeeFrac::Div. + auto res = FeeFrac::Div(FeeFrac::Mul(mul64, mul32), div); + assert((res < 0) == is_negative); + assert(Abs256(res) == quot_abs); + + // Verify the behavior of FeeFrac::MulFallback + FeeFrac::DivFallback. + auto res_fallback = FeeFrac::DivFallback(FeeFrac::MulFallback(mul64, mul32), div); + assert(res == res_fallback); + + // Compare approximately with floating-point. + long double expect = std::floorl(static_cast(mul32) * mul64 / div); + // Expect to be accurate within 50 bits of precision, +- 1 sat. + if (expect == 0.0L) { + assert(res >= -1 && res <= 1); + } else if (expect > 0.0L) { + assert(res >= expect * 0.999999999999999L - 1.0L); + assert(res <= expect * 1.000000000000001L + 1.0L); + } else { + assert(res >= expect * 1.000000000000001L - 1.0L); + assert(res <= expect * 0.999999999999999L + 1.0L); + } +} diff --git a/src/util/feefrac.h b/src/util/feefrac.h index 00a6a427e68..0f32cd50ed2 100644 --- a/src/util/feefrac.h +++ b/src/util/feefrac.h @@ -47,6 +47,32 @@ struct FeeFrac return {high + (low >> 32), static_cast(low)}; } + /** Helper function for 96/32 signed division, rounding towards negative infinity. This is a + * fallback version, separate so it can be tested on platforms where it isn't actually needed. + * + * The exact behavior with negative n does not really matter, but this implementation chooses + * to always round down, for consistency and testability. + * + * The result must fit in an int64_t, and d must be strictly positive. */ + static inline int64_t DivFallback(std::pair n, int32_t d) noexcept + { + Assume(d > 0); + // Compute quot_high = n.first / d, so the result becomes + // (n.second + (n.first - quot_high * d) * 2**32) / d + (quot_high * 2**32), or + // (n.second + (n.first % d) * 2**32) / d + (quot_high * 2**32). + int64_t quot_high = n.first / d; + // Evaluate the parenthesized expression above, so the result becomes + // n_low / d + (quot_high * 2**32) + int64_t n_low = ((n.first % d) << 32) + n.second; + // Evaluate the division so the result becomes quot_low + quot_high * 2**32. We need this + // division to round down however, while the / operator rounds towards zero. In case n_low + // is negative and not a multiple of size, we thus need a correction. + int64_t quot_low = n_low / d; + quot_low -= (n_low % d) < 0; + // Combine and return the result + return (quot_high << 32) + quot_low; + } + #ifdef __SIZEOF_INT128__ /** Helper function for 32*64 signed multiplication, returning an unspecified but totally * ordered type. This is a version relying on __int128. */ @@ -54,8 +80,22 @@ struct FeeFrac { return __int128{a} * b; } + + /** Helper function for 96/32 signed division, rounding towards negative infinity. This is a + * version relying on __int128. + * + * The result must fit in an int64_t, and d must be strictly positive. */ + static inline int64_t Div(__int128 n, int32_t d) noexcept + { + Assume(d > 0); + // Compute the division. + int64_t quot = n / d; + // Make it round down. + return quot - ((n % d) < 0); + } #else static constexpr auto Mul = MulFallback; + static constexpr auto Div = DivFallback; #endif int64_t fee;