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).

Author: sipa

src/test/fuzz/feefrac.cpp

@@ -9,6 +9,7 @@
 #include <test/fuzz/util.h>
 
 #include <compare>
+#include <cmath>
 #include <cstdint>
 #include <iostream>
 
@@ -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<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);
+    }
+}
+
 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<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());
+
+    // 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<int32_t>();
+    auto mul64 = provider.ConsumeIntegral<int64_t>();
+    auto div = provider.ConsumeIntegralInRange<int32_t>(1, std::numeric_limits<int32_t>::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<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);
+    }
+}

src/util/feefrac.h

@@ -47,15 +47,55 @@ struct FeeFrac
         return {high + (low >> 32), static_cast<uint32_t>(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<int64_t, uint32_t> 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. */
     static inline __int128 Mul(int64_t a, int32_t b) noexcept
     {
         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;