Merge bitcoin/bitcoin#30535: feefrac: add support for evaluating at given size

58914ab459 fuzz: assert min diff between FeeFrac and CFeeRate (Pieter Wuille)
0c6bcfd8f7 feefrac: support both rounding up and down for Evaluate (Pieter Wuille)
ecf956ec9d feefrac: add support for evaluating at given size (Pieter Wuille)
7963aecead feefrac: add helper functions for 96-bit division (Pieter Wuille)
800c0dea9a feefrac: rework comments around Mul/MulFallback (Pieter Wuille)
fcfe008db2 feefrac fuzz: use arith_uint256 instead of ad-hoc multiply (Pieter Wuille)
46ff4220bf arith_uint256: modernize comparison operators (Pieter Wuille)

Pull request description:

  The `FeeFrac` type represents a fraction, intended to be used for sats/vbyte or sats/WU. This PR adds functionality to evaluate that feerate for a given size, in order to obtain the fee it corresponds with (rounding down, or rounding up).

  The motivation here is being able to do accurate feerate evaluations in cluster mempool block building heuristics (where rounding down is needed), but in principle this makes it possible to use `FeeFrac` as a more accurate replacement for `CFeeRate` (where for feerate estimation rounding up is desirable). Because of this, both rounding modes are implemented.

  Unit tests are included for known-correct values, plus a fuzz test that verifies the result using `arith_uint256`.

ACKs for top commit:
  l0rinc:
    ACK 58914ab459
  ismaelsadeeq:
    reACK  58914ab459
  glozow:
    light code review ACK 58914ab459

Tree-SHA512: 362b88454bf355cae1f12d6430b1bb9ab66824140e12b27db7c48385f1e8db936da7d0694fb5aad2a00eb9e5fe3083a3a2c0cc40b2a68e2d37e07b3481d4eeae

src/arith_uint256.h

@@ -6,6 +6,7 @@
 #ifndef BITCOIN_ARITH_UINT256_H
 #define BITCOIN_ARITH_UINT256_H
 
+#include <compare>
 #include <cstdint>
 #include <cstring>
 #include <limits>
@@ -212,13 +213,8 @@ public:
     friend inline base_uint operator<<(const base_uint& a, int shift) { return base_uint(a) <<= shift; }
     friend inline base_uint operator*(const base_uint& a, uint32_t b) { return base_uint(a) *= b; }
     friend inline bool operator==(const base_uint& a, const base_uint& b) { return memcmp(a.pn, b.pn, sizeof(a.pn)) == 0; }
-    friend inline bool operator!=(const base_uint& a, const base_uint& b) { return memcmp(a.pn, b.pn, sizeof(a.pn)) != 0; }
-    friend inline bool operator>(const base_uint& a, const base_uint& b) { return a.CompareTo(b) > 0; }
-    friend inline bool operator<(const base_uint& a, const base_uint& b) { return a.CompareTo(b) < 0; }
-    friend inline bool operator>=(const base_uint& a, const base_uint& b) { return a.CompareTo(b) >= 0; }
-    friend inline bool operator<=(const base_uint& a, const base_uint& b) { return a.CompareTo(b) <= 0; }
+    friend inline std::strong_ordering operator<=>(const base_uint& a, const base_uint& b) { return a.CompareTo(b) <=> 0; }
     friend inline bool operator==(const base_uint& a, uint64_t b) { return a.EqualTo(b); }
-    friend inline bool operator!=(const base_uint& a, uint64_t b) { return !a.EqualTo(b); }
 
     /** Hex encoding of the number (with the most significant digits first). */
     std::string GetHex() const;

src/test/feefrac_tests.cpp

@@ -15,7 +15,45 @@ BOOST_AUTO_TEST_CASE(feefrac_operators)
     FeeFrac sum{1500, 400};
     FeeFrac diff{500, -200};
     FeeFrac empty{0, 0};
-    [[maybe_unused]] FeeFrac zero_fee{0, 1}; // zero-fee allowed
+    FeeFrac zero_fee{0, 1}; // zero-fee allowed
+
+    BOOST_CHECK_EQUAL(zero_fee.EvaluateFeeDown(0), 0);
+    BOOST_CHECK_EQUAL(zero_fee.EvaluateFeeDown(1), 0);
+    BOOST_CHECK_EQUAL(zero_fee.EvaluateFeeDown(1000000), 0);
+    BOOST_CHECK_EQUAL(zero_fee.EvaluateFeeDown(0x7fffffff), 0);
+    BOOST_CHECK_EQUAL(zero_fee.EvaluateFeeUp(0), 0);
+    BOOST_CHECK_EQUAL(zero_fee.EvaluateFeeUp(1), 0);
+    BOOST_CHECK_EQUAL(zero_fee.EvaluateFeeUp(1000000), 0);
+    BOOST_CHECK_EQUAL(zero_fee.EvaluateFeeUp(0x7fffffff), 0);
+
+    BOOST_CHECK_EQUAL(p1.EvaluateFeeDown(0), 0);
+    BOOST_CHECK_EQUAL(p1.EvaluateFeeDown(1), 10);
+    BOOST_CHECK_EQUAL(p1.EvaluateFeeDown(100000000), 1000000000);
+    BOOST_CHECK_EQUAL(p1.EvaluateFeeDown(0x7fffffff), int64_t(0x7fffffff) * 10);
+    BOOST_CHECK_EQUAL(p1.EvaluateFeeUp(0), 0);
+    BOOST_CHECK_EQUAL(p1.EvaluateFeeUp(1), 10);
+    BOOST_CHECK_EQUAL(p1.EvaluateFeeUp(100000000), 1000000000);
+    BOOST_CHECK_EQUAL(p1.EvaluateFeeUp(0x7fffffff), int64_t(0x7fffffff) * 10);
+
+    FeeFrac neg{-1001, 100};
+    BOOST_CHECK_EQUAL(neg.EvaluateFeeDown(0), 0);
+    BOOST_CHECK_EQUAL(neg.EvaluateFeeDown(1), -11);
+    BOOST_CHECK_EQUAL(neg.EvaluateFeeDown(2), -21);
+    BOOST_CHECK_EQUAL(neg.EvaluateFeeDown(3), -31);
+    BOOST_CHECK_EQUAL(neg.EvaluateFeeDown(100), -1001);
+    BOOST_CHECK_EQUAL(neg.EvaluateFeeDown(101), -1012);
+    BOOST_CHECK_EQUAL(neg.EvaluateFeeDown(100000000), -1001000000);
+    BOOST_CHECK_EQUAL(neg.EvaluateFeeDown(100000001), -1001000011);
+    BOOST_CHECK_EQUAL(neg.EvaluateFeeDown(0x7fffffff), -21496311307);
+    BOOST_CHECK_EQUAL(neg.EvaluateFeeUp(0), 0);
+    BOOST_CHECK_EQUAL(neg.EvaluateFeeUp(1), -10);
+    BOOST_CHECK_EQUAL(neg.EvaluateFeeUp(2), -20);
+    BOOST_CHECK_EQUAL(neg.EvaluateFeeUp(3), -30);
+    BOOST_CHECK_EQUAL(neg.EvaluateFeeUp(100), -1001);
+    BOOST_CHECK_EQUAL(neg.EvaluateFeeUp(101), -1011);
+    BOOST_CHECK_EQUAL(neg.EvaluateFeeUp(100000000), -1001000000);
+    BOOST_CHECK_EQUAL(neg.EvaluateFeeUp(100000001), -1001000010);
+    BOOST_CHECK_EQUAL(neg.EvaluateFeeUp(0x7fffffff), -21496311306);
 
     BOOST_CHECK(empty == FeeFrac{}); // same as no-args
 
@@ -67,6 +105,23 @@ BOOST_AUTO_TEST_CASE(feefrac_operators)
     BOOST_CHECK(oversized_1 << oversized_2);
     BOOST_CHECK(oversized_1 != oversized_2);
 
+    BOOST_CHECK_EQUAL(oversized_1.EvaluateFeeDown(0), 0);
+    BOOST_CHECK_EQUAL(oversized_1.EvaluateFeeDown(1), 1152921);
+    BOOST_CHECK_EQUAL(oversized_1.EvaluateFeeDown(2), 2305843);
+    BOOST_CHECK_EQUAL(oversized_1.EvaluateFeeDown(1548031267), 1784758530396540);
+    BOOST_CHECK_EQUAL(oversized_1.EvaluateFeeUp(0), 0);
+    BOOST_CHECK_EQUAL(oversized_1.EvaluateFeeUp(1), 1152922);
+    BOOST_CHECK_EQUAL(oversized_1.EvaluateFeeUp(2), 2305843);
+    BOOST_CHECK_EQUAL(oversized_1.EvaluateFeeUp(1548031267), 1784758530396541);
+
+    // Test cases on the threshold where FeeFrac::Evaluate start using Mul/Div.
+    BOOST_CHECK_EQUAL(FeeFrac(0x1ffffffff, 123456789).EvaluateFeeDown(98765432), 6871947728);
+    BOOST_CHECK_EQUAL(FeeFrac(0x200000000, 123456789).EvaluateFeeDown(98765432), 6871947729);
+    BOOST_CHECK_EQUAL(FeeFrac(0x200000001, 123456789).EvaluateFeeDown(98765432), 6871947730);
+    BOOST_CHECK_EQUAL(FeeFrac(0x1ffffffff, 123456789).EvaluateFeeUp(98765432), 6871947729);
+    BOOST_CHECK_EQUAL(FeeFrac(0x200000000, 123456789).EvaluateFeeUp(98765432), 6871947730);
+    BOOST_CHECK_EQUAL(FeeFrac(0x200000001, 123456789).EvaluateFeeUp(98765432), 6871947731);
+
     // Tests paths that use double arithmetic
     FeeFrac busted{(static_cast<int64_t>(INT32_MAX)) + 1, INT32_MAX};
     BOOST_CHECK(!(busted < busted));
@@ -77,6 +132,19 @@ BOOST_AUTO_TEST_CASE(feefrac_operators)
     BOOST_CHECK(max_fee <= max_fee);
     BOOST_CHECK(max_fee >= max_fee);
 
+    BOOST_CHECK_EQUAL(max_fee.EvaluateFeeDown(0), 0);
+    BOOST_CHECK_EQUAL(max_fee.EvaluateFeeDown(1), 977888);
+    BOOST_CHECK_EQUAL(max_fee.EvaluateFeeDown(2), 1955777);
+    BOOST_CHECK_EQUAL(max_fee.EvaluateFeeDown(3), 2933666);
+    BOOST_CHECK_EQUAL(max_fee.EvaluateFeeDown(1256796054), 1229006664189047);
+    BOOST_CHECK_EQUAL(max_fee.EvaluateFeeDown(INT32_MAX), 2100000000000000);
+    BOOST_CHECK_EQUAL(max_fee.EvaluateFeeUp(0), 0);
+    BOOST_CHECK_EQUAL(max_fee.EvaluateFeeUp(1), 977889);
+    BOOST_CHECK_EQUAL(max_fee.EvaluateFeeUp(2), 1955778);
+    BOOST_CHECK_EQUAL(max_fee.EvaluateFeeUp(3), 2933667);
+    BOOST_CHECK_EQUAL(max_fee.EvaluateFeeUp(1256796054), 1229006664189048);
+    BOOST_CHECK_EQUAL(max_fee.EvaluateFeeUp(INT32_MAX), 2100000000000000);
+
     FeeFrac max_fee2{1, 1};
     BOOST_CHECK(max_fee >= max_fee2);
 

src/test/fuzz/feefrac.cpp

@@ -2,54 +2,48 @@
 // Distributed under the MIT software license, see the accompanying
 // file COPYING or http://www.opensource.org/licenses/mit-license.php.
 
+#include <arith_uint256.h>
+#include <policy/feerate.h>
 #include <util/feefrac.h>
 #include <test/fuzz/FuzzedDataProvider.h>
 #include <test/fuzz/fuzz.h>
 #include <test/fuzz/util.h>
 
 #include <compare>
+#include <cmath>
 #include <cstdint>
 #include <iostream>
 
 namespace {
 
-/** Compute a * b, represented in 4x32 bits, highest limb first. */
-std::array<uint32_t, 4> Mul128(uint64_t a, uint64_t b)
+/** The maximum absolute value of an int64_t, as an arith_uint256 (2^63). */
+const auto MAX_ABS_INT64 = arith_uint256{1} << 63;
+
+/** Construct an arith_uint256 whose value equals abs(x). */
+arith_uint256 Abs256(int64_t x)
 {
-    std::array<uint32_t, 4> ret{0, 0, 0, 0};
-
-    /** Perform ret += v << (32 * pos), at 128-bit precision. */
-    auto add_fn = [&](uint64_t v, int pos) {
-        uint64_t accum{0};
-        for (int i = 0; i + pos < 4; ++i) {
-            // Add current value at limb pos in ret.
-            accum += ret[3 - pos - i];
-            // Add low or high half of v.
-            if (i == 0) accum += v & 0xffffffff;
-            if (i == 1) accum += v >> 32;
-            // Store lower half of result in limb pos in ret.
-            ret[3 - pos - i] = accum & 0xffffffff;
-            // Leave carry in accum.
-            accum >>= 32;
-        }
-        // Make sure no overflow.
-        assert(accum == 0);
-    };
-
-    // Multiply the 4 individual limbs (schoolbook multiply, with base 2^32).
-    add_fn((a & 0xffffffff) * (b & 0xffffffff), 0);
-    add_fn((a >> 32) * (b & 0xffffffff), 1);
-    add_fn((a & 0xffffffff) * (b >> 32), 1);
-    add_fn((a >> 32) * (b >> 32), 2);
-    return ret;
+    if (x >= 0) {
+        // For positive numbers, pass through the value.
+        return arith_uint256{static_cast<uint64_t>(x)};
+    } else if (x > std::numeric_limits<int64_t>::min()) {
+        // For negative numbers, negate first.
+        return arith_uint256{static_cast<uint64_t>(-x)};
+    } else {
+        // Special case for x == -2^63 (for which -x results in integer overflow).
+        return MAX_ABS_INT64;
+    }
 }
 
-/* comparison helper for std::array */
-std::strong_ordering compare_arrays(const std::array<uint32_t, 4>& a, const std::array<uint32_t, 4>& b) {
-    for (size_t i = 0; i < a.size(); ++i) {
-        if (a[i] != b[i]) return a[i] <=> b[i];
+/** Construct an arith_uint256 whose value equals abs(x), for 96-bit x. */
+arith_uint256 Abs256(std::pair<int64_t, uint32_t> 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);
     }
-    return std::strong_ordering::equal;
 }
 
 std::strong_ordering MulCompare(int64_t a1, int64_t a2, int64_t b1, int64_t b2)
@@ -59,23 +53,14 @@ std::strong_ordering MulCompare(int64_t a1, int64_t a2, int64_t b1, int64_t b2)
     int sign_b = (b1 == 0 ? 0 : b1 < 0 ? -1 : 1) * (b2 == 0 ? 0 : b2 < 0 ? -1 : 1);
     if (sign_a != sign_b) return sign_a <=> sign_b;
 
-    // Compute absolute values.
-    uint64_t abs_a1 = static_cast<uint64_t>(a1), abs_a2 = static_cast<uint64_t>(a2);
-    uint64_t abs_b1 = static_cast<uint64_t>(b1), abs_b2 = static_cast<uint64_t>(b2);
-    // Use (~x + 1) instead of the equivalent (-x) to silence the linter; mod 2^64 behavior is
-    // intentional here.
-    if (a1 < 0) abs_a1 = ~abs_a1 + 1;
-    if (a2 < 0) abs_a2 = ~abs_a2 + 1;
-    if (b1 < 0) abs_b1 = ~abs_b1 + 1;
-    if (b2 < 0) abs_b2 = ~abs_b2 + 1;
+    // Compute absolute values of products.
+    auto mul_abs_a = Abs256(a1) * Abs256(a2), mul_abs_b = Abs256(b1) * Abs256(b2);
 
     // Compute products of absolute values.
-    auto mul_abs_a = Mul128(abs_a1, abs_a2);
-    auto mul_abs_b = Mul128(abs_b1, abs_b2);
     if (sign_a < 0) {
-        return compare_arrays(mul_abs_b, mul_abs_a);
+        return mul_abs_b <=> mul_abs_a;
     } else {
-        return compare_arrays(mul_abs_a, mul_abs_b);
+        return mul_abs_a <=> mul_abs_b;
     }
 }
 
@@ -121,3 +106,128 @@ 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, a positive 31-bit value den, and rounding mode.
+    FuzzedDataProvider provider(buffer.data(), buffer.size());
+    auto num_high = provider.ConsumeIntegral<int64_t>();
+    auto num_low = provider.ConsumeIntegral<uint32_t>();
+    std::pair<int64_t, uint32_t> num{num_high, num_low};
+    auto den = provider.ConsumeIntegralInRange<int32_t>(1, std::numeric_limits<int32_t>::max());
+    auto round_down = provider.ConsumeBool();
+
+    // 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 and we are
+    // rounding down, or positive and we are rounding up, 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 == round_down) ?
+        (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, round_down);
+    assert(res == 0 || (res < 0) == is_negative);
+    assert(Abs256(res) == quot_abs);
+
+    // Compare approximately with floating-point.
+    long double expect = round_down ? std::floor(num_high * 4294967296.0L + num_low) / den
+                                    : std::ceil(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, a positive 31-bit
+    // divisor, and a rounding mode.
+    FuzzedDataProvider provider(buffer.data(), buffer.size());
+    auto mul32 = provider.ConsumeIntegral<int32_t>();
+    auto mul64 = provider.ConsumeIntegral<int64_t>();
+    auto div = provider.ConsumeIntegralInRange<int32_t>(1, std::numeric_limits<int32_t>::max());
+    auto round_down = provider.ConsumeBool();
+
+    // 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 and we are
+    // rounding down or positive and we rounding up, 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 == round_down) ?
+        (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. In the context
+        // of the Evaluate{Down,Up} calls below, this corresponds to 0 <= at_size <= feefrac.size.
+        assert(mul32 < 0 || mul32 > div);
+        return;
+    }
+
+    // Verify the behavior of FeeFrac::Mul + FeeFrac::Div.
+    auto res = FeeFrac::Div(FeeFrac::Mul(mul64, mul32), div, round_down);
+    assert(res == 0 || (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, round_down);
+    assert(res == res_fallback);
+
+    // Compare approximately with floating-point.
+    long double expect = round_down ? std::floor(static_cast<long double>(mul32) * mul64 / div)
+                                    : std::ceil(static_cast<long double>(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);
+    }
+
+    // Verify the behavior of FeeFrac::Evaluate{Down,Up}.
+    if (mul32 >= 0) {
+        auto res_fee = round_down ?
+            FeeFrac{mul64, div}.EvaluateFeeDown(mul32) :
+            FeeFrac{mul64, div}.EvaluateFeeUp(mul32);
+        assert(res == res_fee);
+
+        // Compare approximately with CFeeRate.
+        if (mul64 <= std::numeric_limits<int64_t>::max() / 1000 &&
+            mul64 >= std::numeric_limits<int64_t>::min() / 1000 &&
+            quot_abs <= arith_uint256{std::numeric_limits<int64_t>::max() / 1000}) {
+            CFeeRate feerate(mul64, (uint32_t)div);
+            CAmount feerate_fee{feerate.GetFee(mul32)};
+            auto allowed_gap = static_cast<int64_t>(mul32 / 1000 + 3 + round_down);
+            assert(feerate_fee - res_fee >= -allowed_gap);
+            assert(feerate_fee - res_fee <= allowed_gap);
+        }
+    }
+}

src/util/feefrac.h

@@ -37,27 +37,70 @@
  */
 struct FeeFrac
 {
-    /** Fallback version for Mul (see below).
-     *
-     * Separate to permit testing on platforms where it isn't actually needed.
-     */
+    /** Helper function for 32*64 signed multiplication, returning an unspecified but totally
+     *  ordered type. This is a fallback version, separate so it can be tested on platforms where
+     *  it isn't actually needed. */
     static inline std::pair<int64_t, uint32_t> MulFallback(int64_t a, int32_t b) noexcept
     {
-        // Otherwise, emulate 96-bit multiplication using two 64-bit multiplies.
         int64_t low = int64_t{static_cast<uint32_t>(a)} * b;
         int64_t high = (a >> 32) * b;
         return {high + (low >> 32), static_cast<uint32_t>(low)};
     }
 
-    // Compute a * b, returning an unspecified but totally ordered type.
+    /** Helper function for 96/32 signed division, rounding towards negative infinity (if
+     *  round_down) or positive infinity (if !round_down). This is a fallback version, separate so
+     *  that 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 be consistent for testability reasons.
+     *
+     * The result must fit in an int64_t, and d must be strictly positive. */
+    static inline int64_t DivFallback(std::pair<int64_t, uint32_t> n, int32_t d, bool round_down) 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. It is possible
+        // that the / operator here rounds in the wrong direction (if n_low is not a multiple of
+        // size, and is (if round_down) negative, or (if !round_down) positive). If so, make a
+        // correction.
+        int64_t quot_low = n_low / d;
+        int32_t mod_low = n_low % d;
+        quot_low += (mod_low > 0) - (mod_low && round_down);
+        // 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. */
     static inline __int128 Mul(int64_t a, int32_t b) noexcept
     {
-        // If __int128 is available, use 128-bit wide multiply.
         return __int128{a} * b;
     }
+
+    /** Helper function for 96/32 signed division, rounding towards negative infinity (if
+     *  round_down), or towards positive infinity (if !round_down). 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, bool round_down) noexcept
+    {
+        Assume(d > 0);
+        // Compute the division.
+        int64_t quot = n / d;
+        int32_t mod = n % d;
+        // Correct result if the / operator above rounded in the wrong direction.
+        return quot + (mod > 0) - (mod && round_down);
+    }
 #else
     static constexpr auto Mul = MulFallback;
+    static constexpr auto Div = DivFallback;
 #endif
 
     int64_t fee;
@@ -144,6 +187,39 @@ struct FeeFrac
         std::swap(a.fee, b.fee);
         std::swap(a.size, b.size);
     }
+
+    /** Compute the fee for a given size `at_size` using this object's feerate.
+     *
+     * This effectively corresponds to evaluating (this->fee * at_size) / this->size, with the
+     * result rounded towards negative infinity (if RoundDown) or towards positive infinity
+     * (if !RoundDown).
+     *
+     * Requires this->size > 0, at_size >= 0, and that the correct result fits in a int64_t. This
+     * is guaranteed to be the case when 0 <= at_size <= this->size.
+     */
+    template<bool RoundDown>
+    int64_t EvaluateFee(int32_t at_size) const noexcept
+    {
+        Assume(size > 0);
+        Assume(at_size >= 0);
+        if (fee >= 0 && fee < 0x200000000) [[likely]] {
+            // Common case where (this->fee * at_size) is guaranteed to fit in a uint64_t.
+            if constexpr (RoundDown) {
+                return (uint64_t(fee) * at_size) / uint32_t(size);
+            } else {
+                return (uint64_t(fee) * at_size + size - 1U) / uint32_t(size);
+            }
+        } else {
+            // Otherwise, use Mul and Div.
+            return Div(Mul(fee, at_size), size, RoundDown);
+        }
+    }
+
+public:
+    /** Compute the fee for a given size `at_size` using this object's feerate, rounding down. */
+    int64_t EvaluateFeeDown(int32_t at_size) const noexcept { return EvaluateFee<true>(at_size); }
+    /** Compute the fee for a given size `at_size` using this object's feerate, rounding up. */
+    int64_t EvaluateFeeUp(int32_t at_size) const noexcept { return EvaluateFee<false>(at_size); }
 };
 
 /** Compare the feerate diagrams implied by the provided sorted chunks data.