From d6a544877fb84c93050fb00eeb9071e7c3554f4d Mon Sep 17 00:00:00 2001 From: York Cao <978176728@qq.com> Date: Mon, 3 Aug 2026 20:13:53 +0800 Subject: [PATCH] [opt](decimal) fast paths for wide-integer (Decimal256) division ### What problem does this PR solve? Issue Number: no issue Problem Summary: `wide::integer` division (backing Decimal256 and other >128-bit integer types) always fell back to a generic bit-by-bit binary long-division loop that iterates ~Bits times (256 shift/compare/subtract rounds for a 256-bit value), even when the operands are small. Decimal256 arithmetic and the CEIL/ROUND `x / 10^k` rounding paths hit this hot loop constantly with operands that are far narrower than 256 bits, so the general algorithm dominates the cost. This adds three stacked fast paths in front of the generic loop, each bit-exact with it (verified against native `__int128` oracles and via q*d+r==n identities), and each writing the remainder back into `numerator` so `operator%` stays correct: 1. Both operands fit in 128 bits (the common money/count magnitude): perform a single native `unsigned __int128` divide. Placed first because it is the cheapest and most frequently hit. 2. Divisor fits in a single 64-bit limb (e.g. `x / 10^k`): schoolbook word-by-word long division, one hardware 128/64 divide per limb -- O(item_count) divides instead of ~Bits iterations. 3. Divisor fits in two 64-bit limbs (65..128-bit divisor): route through a new `divide_knuth()` helper implementing Knuth's Algorithm D (Hacker's Delight `divmnu`) in base 2^32, which keeps every intermediate product within a uint64_t and stays overflow-safe. Divisors wider than 128 bits and the zero-divisor throw are unchanged and fall through to the existing generic path. On a fast-path miss the only added cost is a short limb scan. ### Release note None ### Check List (For Author) - Test: Unit Test - Added/extended `be/test/core/wide_integer_test.cpp` (18 new cases: single-limb, two-limb Knuth, both-fit-128, signed, boundary, divide-by-zero, and randomized differential/ground-truth fuzz against native __int128). All 23 WideInteger tests pass locally (ASAN build). - Behavior changed: No (pure performance optimization; results are bit-exact with the previous slow path) - Does this need documentation: No --- be/src/core/wide_integer_impl.h | 234 ++++++++++++++ be/test/core/wide_integer_test.cpp | 469 +++++++++++++++++++++++++++++ 2 files changed, 703 insertions(+) diff --git a/be/src/core/wide_integer_impl.h b/be/src/core/wide_integer_impl.h index a5ad6ae041b0b7..e493ec19939cf9 100644 --- a/be/src/core/wide_integer_impl.h +++ b/be/src/core/wide_integer_impl.h @@ -831,6 +831,136 @@ struct integer::_impl { return is_zero; } + /// Multi-word unsigned division via Knuth's Algorithm D, a faithful port of + /// Hacker's Delight `divmnu` (Fig. 9-2, Henry S. Warren) run in base 2^32 over + /// 32-bit half-limbs. Working in base 2^32 keeps every intermediate product + /// (qhat * vn[i], with both factors < 2^32) within a uint64_t, so no 128-bit + /// arithmetic is needed and the estimate-correction stays overflow-safe. + /// Preconditions: denominator is non-zero and occupies at least two 32-bit + /// half-limbs (i.e. >= 2^32). Returns quotient; leaves remainder in numerator, + /// matching the divide() contract. + template + constexpr static integer divide_knuth( + integer& numerator, const integer& denominator) { + constexpr unsigned half_count = Bits2 / 32; // 32-bit half-limbs (8 for 256-bit) + constexpr uint64_t base = uint64_t(1) << 32; + + // Unpack the 64-bit limbs into little-endian 32-bit half-limbs. + uint32_t u[half_count]; + uint32_t v[half_count]; + for (unsigned i = 0; i < item_count; ++i) { + const uint64_t nu = numerator.items[little(i)]; + const uint64_t de = denominator.items[little(i)]; + u[2 * i] = static_cast(nu); + u[2 * i + 1] = static_cast(nu >> 32); + v[2 * i] = static_cast(de); + v[2 * i + 1] = static_cast(de >> 32); + } + + int n = half_count; + while (n > 0 && v[n - 1] == 0) { + --n; + } + int m = half_count; + while (m > 0 && u[m - 1] == 0) { + --m; + } + + integer quotient = 0; + // Dividend shorter than divisor: quotient is 0 and remainder is the dividend + // (already sitting in numerator), so nothing else to do. + if (m < n) { + return quotient; + } + + // Normalize: left-shift so the divisor's top half-limb has its high bit set. + int s = 0; + for (uint32_t hi = v[n - 1]; (hi & 0x80000000u) == 0u; hi <<= 1) { + ++s; + } + + uint32_t vn[half_count]; + for (int i = n - 1; i > 0; --i) { + vn[i] = (v[i] << s) | (s == 0 ? 0u : (v[i - 1] >> (32 - s))); + } + vn[0] = v[0] << s; + + uint32_t un[half_count + 1]; + un[m] = (s == 0) ? 0u : (u[m - 1] >> (32 - s)); + for (int i = m - 1; i > 0; --i) { + un[i] = (u[i] << s) | (s == 0 ? 0u : (u[i - 1] >> (32 - s))); + } + un[0] = u[0] << s; + + uint32_t q[half_count]; + for (unsigned i = 0; i < half_count; ++i) { + q[i] = 0; + } + + for (int j = m - n; j >= 0; --j) { + // Estimate quotient digit qhat and its remainder rhat, then correct a + // possible overestimate. The short-circuit keeps qhat * vn[n-2] from + // being evaluated (and overflowing) while qhat >= base. + const uint64_t num_top = static_cast(un[j + n]) * base + un[j + n - 1]; + uint64_t qhat = num_top / vn[n - 1]; + uint64_t rhat = num_top - qhat * vn[n - 1]; + while (qhat >= base || qhat * vn[n - 2] > base * rhat + un[j + n - 2]) { + --qhat; + rhat += vn[n - 1]; + if (rhat >= base) { + break; + } + } + + // Multiply the divisor by qhat and subtract from the running dividend, + // propagating the borrow through k (signed). + int64_t k = 0; + int64_t t = 0; + for (int i = 0; i < n; ++i) { + const uint64_t p = qhat * vn[i]; + t = static_cast(un[i + j]) - k - static_cast(p & 0xFFFFFFFFu); + un[i + j] = static_cast(t); + k = static_cast(p >> 32) - (t >> 32); + } + t = static_cast(un[j + n]) - k; + un[j + n] = static_cast(t); + + q[j] = static_cast(qhat); + if (t < 0) { + // qhat was one too large: decrement and add the divisor back. + --q[j]; + uint64_t carry = 0; + for (int i = 0; i < n; ++i) { + const uint64_t sum = static_cast(un[i + j]) + vn[i] + carry; + un[i + j] = static_cast(sum); + carry = sum >> 32; + } + un[j + n] = static_cast(static_cast(un[j + n]) + carry); + } + } + + // Repack quotient half-limbs into 64-bit limbs. + for (unsigned i = 0; i < item_count; ++i) { + quotient.items[little(i)] = + static_cast(q[2 * i]) | (static_cast(q[2 * i + 1]) << 32); + } + + // De-normalize the remainder (right-shift by s) back into numerator. + uint32_t rem[half_count]; + for (unsigned i = 0; i < half_count; ++i) { + rem[i] = 0; + } + for (int i = 0; i < n; ++i) { + rem[i] = (un[i] >> s) | (s == 0 ? 0u : (un[i + 1] << (32 - s))); + } + for (unsigned i = 0; i < item_count; ++i) { + numerator.items[little(i)] = static_cast(rem[2 * i]) | + (static_cast(rem[2 * i + 1]) << 32); + } + + return quotient; + } + /// returns quotient as result and remainder in numerator. template constexpr static integer divide(integer& numerator, @@ -861,10 +991,114 @@ struct integer::_impl { return res; } + // Fast path for wider integers (Decimal256) whose operands BOTH fit in 128 + // bits: money/count magnitudes almost always do. A single hardware __int128 + // divide suffices. This layer sits ahead of the divisor-width fast paths + // below on purpose: when the divisor is 65..128 bits those paths route to + // Knuth Algorithm D (a multi-word estimate/correct loop), whereas here -- if + // the dividend is also <= 128 bits -- one native divide wins outright. + // Bit-exact with the slow path, and the remainder written back to numerator + // keeps operator_percent correct. On a miss the extra cost is only the limb + // scan, then control falls through to the divisor-width fast paths and the + // generic loop below (no regression). + if constexpr (Bits > 128 && sizeof(base_type) == 8) { + bool operands_fit_128 = true; + for (unsigned i = 2; i < item_count; ++i) { + if (numerator.items[little(i)] != 0 || denominator.items[little(i)] != 0) { + operands_fit_128 = false; + break; + } + } + if (operands_fit_128) { + using CompilerUInt128 = unsigned __int128; + CompilerUInt128 b = (CompilerUInt128(denominator.items[little(1)]) << 64) + + denominator.items[little(0)]; + // A zero denominator falls through to the throw below; never divide by it. + if (b != 0) { + CompilerUInt128 a = (CompilerUInt128(numerator.items[little(1)]) << 64) + + numerator.items[little(0)]; + CompilerUInt128 c = a / b; + + integer res; + res.items[little(0)] = static_cast(c); + res.items[little(1)] = static_cast(c >> 64); + for (unsigned i = 2; i < item_count; ++i) { + res.items[little(i)] = 0; + } + + CompilerUInt128 remainder = a - b * c; + numerator.items[little(0)] = static_cast(remainder); + numerator.items[little(1)] = static_cast(remainder >> 64); + for (unsigned i = 2; i < item_count; ++i) { + numerator.items[little(i)] = 0; + } + + return res; + } + } + } + if (is_zero(denominator)) { throw doris::Exception(doris::ErrorCode::INVALID_ARGUMENT, "Division by zero"); } + /// Fast path for a divisor that fits in a single 64-bit limb: schoolbook long + /// division word-by-word from the most- to least-significant limb. Each step + /// divides {remainder : current_word} (128 bits) by the 64-bit divisor, which + /// lowers to a single hardware divide (e.g. `divq` on x86-64). O(item_count) + /// divides vs ~Bits iterations of the generic binary long division below. + { + bool divisor_is_single_limb = true; + for (unsigned i = 1; i < item_count; ++i) { + if (denominator.items[little(i)] != 0) { + divisor_is_single_limb = false; + break; + } + } + + if (divisor_is_single_limb) { + using CompilerUInt128 = unsigned __int128; + // Non-zero divisor guaranteed: denominator != 0 (checked) and high limbs are 0. + const CompilerUInt128 d = denominator.items[little(0)]; + + integer quotient = 0; + CompilerUInt128 remainder = 0; + for (int i = static_cast(item_count) - 1; i >= 0; --i) { + // remainder < d < 2^64, so cur < 2^128 and cur / d < 2^64 (fits one limb). + const CompilerUInt128 cur = + (remainder << base_bits) | + static_cast(numerator.items[little(i)]); + quotient.items[little(i)] = static_cast(cur / d); + remainder = cur % d; + } + + // Contract: return quotient, leave remainder (single limb) in numerator. + for (unsigned i = 0; i < item_count; ++i) { + numerator.items[little(i)] = 0; + } + numerator.items[little(0)] = static_cast(remainder); + return quotient; + } + } + + /// Fast path for a divisor that fits in two 64-bit limbs (128-bit): route it + /// through Knuth's Algorithm D (base 2^32). Reaching here means the divisor + /// did not fit a single limb, so it has >= 3 significant 32-bit half-limbs. + /// Divisors wider than 128 bits fall through to the generic binary long + /// division below. + if constexpr (item_count > 2) { + bool divisor_fits_two_limbs = true; + for (unsigned i = 2; i < item_count; ++i) { + if (denominator.items[little(i)] != 0) { + divisor_fits_two_limbs = false; + break; + } + } + if (divisor_fits_two_limbs) { + return divide_knuth(numerator, denominator); + } + } + integer x = 1; integer quotient = 0; diff --git a/be/test/core/wide_integer_test.cpp b/be/test/core/wide_integer_test.cpp index 6db61bcba4a7df..46bf7ec24ceb82 100644 --- a/be/test/core/wide_integer_test.cpp +++ b/be/test/core/wide_integer_test.cpp @@ -20,6 +20,9 @@ #include +#include +#include + #include "core/types.h" #include "core/uint128.h" @@ -194,4 +197,470 @@ TEST(WideInteger, Shift) { #endif } +TEST(WideInteger, SingleLimbDivisorFastPath) { + // A 256-bit dividend with all four 64-bit limbs populated. + const UInt256 n = (UInt256(0xFEDCBA9876543210ULL) << 192) | + (UInt256(0x1122334455667788ULL) << 128) | + (UInt256(0x99AABBCCDDEEFF00ULL) << 64) | UInt256(0x0123456789ABCDEFULL); + + // Divisors that all fit in a single 64-bit limb (exercise the fast path). + const UInt256 single_limb_divisors[] = { + UInt256(1ULL), + UInt256(2ULL), + UInt256(3ULL), + UInt256(7ULL), + UInt256(10ULL), + UInt256(1000000000ULL), + UInt256(1000000000000000000ULL), // 10^18 + UInt256(0x8000000000000000ULL), // 2^63 + UInt256(0xFFFFFFFFFFFFFFFFULL), // 2^64 - 1, the largest single limb + }; + for (const UInt256& d : single_limb_divisors) { + const UInt256 q = n / d; + const UInt256 r = n % d; + // q * d + r must reconstruct n, and the remainder must be strictly less than d. + ASSERT_EQ(q * d + r, n); + ASSERT_TRUE(r < d); + } +} + +TEST(WideInteger, MultiLimbDivisorBoundary) { + const UInt256 n = (UInt256(0xFEDCBA9876543210ULL) << 192) | + (UInt256(0x1122334455667788ULL) << 128) | + (UInt256(0x99AABBCCDDEEFF00ULL) << 64) | UInt256(0x0123456789ABCDEFULL); + + // Divisors that need two or more limbs must take the general path, not the fast path. + const UInt256 multi_limb_divisors[] = { + UInt256(1ULL) << 64, // 2^64: first value past single limb + (UInt256(1ULL) << 64) + UInt256(1ULL), // 2^64 + 1 + (UInt256(1ULL) << 100) + UInt256(12345ULL), + (UInt256(1ULL) << 192) + UInt256(0xDEADBEEFULL), + }; + for (const UInt256& d : multi_limb_divisors) { + const UInt256 q = n / d; + const UInt256 r = n % d; + ASSERT_EQ(q * d + r, n); + ASSERT_TRUE(r < d); + } +} + +TEST(WideInteger, SingleLimbKnownAnswers) { + // 10^38 / 10^19 == 10^19 exactly (10^19 fits in a single 64-bit limb). + const UInt256 p19 = UInt256(10000000000000000000ULL); // 10^19 + const UInt256 p38 = p19 * p19; // 10^38 + ASSERT_EQ(p38 / p19, p19); + ASSERT_EQ(p38 % p19, UInt256(0ULL)); + + // Exact division and non-zero remainder with a small divisor. + ASSERT_EQ(UInt256(1000ULL) / UInt256(7ULL), UInt256(142ULL)); + ASSERT_EQ(UInt256(1000ULL) % UInt256(7ULL), UInt256(6ULL)); + + // Dividend smaller than divisor -> quotient 0, remainder is the dividend. + ASSERT_EQ(UInt256(5ULL) / UInt256(9999999967ULL), UInt256(0ULL)); + ASSERT_EQ(UInt256(5ULL) % UInt256(9999999967ULL), UInt256(5ULL)); +} + +TEST(WideInteger, SingleLimbSignedDivision) { + // The fast path runs on the unsigned magnitudes; sign handling stays in the wrappers. + const Int256 n = (Int256(0x0011223344556677LL) << 128) | Int256(0x8899AABBCCDDEEFFLL); + const Int256 d = 1000000007; // prime, single limb + + ASSERT_EQ((-n) / d, -(n / d)); + ASSERT_EQ(n / (-d), -(n / d)); + ASSERT_EQ((-n) / (-d), n / d); + + // Truncation-toward-zero semantics for the remainder sign. + ASSERT_EQ((-n) % d, -(n % d)); + ASSERT_EQ(n % (-d), n % d); +} + +TEST(WideInteger, TwoLimbDivisorDifferential) { + // Cross-check the 128-bit-divisor Knuth path against the compiler's native + // unsigned __int128 division. Both operands fit in 128 bits so the quotient + // and remainder are exactly representable and comparable. + std::mt19937_64 rng(0xC0FFEE1234ULL); + auto rnd = [&rng]() { return rng(); }; + for (int iter = 0; iter < 20000; ++iter) { + const unsigned __int128 num = (static_cast(rnd()) << 64) | rnd(); + // Force a genuine 2-limb divisor: the high 64 bits must be non-zero. + unsigned __int128 den = (static_cast(rnd()) << 64) | rnd(); + if ((den >> 64) == 0) { + den |= (static_cast(1) << 64); + } + + const unsigned __int128 expected_q = num / den; + const unsigned __int128 expected_r = num % den; + + const UInt256 q = UInt256(num) / UInt256(den); + const UInt256 r = UInt256(num) % UInt256(den); + + ASSERT_EQ(static_cast(q), expected_q) << "iter " << iter; + ASSERT_EQ(static_cast(r), expected_r) << "iter " << iter; + } +} + +TEST(WideInteger, TwoLimbDivisorIdentity) { + // Full 256-bit dividends divided by random 2-limb (128-bit) divisors. No native + // oracle exists for 256-bit, so verify q*d + r == n with 0 <= r < d. q*d <= n < 2^256 + // so the multiplication does not overflow. + std::mt19937_64 rng(0xBADC0DE99ULL); + auto rnd = [&rng]() { return rng(); }; + auto make256 = [](uint64_t w3, uint64_t w2, uint64_t w1, uint64_t w0) { + return (UInt256(w3) << 192) | (UInt256(w2) << 128) | (UInt256(w1) << 64) | UInt256(w0); + }; + for (int iter = 0; iter < 20000; ++iter) { + const UInt256 n = make256(rnd(), rnd(), rnd(), rnd()); + uint64_t hi = rnd(); + if (hi == 0) { + hi = 1; // ensure the divisor spans two limbs (>= 2^64) + } + const UInt256 d = (UInt256(hi) << 64) | UInt256(rnd()); + + const UInt256 q = n / d; + const UInt256 r = n % d; + ASSERT_EQ(q * d + r, n) << "iter " << iter; + ASSERT_TRUE(r < d) << "iter " << iter; + } +} + +TEST(WideInteger, TwoLimbDivisorBoundaries) { + const UInt256 n = (UInt256(0xFEDCBA9876543210ULL) << 192) | + (UInt256(0x1122334455667788ULL) << 128) | + (UInt256(0x99AABBCCDDEEFF00ULL) << 64) | UInt256(0x0123456789ABCDEFULL); + + const UInt256 two_limb_divisors[] = { + UInt256(1ULL) << 64, // 2^64: smallest 2-limb divisor + (UInt256(1ULL) << 64) + UInt256(1ULL), // 2^64 + 1 + UInt256(1ULL) << 127, // high bit only (s == 0 after normalize) + (UInt256(1ULL) << 127) + UInt256(0x123456789AULL), + (UInt256(1ULL) << 128) - UInt256(1ULL), // 2^128 - 1: largest 2-limb divisor + (UInt256(10000000000000000000ULL) * UInt256(100ULL)), // 10^21, needs two limbs + }; + for (const UInt256& d : two_limb_divisors) { + const UInt256 q = n / d; + const UInt256 r = n % d; + ASSERT_EQ(q * d + r, n); + ASSERT_TRUE(r < d); + } + + // Known answer: 10^38 / 10^20 == 10^18, remainder 0 (10^20 is a 2-limb divisor). + const UInt256 p19 = UInt256(10000000000000000000ULL); + const UInt256 p38 = p19 * p19; + const UInt256 p20 = p19 * UInt256(10ULL); + const UInt256 p18 = UInt256(1000000000000000000ULL); + ASSERT_EQ(p38 / p20, p18); + ASSERT_EQ(p38 % p20, UInt256(0ULL)); +} + +TEST(WideInteger, WideDivisorRegression) { + // Divisors wider than 128 bits take the generic binary long-division fallback; + // confirm that path still satisfies the division identity. + std::mt19937_64 rng(0x5EED1234ULL); + auto rnd = [&rng]() { return rng(); }; + auto make256 = [](uint64_t w3, uint64_t w2, uint64_t w1, uint64_t w0) { + return (UInt256(w3) << 192) | (UInt256(w2) << 128) | (UInt256(w1) << 64) | UInt256(w0); + }; + for (int iter = 0; iter < 5000; ++iter) { + const UInt256 n = make256(rnd(), rnd(), rnd(), rnd()); + uint64_t w2 = rnd(); + if (w2 == 0) { + w2 = 1; // force the divisor past 128 bits (third limb non-zero) + } + const UInt256 d = make256(0, w2, rnd(), rnd()); + + const UInt256 q = n / d; + const UInt256 r = n % d; + ASSERT_EQ(q * d + r, n) << "iter " << iter; + ASSERT_TRUE(r < d) << "iter " << iter; + } +} + +TEST(WideInteger, TwoLimbSignedDivision) { + const Int256 n = (Int256(0x0011223344556677LL) << 192) | (Int256(0x18293A4B5C6D7E8FLL) << 128) | + (Int256(0x1122334455667788LL) << 64) | Int256(0x99AABBCCDDEEFF01LL); + const Int256 d = (Int256(0x0000000100000002LL) << 64) | Int256(0x0000000300000005LL); + + ASSERT_EQ((-n) / d, -(n / d)); + ASSERT_EQ(n / (-d), -(n / d)); + ASSERT_EQ((-n) / (-d), n / d); + ASSERT_EQ((-n) % d, -(n % d)); + ASSERT_EQ(n % (-d), n % d); +} + +// The 256-bit divide() fast path (Decimal256): when both operands fit in 128 +// bits it takes a single hardware __int128 divide instead of the shift-subtract +// software loop. These tests pin it to be bit-exact with the slow path and with +// native __int128, and confirm the slow path (high limbs set) still works. +TEST(WideInteger, Divide256FastPathFits128) { + // Values whose magnitudes fit in 128 bits -> fast path. Oracle = __int128. + // Cover both limbs (values > 2^64), exact and inexact quotients, r==0, a cases[] = { + {(unsigned __int128)1000000000000000000ULL, 7}, + {(hi << 64) | 0x0123456789abcdefULL, 1000000007ULL}, + {(hi << 60), (unsigned __int128)999983ULL}, + {(unsigned __int128)42, (unsigned __int128)100}, // a < b -> q=0, r=a + {(unsigned __int128)1000000, (unsigned __int128)1000}, // exact -> r=0 + {~(unsigned __int128)0, (unsigned __int128)3}, // max 128-bit numerator + {(hi << 64) | hi, hi}, + }; + for (auto [a128, b128] : cases) { + UInt256 a(a128); + UInt256 b(b128); + UInt256 q = a / b; + UInt256 r = a % b; + ASSERT_EQ(__uint128_t(q), a128 / b128); + ASSERT_EQ(__uint128_t(r), a128 % b128); + // Division identity within 256-bit arithmetic. + ASSERT_EQ(q * b + r, a); + } +} + +TEST(WideInteger, Divide256SlowPathHighLimbs) { + // Numerator needs bits above 128. The small (<=64-bit) divisors here now route + // through the layer-2 word-wise fast path; only the wide_divisor case below + // (divisor > 128 bits) genuinely exercises the shift-subtract software loop. + // Either way the results must reconstruct exactly. + UInt256 big = UInt256(12345678901234567890ULL) << 130; // bits set above 128 + big += UInt256(0xdeadbeefcafef00dULL); + + for (uint64_t d : + {uint64_t(1), uint64_t(2), uint64_t(1000000007ULL), uint64_t(9999999999999999999ULL)}) { + UInt256 divisor(d); + UInt256 q = big / divisor; + UInt256 r = big % divisor; + ASSERT_TRUE(r < divisor); + ASSERT_EQ(q * divisor + r, big); // reconstruct exactly + } + + // Divisor exceeds 128 bits -> both operands wide, so this is the true slow path. + UInt256 wide_divisor = UInt256(1000000009ULL) << 100; + UInt256 q = big / wide_divisor; + UInt256 r = big % wide_divisor; + ASSERT_TRUE(r < wide_divisor); + ASSERT_EQ(q * wide_divisor + r, big); +} + +TEST(WideInteger, Divide256WideNumeratorSmallDivisor) { + // The CEIL/ROUND `x / 10^k` hot shape: a full-width Decimal256 numerator with + // a small (<=64-bit) divisor. This exercises the word-wise long-division fast + // path (layer 2), which the "both fit 128" path (layer 1) does not cover. + // Build a numerator that genuinely needs bits above 128. + UInt256 x = UInt256(0xABCDEF0123456789ULL); + x <<= 64; + x += UInt256(0x1122334455667788ULL); + x <<= 64; + x += UInt256(0x99AABBCCDDEEFF00ULL); // ~192 significant bits + + // Powers of ten (the real rounding divisors) plus a couple of primes. + for (uint64_t d : {uint64_t(1), uint64_t(10), uint64_t(1000), uint64_t(1000000000ULL), + uint64_t(1000000000000000000ULL), uint64_t(1000000007ULL), + uint64_t(9999999999999999999ULL)}) { + UInt256 divisor(d); + UInt256 q = x / divisor; + UInt256 r = x % divisor; + ASSERT_TRUE(r < divisor); + ASSERT_EQ(q * divisor + r, x); // bit-exact reconstruction + } + + // The quotient of a wide value by a small divisor still needs high limbs, so + // this really is the wide path, not an accidental fit-128. + UInt256 q10 = x / UInt256(10); + ASSERT_TRUE(q10 > (UInt256(1) << 128)); +} + +TEST(WideInteger, Divide256SignedNegative) { + // Sign handling is applied by operator_slash around the unsigned divide; the + // fast path must not disturb it. Both operands fit 128 bits. + Int256 a = Int256(1000000000000000000LL) * Int256(1000000000LL); // 1e27, fits 128b + Int256 b = 999983; + ASSERT_EQ(a / b, Int256(1000000000000000000LL) * Int256(1000000000LL) / Int256(999983)); + ASSERT_EQ((-a) / b, -(a / b)); + ASSERT_EQ(a / (-b), -(a / b)); + ASSERT_EQ((-a) / (-b), a / b); + // Remainder sign follows the dividend (truncated division). + ASSERT_EQ((-a) % b, -(a % b)); + ASSERT_EQ(a % b + (a / b) * b, a); +} + +TEST(WideInteger, Divide256BoundaryAt128Bits) { + // Exactly 2^128 in the numerator forces a high limb -> slow path; 2^128 - 1 + // is the largest fast-path numerator. Both must be correct. + UInt256 two_128 = UInt256(1) << 128; + UInt256 max_128 = two_128 - 1; // all low 128 bits set, high limbs zero + UInt256 d(1000000007ULL); + + UInt256 q1 = max_128 / d, r1 = max_128 % d; + ASSERT_EQ(q1 * d + r1, max_128); + ASSERT_TRUE(r1 < d); + + UInt256 q2 = two_128 / d, r2 = two_128 % d; + ASSERT_EQ(q2 * d + r2, two_128); + ASSERT_TRUE(r2 < d); + + // The two results differ by exactly the extra unit. + ASSERT_EQ(two_128 - max_128, UInt256(1)); +} + +TEST(WideInteger, Divide256ByZeroThrows) { + // A zero denominator must still throw on the wide path (fast path skips it). + UInt256 a(123456789ULL); + UInt256 zero(0); + bool threw = false; + try { + UInt256 q = a / zero; + (void)q; + } catch (...) { + threw = true; + } + ASSERT_TRUE(threw); +} + +// Assemble a 256-bit value from four explicit 64-bit limbs (limb0 = least +// significant). Lets tests place bits in any limb, including the top one. +static UInt256 u256_from_limbs(uint64_t l3, uint64_t l2, uint64_t l1, uint64_t l0) { + UInt256 v(l3); + v <<= 64; + v += UInt256(l2); + v <<= 64; + v += UInt256(l1); + v <<= 64; + v += UInt256(l0); + return v; +} + +// Randomized layer-1 fast path ("both operands fit 128 bits") against a native +// __int128 oracle -- the strongest possible reference, hardware division. Every +// generated (a, b) has zero high limbs so it is guaranteed to take layer 1. +TEST(WideInteger, Divide256FastPathFits128Fuzz) { + std::mt19937_64 rng(0xD1A5C0DE12345678ULL); + std::uniform_int_distribution anybits; + + for (int iter = 0; iter < 50000; ++iter) { + unsigned __int128 a = (static_cast(anybits(rng)) << 64) | anybits(rng); + unsigned __int128 b = (static_cast(anybits(rng)) << 64) | anybits(rng); + if (b == 0) { + b = 1; + } + // Occasionally shrink one side so a < b, a fits 64, exact multiples, etc. + switch (iter & 7) { + case 0: + a &= 0xFFFFFFFFFFFFFFFFULL; + break; // a fits 64 bits + case 1: + b &= 0xFFFFFFFFULL; + if (b == 0) b = 1; + break; // small divisor + case 2: + a = b + (a % (b ? b : 1)); + break; // a slightly >= b + case 3: + b = a ? (a / 2 + 1) : 1; + break; // q around 1..2 + default: + break; + } + + UInt256 wa(a), wb(b); + UInt256 q = wa / wb; + UInt256 r = wa % wb; + ASSERT_EQ(__uint128_t(q), a / b) << "iter=" << iter; + ASSERT_EQ(__uint128_t(r), a % b) << "iter=" << iter; + } +} + +// Randomized layer-2 fast path (wide dividend / <=64-bit divisor) with an +// INDEPENDENT ground-truth oracle. Instead of computing q,r from the operation +// and checking q*d+r==x (which a wrapped-around wrong q could still satisfy), we +// pick q and r FIRST, build x = q*d + r under bounds that cannot overflow 2^256, +// then require the divide to recover exactly that q and r. This is the check the +// existing reconstruction tests cannot make. +TEST(WideInteger, Divide256WideBySmallGroundTruth) { + std::mt19937_64 rng(0xC0FFEE5EED9901ULL); + std::uniform_int_distribution anybits; + + const uint64_t edge_divisors[] = {1ULL, + 2ULL, + 10ULL, + 1000000000000000000ULL, // 10^18 + 10000000000000000000ULL, // 10^19 + (1ULL << 63), // 2^63 + 0xFFFFFFFFFFFFFFFFULL, // 2^64 - 1 + 1000000007ULL}; + const size_t edge_count = sizeof(edge_divisors) / sizeof(uint64_t); + + for (int iter = 0; iter < 50000; ++iter) { + // q uses limbs 0..2 (up to 192 bits). With d < 2^64, q*d < 2^256, so + // x = q*d + r (r < d < 2^64) never overflows -> division is exact. + UInt256 q = u256_from_limbs(0, anybits(rng), anybits(rng), anybits(rng)); + + uint64_t d = ((iter & 3) == 0) ? edge_divisors[(iter >> 2) % edge_count] : anybits(rng); + if (d == 0) { + d = 1; + } + uint64_t r = (d == 1) ? 0 : (anybits(rng) % d); // 0 <= r < d + + UInt256 dv(d); + UInt256 x = q * dv + UInt256(r); + + ASSERT_EQ(x / dv, q) << "iter=" << iter << " d=" << d; + ASSERT_EQ(x % dv, UInt256(r)) << "iter=" << iter << " d=" << d; + } + + // Also pin the all-limbs-set extreme: numerator = 2^256 - 1. We cannot build + // this from a known q, so we fall back to reconstruction + r < d as a smoke + // test for the largest possible dividend (the randomized cases above are the + // real ground-truth guarantee). + UInt256 maxv = ~UInt256(0); + for (uint64_t d : edge_divisors) { + UInt256 dv(d); + UInt256 q = maxv / dv; + UInt256 r = maxv % dv; + ASSERT_TRUE(r < dv) << "d=" << d; + ASSERT_EQ(q * dv + r, maxv) << "d=" << d; + } +} + +// Signed wide dividend / small divisor with an INDEPENDENT oracle for both sign +// and magnitude. Truncated division: q truncates toward zero, and the remainder +// takes the sign of the dividend. We construct a positive (q_abs, r_abs, d), then +// check all four sign combinations against values derived without calling divide. +TEST(WideInteger, Divide256SignedWideGroundTruth) { + std::mt19937_64 rng(0x5160EDA7A1234567ULL); + std::uniform_int_distribution anybits; + + for (int iter = 0; iter < 20000; ++iter) { + // Cap q_abs to ~189 bits so the signed magnitude X = q_abs*d + r_abs stays + // well below 2^255 for any d < 2^64 (no signed overflow), while still + // needing bits above 128 -> genuinely the wide path. + uint64_t l2 = anybits(rng) & ((1ULL << 61) - 1); + UInt256 q_abs_u = u256_from_limbs(0, l2, anybits(rng), anybits(rng)); + + uint64_t d = anybits(rng); + if (d < 2) { + d = 2; // keep room for a nonzero remainder + } + uint64_t r_abs = anybits(rng) % d; // 0 <= r_abs < d + + UInt256 x_u = q_abs_u * UInt256(d) + UInt256(r_abs); + Int256 X(x_u); + Int256 D(d); + Int256 Q(q_abs_u); + Int256 R(r_abs); + + // ++ : x/d + ASSERT_EQ(X / D, Q) << "iter=" << iter; + ASSERT_EQ(X % D, R) << "iter=" << iter; + // -+ : (-x)/d = -q, remainder follows dividend sign -> -r + ASSERT_EQ((-X) / D, -Q) << "iter=" << iter; + ASSERT_EQ((-X) % D, -R) << "iter=" << iter; + // +- : x/(-d) = -q, remainder still follows dividend sign -> +r + ASSERT_EQ(X / (-D), -Q) << "iter=" << iter; + ASSERT_EQ(X % (-D), R) << "iter=" << iter; + // -- : (-x)/(-d) = q, remainder -> -r + ASSERT_EQ((-X) / (-D), Q) << "iter=" << iter; + ASSERT_EQ((-X) % (-D), -R) << "iter=" << iter; + } +} + } // namespace doris