From ab76d46fac552467ebca002dfa127b7312b60ea7 Mon Sep 17 00:00:00 2001 From: cong-or Date: Tue, 28 Jul 2026 10:37:03 +0100 Subject: [PATCH 1/2] Add wide division: div_rem_wide and wrapping_div_wide Closes #1315 --- src/uint/div.rs | 262 ++++++++++++++++++++++++++++- src/uint/ref_type/div.rs | 346 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 607 insertions(+), 1 deletion(-) diff --git a/src/uint/div.rs b/src/uint/div.rs index 0cb654421..149043f9a 100644 --- a/src/uint/div.rs +++ b/src/uint/div.rs @@ -168,6 +168,126 @@ impl Uint { y } + /// Computes `(lo + hi * 2^Self::BITS) / rhs` for a double-width dividend, returning the + /// wrapped quotient and the remainder. + /// + /// The quotient of such a dividend may exceed `Self::BITS`; only its low `Self::BITS` bits are + /// returned (i.e. the quotient is reduced modulo `2^Self::BITS`). This is the quotient-tracking + /// counterpart of [`Uint::rem_wide`], and avoids widening the operands via + /// [`Concat`][`crate::Concat`], so it is available for any limb count. + /// + /// ### Usage: + /// ``` + /// use crypto_bigint::{U256, NonZero}; + /// + /// // dividend = 3 * 2^256 + 5, so dividing by 3 gives quotient 2^256 + 1, remainder 2 + /// let lo = U256::from(5u64); + /// let hi = U256::from(3u64); + /// let rhs = NonZero::new(U256::from(3u64)).unwrap(); + /// let (quo, rem) = U256::div_rem_wide((lo, hi), &rhs); + /// + /// // the true quotient 2^256 + 1 doesn't fit in 256 bits, so it wraps down to 1 + /// assert_eq!(quo, U256::ONE); + /// assert_eq!(rem, U256::from(2u64)); + /// ``` + #[inline] + #[must_use] + pub const fn div_rem_wide(lower_upper: (Self, Self), rhs: &NonZero) -> (Self, Self) { + let (mut lo, mut hi) = lower_upper; + let mut y = *rhs.as_ref(); + let mut quo = Self::ZERO; + UintRef::div_rem_wide( + (lo.as_mut_uint_ref(), hi.as_mut_uint_ref()), + y.as_mut_uint_ref(), + quo.as_mut_uint_ref(), + ); + (quo, y) + } + + /// Computes the wrapped quotient `(lo + hi * 2^Self::BITS) / rhs`, reduced modulo + /// `2^Self::BITS`. + /// + /// The quotient-only counterpart of [`Uint::rem_wide`]; see [`Uint::div_rem_wide`] for + /// details. + #[inline] + #[must_use] + pub const fn wrapping_div_wide(lower_upper: (Self, Self), rhs: &NonZero) -> Self { + Self::div_rem_wide(lower_upper, rhs).0 + } + + /// Exactly divides the double-width dividend `(lo, hi)` by `rhs`, returning the quotient if + /// the division is exact and [`CtOption::none()`] if `rhs` does not divide the dividend. + /// + /// When the division is exact the quotient is guaranteed to fit in `Self`, so no wrapping + /// occurs. This is the wide counterpart of [`Uint::div_exact`]. + /// + /// ### Usage: + /// ``` + /// use crypto_bigint::{U256, NonZero}; + /// + /// let rhs = NonZero::new(U256::from(3u64)).unwrap(); + /// + /// // 15 = 3 * 5 exactly + /// let quo = U256::div_wide_exact((U256::from(15u64), U256::ZERO), &rhs).unwrap(); + /// assert_eq!(quo, U256::from(5u64)); + /// + /// // 16 is not divisible by 3 + /// let not_exact = U256::div_wide_exact((U256::from(16u64), U256::ZERO), &rhs); + /// assert!(bool::from(not_exact.is_none())); + /// ``` + #[inline] + #[must_use] + pub const fn div_wide_exact(lower_upper: (Self, Self), rhs: &NonZero) -> CtOption { + let (quo, rem) = Self::div_rem_wide(lower_upper, rhs); + CtOption::new(quo, rem.is_zero()) + } + + /// Computes `(lo + hi * 2^Self::BITS) / rhs` for a double-width dividend, returning the + /// wrapped quotient and the remainder. + /// + /// This is variable-time only with respect to `rhs`. When used with a fixed `rhs`, it is + /// constant-time with respect to the dividend. See [`Uint::div_rem_wide`] for details. + #[inline] + #[must_use] + pub const fn div_rem_wide_vartime( + lower_upper: (Self, Self), + rhs: &NonZero, + ) -> (Self, Self) { + let (mut lo, mut hi) = lower_upper; + let mut y = *rhs.as_ref(); + let mut quo = Self::ZERO; + UintRef::div_rem_wide_vartime( + (lo.as_mut_uint_ref(), hi.as_mut_uint_ref()), + y.as_mut_uint_ref(), + quo.as_mut_uint_ref(), + ); + (quo, y) + } + + /// Computes the wrapped quotient `(lo + hi * 2^Self::BITS) / rhs`, reduced modulo + /// `2^Self::BITS`. + /// + /// This is variable-time only with respect to `rhs`. See [`Uint::wrapping_div_wide`]. + #[inline] + #[must_use] + pub const fn wrapping_div_wide_vartime(lower_upper: (Self, Self), rhs: &NonZero) -> Self { + Self::div_rem_wide_vartime(lower_upper, rhs).0 + } + + /// Exactly divides the double-width dividend `(lo, hi)` by `rhs`, returning the quotient if the + /// division is exact and [`CtOption::none()`] otherwise. + /// + /// This is variable-time only with respect to `rhs`. See [`Uint::div_wide_exact`]. + #[inline] + #[must_use] + pub const fn div_wide_exact_vartime( + lower_upper: (Self, Self), + rhs: &NonZero, + ) -> CtOption { + let (quo, rem) = Self::div_rem_wide_vartime(lower_upper, rhs); + CtOption::new(quo, rem.is_zero()) + } + /// Computes `self` % 2^k. Faster than reduce since its a power of 2. /// Limited to 2^16-1 since Uint doesn't support higher. /// @@ -544,7 +664,11 @@ mod tests { }; #[cfg(feature = "rand_core")] - use {crate::Random, chacha20::ChaCha8Rng, rand_core::Rng, rand_core::SeedableRng}; + use { + crate::{Random, U192, U384}, + chacha20::ChaCha8Rng, + rand_core::{Rng, SeedableRng}, + }; #[test] fn div_word() { @@ -998,4 +1122,140 @@ mod tests { }; assert_eq!(a.divide_x_by_y(), U1024::from(2707385u64)); } + + /// Check the wide-division methods (constant-time and variable-time) for a dividend + /// `lo + hi * 2^(L * Limb::BITS)` against the trusted `div_rem` reference, which divides the + /// same value widened into `Uint` (with `W == 2 * L`). + fn check(lo: Uint, hi: Uint, den: Uint) { + let nz = den.to_nz().unwrap(); + + // Reference: build the wide dividend and divide it with the trusted `div_rem`. + let wide: Uint = lo.concat_resize(&hi); + let (full_q, full_r) = wide.div_rem(&den.resize::().to_nz().unwrap()); + let exp_q = full_q.resize::(); + let exp_r = full_r.resize::(); + let exact = exp_r == Uint::::ZERO; + + // Both the constant-time and variable-time paths must match the reference. + for (q, r) in [ + Uint::::div_rem_wide((lo, hi), &nz), + Uint::::div_rem_wide_vartime((lo, hi), &nz), + ] { + assert_eq!(q, exp_q, "div_rem_wide quotient: ({lo}, {hi}) / {den}"); + assert_eq!(r, exp_r, "div_rem_wide remainder: ({lo}, {hi}) / {den}"); + } + assert_eq!(Uint::::wrapping_div_wide((lo, hi), &nz), exp_q); + assert_eq!(Uint::::wrapping_div_wide_vartime((lo, hi), &nz), exp_q); + + for maybe_quo in [ + Uint::::div_wide_exact((lo, hi), &nz), + Uint::::div_wide_exact_vartime((lo, hi), &nz), + ] { + assert_eq!(bool::from(maybe_quo.is_some()), exact); + if exact { + assert_eq!(maybe_quo.unwrap(), exp_q); + } + } + } + + #[test] + fn div_rem_wide_edge() { + let two_127 = U128::from_be_hex("80000000000000000000000000000000"); // 2^127 + + // Boundary cases, each checked against the concat + `div_rem` reference. + for (lo, hi, den) in [ + (U128::ZERO, U128::ZERO, U128::from(7u64)), // both halves zero + (U128::from(100u64), U128::ZERO, U128::from(7u64)), // high half zero + (U128::ZERO, U128::ONE, U128::ONE), // divisor 1: quotient wraps + (U128::MAX, U128::MAX, U128::MAX), // 2^128 + 1 wraps to 1 + (U128::MAX, U128::ZERO, U128::from(2u64)), // half-word divisor + (two_127, U128::ONE, U128::from(3u64)), // 3 * 2^127, single-word divisor + ] { + check::<{ U128::LIMBS }, { U256::LIMBS }>(lo, hi, den); + } + // Single-word operands (LIMBS == 1). + for (lo, hi, den) in [ + (U64::MAX, U64::MAX, U64::MAX), + (U64::from(5u64), U64::from(9u64), U64::from(4u64)), + ] { + check::<{ U64::LIMBS }, { U128::LIMBS }>(lo, hi, den); + } + + // A couple of hand-computed values for extra confidence. + // (2^256 - 1) / (2^128 - 1) = 2^128 + 1, wrapped mod 2^128 = 1, remainder 0. + let (q, r) = U128::div_rem_wide((U128::MAX, U128::MAX), &U128::MAX.to_nz().unwrap()); + assert_eq!(q, U128::ONE); + assert_eq!(r, U128::ZERO); + + // 3 * 2^127 = 2^128 + 2^127 => hi = 1, lo = 2^127; the exact quotient is 2^127. + let three = U128::from(3u64).to_nz().unwrap(); + let (q, r) = U128::div_rem_wide((two_127, U128::ONE), &three); + assert_eq!(q, two_127); + assert_eq!(r, U128::ZERO); + assert_eq!( + U128::div_wide_exact((two_127, U128::ONE), &three).unwrap(), + two_127 + ); + } + + /// Force a randomly generated value to be non-zero so it is a valid divisor. + #[cfg(feature = "rand_core")] + fn nz(v: Uint) -> Uint { + if bool::from(v.is_zero()) { + Uint::ONE + } else { + v + } + } + + #[cfg(feature = "rand_core")] + #[test] + fn div_rem_wide_vs_concat() { + let mut rng = ChaCha8Rng::from_seed([9u8; 32]); + for _ in 0..300 { + // Single-word path (U64 operands). + let lo64 = U64::random_from_rng(&mut rng); + let hi64 = U64::random_from_rng(&mut rng); + check::<{ U64::LIMBS }, { U128::LIMBS }>( + lo64, + hi64, + nz(U64::random_from_rng(&mut rng)), + ); + + // Multi-limb path (U128): a full-width divisor, then a single-word-value divisor + // (ywords == 1, which hits the div2by1 correction). + let lo = U128::random_from_rng(&mut rng); + let hi = U128::random_from_rng(&mut rng); + check::<{ U128::LIMBS }, { U256::LIMBS }>(lo, hi, nz(U128::random_from_rng(&mut rng))); + let d_small = U64::random_from_rng(&mut rng).resize::<{ U128::LIMBS }>(); + check::<{ U128::LIMBS }, { U256::LIMBS }>(lo, hi, nz(d_small)); + + // Odd limb count (U192): a full-width divisor (>= 3 significant limbs, so the tail + // `done` masking / vartime early-break fires) and a narrower divisor. + let lo192 = U192::random_from_rng(&mut rng); + let hi192 = U192::random_from_rng(&mut rng); + check::<{ U192::LIMBS }, { U384::LIMBS }>( + lo192, + hi192, + nz(U192::random_from_rng(&mut rng)), + ); + let d192_narrow = U128::random_from_rng(&mut rng).resize::<{ U192::LIMBS }>(); + check::<{ U192::LIMBS }, { U384::LIMBS }>(lo192, hi192, nz(d192_narrow)); + + // Wider operands (U256): a full-width divisor, a single-word-value divisor, and one + // with exactly LIMBS-1 significant limbs (fires the vartime early-break while the + // divisor is still narrower than the dividend's high half). + let lo256 = U256::random_from_rng(&mut rng); + let hi256 = U256::random_from_rng(&mut rng); + check::<{ U256::LIMBS }, { U512::LIMBS }>( + lo256, + hi256, + nz(U256::random_from_rng(&mut rng)), + ); + let d256_narrow = U128::random_from_rng(&mut rng).resize::<{ U256::LIMBS }>(); + check::<{ U256::LIMBS }, { U512::LIMBS }>(lo256, hi256, nz(d256_narrow)); + let d256_3limb = U192::random_from_rng(&mut rng).resize::<{ U256::LIMBS }>(); + check::<{ U256::LIMBS }, { U512::LIMBS }>(lo256, hi256, nz(d256_3limb)); + } + } } diff --git a/src/uint/ref_type/div.rs b/src/uint/ref_type/div.rs index 2205b098d..fd85c5a64 100644 --- a/src/uint/ref_type/div.rs +++ b/src/uint/ref_type/div.rs @@ -204,6 +204,352 @@ impl UintRef { y.shr_assign_limb_vartime(lshift); } + /// Computes `x_lower_upper` / `rhs`, returning the wrapped quotient in `quo` and the + /// remainder in `rhs`. + /// + /// The `x_lower_upper` tuple represents a wide (double-width) dividend `x_lo + x_hi * B`, + /// where `B = 2^(x_lo.bits_precision())`. The size of `x_lower_upper.1` and of `quo` must each + /// be at least as large as `rhs`. `x_lower_upper` is left in an indeterminate state. + /// + /// The true quotient may be up to twice the width of `rhs`; only its low `quo.nlimbs()` limbs + /// are retained (i.e. the quotient is reduced modulo `2^(quo.nlimbs() * Limb::BITS)`). This is + /// the quotient-tracking counterpart of [`UintRef::rem_wide`]. + /// + /// # Panics + /// If the divisor is zero. + #[inline(always)] + pub(crate) const fn div_rem_wide( + x_lower_upper: (&mut Self, &mut Self), + rhs: &mut Self, + quo: &mut Self, + ) { + let (x_lo, x) = x_lower_upper; + let y = rhs; + + // Short circuit for single-word divisor (only reachable when the operands are one limb + // wide, since `div3by2` in the main path requires a two-limb divisor). + if y.nlimbs() == 1 { + let reciprocal = Reciprocal::new(y.limbs[0].to_nz().expect_copied("zero divisor")); + let sh = reciprocal.shift(); + + // Left-shift the wide dividend so that the divisor is normalized (high bit set). + let lo_carry = x_lo.shl_assign_limb(sh); + let mut hi = x.shl_assign_limb(sh); + x.limbs[0] = x.limbs[0].bitor(lo_carry); + + // Long division by a single limb, most-significant limb first. The high quotient limb + // overflows the wrapped result and is discarded; the low limb is the wrapped quotient. + (x.limbs[0].0, hi.0) = div2by1(x.limbs[0].0, hi.0, &reciprocal); + (x_lo.limbs[0].0, hi.0) = div2by1(x_lo.limbs[0].0, hi.0, &reciprocal); + + quo.limbs[0] = x_lo.limbs[0]; + y.limbs[0] = hi.shr(sh); + return; + } + + // Compute the size of the divisor + let ybits = y.bits(); + assert!(ybits > 0, "zero divisor"); + let ywords = ybits.div_ceil(Limb::BITS); + + // Shift the entire divisor such that the high bit is set + let yz = y.bits_precision() - ybits; + y.unbounded_shl_assign(yz); + + // Shift the dividend to align the words + let lshift = yz & (Limb::BITS - 1); + let x_lo_carry = x_lo.shl_assign_limb(lshift); + let x_hi = x.shl_assign_limb(lshift); + x.limbs[0] = x.limbs[0].bitor(x_lo_carry); + + // Perform the core division algorithm + Self::div_rem_wide_shifted((x_lo, x), x_hi, y, ywords, quo); + + // Unshift the remainder from the earlier adjustment + y.shr_assign_limb(lshift); + } + + /// Computes `x_lower_upper` / `rhs`, returning the wrapped quotient in `quo` and the + /// remainder in `rhs`. + /// + /// This function operates in variable-time with respect to `rhs`. For a fixed divisor, it + /// operates in constant-time. + /// + /// The `x_lower_upper` tuple represents a wide (double-width) dividend. The size of + /// `x_lower_upper.1` and of `quo` must each be at least as large as `rhs`. `x_lower_upper` is + /// left in an indeterminate state. See [`UintRef::div_rem_wide`] for the wrapping semantics. + /// + /// # Panics + /// If the divisor is zero. + #[inline(always)] + pub(crate) const fn div_rem_wide_vartime( + x_lower_upper: (&mut Self, &mut Self), + rhs: &mut Self, + quo: &mut Self, + ) { + let (x_lo, x) = x_lower_upper; + let xsize = x.nlimbs(); + let ysize = bitlen::to_limbs(rhs.bits_vartime()); + let y = rhs.leading_mut(ysize); + + match (xsize, ysize) { + (_, 0) => panic!("zero divisor"), + (0, _) => { + // Empty dividend: both quotient and remainder are zero. + y.fill(Limb::ZERO); + quo.fill(Limb::ZERO); + return; + } + (_, 1) => { + // Single-word divisor: long division by one limb, most-significant limb first, + // keeping the quotient. The high half's quotient overflows the wrapped result and + // is discarded; the low half's quotient is what we keep. + let reciprocal = Reciprocal::new(y.limbs[0].to_nz().expect_copied("zero divisor")); + let sh = reciprocal.shift(); + + // Left-shift the wide dividend so the divisor is normalized (high bit set). + let lo_carry = x_lo.shl_assign_limb_vartime(sh); + let mut hi = x.shl_assign_limb_vartime(sh); + x.limbs[0] = x.limbs[0].bitor(lo_carry); + + let mut j = xsize; + while j > 0 { + j -= 1; + (x.limbs[j].0, hi.0) = div2by1(x.limbs[j].0, hi.0, &reciprocal); + } + let mut j = x_lo.nlimbs(); + while j > 0 { + j -= 1; + (x_lo.limbs[j].0, hi.0) = div2by1(x_lo.limbs[j].0, hi.0, &reciprocal); + } + + quo.copy_from(x_lo); + y.fill(Limb::ZERO); + y.limbs[0] = hi.shr(sh); + return; + } + _ if ysize > xsize => { + panic!("divisor too large"); + } + _ => (), + } + + let lshift = y.limbs[ysize - 1].leading_zeros(); + + // Shift divisor such that it has no leading zeros + // This means that div2by1 requires no extra shifts, and ensures that the high word >= b/2 + y.shl_assign_limb_vartime(lshift); + + // Shift the dividend to align the words + let x_lo_carry = x_lo.shl_assign_limb_vartime(lshift); + let mut x_hi = x.shl_assign_limb_vartime(lshift); + x.limbs[0] = x.limbs[0].bitor(x_lo_carry); + + // Calculate a reciprocal from the highest word of the divisor + let reciprocal = Reciprocal::new(y.limbs[ysize - 1].to_nz().expect_copied("zero divisor")); + + // Perform the core division algorithm + x_hi = Self::div_rem_wide_large_shifted::( + (x_lo, x), + x_hi, + y, + #[allow(clippy::cast_possible_truncation, reason = "TODO")] + { + ysize as u32 + }, + reciprocal, + quo, + ); + + // Copy the remainder to the divisor + y.leading_mut(ysize - 1).copy_from(x.leading(ysize - 1)); + y.limbs[ysize - 1] = x_hi; + + // Unshift the remainder from the earlier adjustment + y.shr_assign_limb_vartime(lshift); + } + + /// Conditionally shift the limbs one position toward the most-significant end, dropping the + /// top limb and inserting `limb` at the least-significant position. A no-op when `shift` is + /// falsy. + /// + /// Used as a quotient accumulator: feeding quotient limbs most-significant first retains the + /// low `self.nlimbs()` limbs of the full-width quotient. + #[inline(always)] + const fn shift_in_limb(&mut self, limb: Limb, shift: Choice) { + let mut j = self.nlimbs(); + while j > 1 { + j -= 1; + self.limbs[j] = Limb::select(self.limbs[j], self.limbs[j - 1], shift); + } + if self.nlimbs() > 0 { + self.limbs[0] = Limb::select(self.limbs[0], limb, shift); + } + } + + /// Perform in-place wide division for a pre-shifted dividend and divisor, tracking both the + /// quotient and the remainder. + /// + /// The dividend and divisor must be left-shifted such that the high bit of the divisor is set, + /// and `x_hi` holds the top bits of the (high half of the) dividend. + /// + /// The wrapped quotient is written to `quo` and the shifted remainder to `y` (the latter must + /// be unshifted by the caller). `x` is left in an indeterminate state. + #[inline(always)] + #[allow(clippy::cast_possible_truncation)] + const fn div_rem_wide_shifted( + x: (&mut Self, &mut Self), + mut x_hi: Limb, + y: &mut Self, + ywords: u32, + quo: &mut Self, + ) { + let (x_lo, x) = x; + let ysize = y.nlimbs(); + + // Calculate a reciprocal from the highest word of the divisor + let reciprocal = Reciprocal::new(y.limbs[ysize - 1].to_nz().expect_copied("zero divisor")); + debug_assert!(reciprocal.shift() == 0); + + // Perform the core division algorithm + x_hi = + Self::div_rem_wide_large_shifted::((x_lo, x), x_hi, y, ywords, reciprocal, quo); + + // Calculate quotient and remainder for the case where the divisor is a single word. + let limb_div = Choice::from_u32_eq(1, ywords); + // Note that `div2by1()` will panic if `x_hi >= reciprocal.divisor_normalized`, + // but this can only be the case if `limb_div` is falsy, in which case we discard + // the result anyway, so we conditionally set `x_hi` to zero for this branch. + let x_hi_adjusted = Limb::select(Limb::ZERO, x_hi, limb_div); + let (quo2, rem2) = div2by1(x.limbs[0].0, x_hi_adjusted.0, &reciprocal); + + // For a single-word divisor the main loop never computes the least-significant quotient + // limb; inject it here by shifting the quotient up one limb and storing `quo2`. + quo.shift_in_limb(Limb(quo2), limb_div); + + // Copy out the low limb of the remainder + y.limbs[0] = Limb::select(x.limbs[0], Limb(rem2), limb_div); + + // Copy the remainder to divisor + let mut i = 1; + while i < ysize { + y.limbs[i] = Limb::select( + Limb::ZERO, + x.limbs[i], + Choice::from_u32_lt(i as u32, ywords), + ); + y.limbs[i] = Limb::select(y.limbs[i], x_hi, Choice::from_u32_eq(i as u32, ywords - 1)); + i += 1; + } + } + + /// Computes `x` / `y` for a "large" divisor (>1 limbs), returning the shifted remainder in + /// `x.1` and the wrapped quotient in `quo`. + /// + /// Mirrors [`UintRef::rem_wide_large_shifted`], additionally capturing the (borrow-corrected) + /// quotient word at each step. The dividend and divisor must be left-shifted such that the + /// high bit of the divisor is set, and `x_hi` holds the top bits of the dividend. + #[inline(always)] + #[allow(clippy::cast_possible_truncation)] + const fn div_rem_wide_large_shifted( + x: (&Self, &mut Self), + mut x_hi: Limb, + y: &Self, + ywords: u32, + reciprocal: Reciprocal, + quo: &mut Self, + ) -> Limb { + assert!( + y.nlimbs() <= x.1.nlimbs(), + "invalid input sizes for div_rem_wide_large_shifted" + ); + + let (x_lo, x) = x; + let xsize = x.nlimbs(); + let ysize = y.nlimbs(); + let mut extra_limbs = x_lo.nlimbs(); + + let mut xi = xsize - 1; + let mut x_xi = x.limbs[xi]; + let mut i; + let mut carry; + + // Compute the adjusted reciprocal + let v = reciprocal.reciprocal_3by2(y.limbs[ysize - 2].0, y.limbs[ysize - 1].0); + + while xi > 0 { + // Divide high dividend words by the high divisor word to estimate the quotient word + let (mut quotient_word, _) = div3by2( + (x.limbs[xi - 1].0, x_xi.0, x_hi.0), + (y.limbs[ysize - 2].0, y.limbs[ysize - 1].0), + v, + ); + + // This loop is a no-op once xi is smaller than the number of words in the divisor. + // In variable-time mode we can stop as soon as that happens. + let done = Choice::from_u32_lt(xi as u32, ywords - 1); + if VARTIME && done.to_bool_vartime() { + break; + } + quotient_word = word::select(quotient_word, 0, done); + + // Subtract q*divisor from the dividend + let borrow = { + carry = Limb::ZERO; + let mut borrow = Limb::ZERO; + let mut tmp; + i = (xi + 1).saturating_sub(ysize); + while i <= xi { + (tmp, carry) = y.limbs[ysize + i - xi - 1].carrying_mul_add( + Limb(quotient_word), + carry, + Limb::ZERO, + ); + (x.limbs[i], borrow) = x.limbs[i].borrowing_sub(tmp, borrow); + i += 1; + } + (_, borrow) = x_hi.borrowing_sub(carry, borrow); + borrow + }; + + // If the subtraction borrowed, then decrement quo and add back the divisor. + // The probability of this being needed is very low, about 2/(Limb::MAX+1) + quotient_word = { + carry = Limb::ZERO; + i = (xi + 1).saturating_sub(ysize); + while i <= xi { + (x.limbs[i], carry) = + x.limbs[i].carrying_add(y.limbs[ysize + i - xi - 1].bitand(borrow), carry); + i += 1; + } + quotient_word.saturating_sub(borrow.0 & 1) + }; + + // Capture the corrected quotient word (most-significant first). The `done` iterations + // at the tail contribute no quotient word and must not shift the accumulator. + quo.shift_in_limb(Limb(quotient_word), done.not()); + + // If we have lower limbs remaining, shift the dividend words one word left + if extra_limbs > 0 { + x_hi = x.limbs[xi]; + x_xi = x.limbs[xi - 1]; + extra_limbs -= 1; + i = xi; + while i > 0 { + x.limbs[i] = x.limbs[i - 1]; + i -= 1; + } + x.limbs[0] = x_lo.limbs[extra_limbs]; + } else { + x_hi = Limb::select(x.limbs[xi], x_hi, done); + x_xi = Limb::select(x.limbs[xi - 1], x_xi, done); + xi -= 1; + } + } + + x_hi + } + /// Perform in-place division (`self` / `y`) for a pre-shifted dividend and divisor. /// /// The dividend and divisor must be left-shifted such that the high bit of the divisor From 0ac5d1ca7cfbeaaa0f767eb63a107741112a9eeb Mon Sep 17 00:00:00 2001 From: cong-or Date: Wed, 5 Aug 2026 20:15:38 +0100 Subject: [PATCH 2/2] Address review feedback on wide division Make the wide-division methods clearer and safer to use: - Rename div_rem_wide to wrapping_div_rem_wide so the name warns that the quotient can be cut short when it is too big to fit. - Also hand back a yes/no flag saying whether the quotient fit without being cut short. - Fix exact division: a zero remainder on its own is not enough, the quotient has to fit too. For example 2^128 divided by 1 now correctly reports "no exact answer" instead of returning a wrong value. - Reuse the existing limb-shift helper instead of a near-copy of it. - Tidy the tests: clearer name for the checker and simpler random divisors. --- src/uint/div.rs | 203 +++++++++++++++++++++++++-------------- src/uint/ref_type/div.rs | 60 ++++++++---- 2 files changed, 169 insertions(+), 94 deletions(-) diff --git a/src/uint/div.rs b/src/uint/div.rs index 149043f9a..aab49a363 100644 --- a/src/uint/div.rs +++ b/src/uint/div.rs @@ -2,8 +2,8 @@ use super::div_limb::Reciprocal; use crate::{ - CheckedDiv, CtOption, Div, DivAssign, DivRemLimb, DivVartime, Limb, NonZero, Rem, RemAssign, - RemLimb, RemMixed, ToUnsigned, Uint, UintRef, Unsigned, Wrapping, + CheckedDiv, Choice, CtOption, Div, DivAssign, DivRemLimb, DivVartime, Limb, NonZero, Rem, + RemAssign, RemLimb, RemMixed, ToUnsigned, Uint, UintRef, Unsigned, Wrapping, }; impl Uint { @@ -169,11 +169,14 @@ impl Uint { } /// Computes `(lo + hi * 2^Self::BITS) / rhs` for a double-width dividend, returning the - /// wrapped quotient and the remainder. + /// wrapped quotient, the remainder, and a [`Choice`] that is truthy when the quotient fit in + /// `Self` without truncation. /// /// The quotient of such a dividend may exceed `Self::BITS`; only its low `Self::BITS` bits are - /// returned (i.e. the quotient is reduced modulo `2^Self::BITS`). This is the quotient-tracking - /// counterpart of [`Uint::rem_wide`], and avoids widening the operands via + /// returned (i.e. the quotient is reduced modulo `2^Self::BITS`), which is why the name is + /// prefixed with `wrapping`. The returned [`Choice`] is truthy exactly when no wrapping + /// occurred, i.e. when the high half of the dividend is less than `rhs`. This is the + /// quotient-tracking counterpart of [`Uint::rem_wide`], and avoids widening the operands via /// [`Concat`][`crate::Concat`], so it is available for any limb count. /// /// ### Usage: @@ -184,42 +187,49 @@ impl Uint { /// let lo = U256::from(5u64); /// let hi = U256::from(3u64); /// let rhs = NonZero::new(U256::from(3u64)).unwrap(); - /// let (quo, rem) = U256::div_rem_wide((lo, hi), &rhs); + /// let (quo, rem, fits) = U256::wrapping_div_rem_wide((lo, hi), &rhs); /// /// // the true quotient 2^256 + 1 doesn't fit in 256 bits, so it wraps down to 1 /// assert_eq!(quo, U256::ONE); /// assert_eq!(rem, U256::from(2u64)); + /// assert!(!bool::from(fits)); /// ``` #[inline] #[must_use] - pub const fn div_rem_wide(lower_upper: (Self, Self), rhs: &NonZero) -> (Self, Self) { + pub const fn wrapping_div_rem_wide( + lower_upper: (Self, Self), + rhs: &NonZero, + ) -> (Self, Self, Choice) { let (mut lo, mut hi) = lower_upper; let mut y = *rhs.as_ref(); let mut quo = Self::ZERO; - UintRef::div_rem_wide( + let fits = UintRef::wrapping_div_rem_wide( (lo.as_mut_uint_ref(), hi.as_mut_uint_ref()), y.as_mut_uint_ref(), quo.as_mut_uint_ref(), ); - (quo, y) + (quo, y, fits) } /// Computes the wrapped quotient `(lo + hi * 2^Self::BITS) / rhs`, reduced modulo /// `2^Self::BITS`. /// - /// The quotient-only counterpart of [`Uint::rem_wide`]; see [`Uint::div_rem_wide`] for - /// details. + /// The quotient-only counterpart of [`Uint::rem_wide`]; see [`Uint::wrapping_div_rem_wide`] + /// for details. #[inline] #[must_use] pub const fn wrapping_div_wide(lower_upper: (Self, Self), rhs: &NonZero) -> Self { - Self::div_rem_wide(lower_upper, rhs).0 + Self::wrapping_div_rem_wide(lower_upper, rhs).0 } - /// Exactly divides the double-width dividend `(lo, hi)` by `rhs`, returning the quotient if - /// the division is exact and [`CtOption::none()`] if `rhs` does not divide the dividend. + /// Exactly divides the double-width dividend `(lo, hi)` by `rhs`, returning the quotient in a + /// [`CtOption`] that is [`none`][`CtOption::none()`] unless the division is exact *and* the + /// quotient fits in `Self`. /// - /// When the division is exact the quotient is guaranteed to fit in `Self`, so no wrapping - /// occurs. This is the wide counterpart of [`Uint::div_exact`]. + /// The quotient of a double-width dividend may exceed `Self::BITS` even when the division is + /// exact (e.g. `2^Self::BITS / 1`), so a zero remainder alone is not sufficient: the result is + /// only present when the true quotient is also representable in `Self`. This is the wide + /// counterpart of [`Uint::div_exact`]. /// /// ### Usage: /// ``` @@ -234,34 +244,40 @@ impl Uint { /// // 16 is not divisible by 3 /// let not_exact = U256::div_wide_exact((U256::from(16u64), U256::ZERO), &rhs); /// assert!(bool::from(not_exact.is_none())); + /// + /// // 2^256 is divisible by 1, but the quotient 2^256 does not fit in `U256` + /// let one = NonZero::new(U256::ONE).unwrap(); + /// let overflows = U256::div_wide_exact((U256::ZERO, U256::ONE), &one); + /// assert!(bool::from(overflows.is_none())); /// ``` #[inline] #[must_use] pub const fn div_wide_exact(lower_upper: (Self, Self), rhs: &NonZero) -> CtOption { - let (quo, rem) = Self::div_rem_wide(lower_upper, rhs); - CtOption::new(quo, rem.is_zero()) + let (quo, rem, fits) = Self::wrapping_div_rem_wide(lower_upper, rhs); + CtOption::new(quo, rem.is_zero().and(fits)) } /// Computes `(lo + hi * 2^Self::BITS) / rhs` for a double-width dividend, returning the - /// wrapped quotient and the remainder. + /// wrapped quotient, the remainder, and a [`Choice`] that is truthy when the quotient fit in + /// `Self` without truncation. /// /// This is variable-time only with respect to `rhs`. When used with a fixed `rhs`, it is - /// constant-time with respect to the dividend. See [`Uint::div_rem_wide`] for details. + /// constant-time with respect to the dividend. See [`Uint::wrapping_div_rem_wide`] for details. #[inline] #[must_use] - pub const fn div_rem_wide_vartime( + pub const fn wrapping_div_rem_wide_vartime( lower_upper: (Self, Self), rhs: &NonZero, - ) -> (Self, Self) { + ) -> (Self, Self, Choice) { let (mut lo, mut hi) = lower_upper; let mut y = *rhs.as_ref(); let mut quo = Self::ZERO; - UintRef::div_rem_wide_vartime( + let fits = UintRef::wrapping_div_rem_wide_vartime( (lo.as_mut_uint_ref(), hi.as_mut_uint_ref()), y.as_mut_uint_ref(), quo.as_mut_uint_ref(), ); - (quo, y) + (quo, y, fits) } /// Computes the wrapped quotient `(lo + hi * 2^Self::BITS) / rhs`, reduced modulo @@ -271,11 +287,12 @@ impl Uint { #[inline] #[must_use] pub const fn wrapping_div_wide_vartime(lower_upper: (Self, Self), rhs: &NonZero) -> Self { - Self::div_rem_wide_vartime(lower_upper, rhs).0 + Self::wrapping_div_rem_wide_vartime(lower_upper, rhs).0 } - /// Exactly divides the double-width dividend `(lo, hi)` by `rhs`, returning the quotient if the - /// division is exact and [`CtOption::none()`] otherwise. + /// Exactly divides the double-width dividend `(lo, hi)` by `rhs`, returning the quotient in a + /// [`CtOption`] that is [`none`][`CtOption::none()`] unless the division is exact *and* the + /// quotient fits in `Self`. /// /// This is variable-time only with respect to `rhs`. See [`Uint::div_wide_exact`]. #[inline] @@ -284,8 +301,8 @@ impl Uint { lower_upper: (Self, Self), rhs: &NonZero, ) -> CtOption { - let (quo, rem) = Self::div_rem_wide_vartime(lower_upper, rhs); - CtOption::new(quo, rem.is_zero()) + let (quo, rem, fits) = Self::wrapping_div_rem_wide_vartime(lower_upper, rhs); + CtOption::new(quo, rem.is_zero().and(fits)) } /// Computes `self` % 2^k. Faster than reduce since its a power of 2. @@ -1126,33 +1143,47 @@ mod tests { /// Check the wide-division methods (constant-time and variable-time) for a dividend /// `lo + hi * 2^(L * Limb::BITS)` against the trusted `div_rem` reference, which divides the /// same value widened into `Uint` (with `W == 2 * L`). - fn check(lo: Uint, hi: Uint, den: Uint) { - let nz = den.to_nz().unwrap(); + fn check_wide_division( + lo: Uint, + hi: Uint, + den: NonZero>, + ) { + let den_uint = *den.as_ref(); // Reference: build the wide dividend and divide it with the trusted `div_rem`. let wide: Uint = lo.concat_resize(&hi); - let (full_q, full_r) = wide.div_rem(&den.resize::().to_nz().unwrap()); + let (full_q, full_r) = wide.div_rem(&den_uint.resize::().to_nz().unwrap()); let exp_q = full_q.resize::(); let exp_r = full_r.resize::(); let exact = exp_r == Uint::::ZERO; + // The quotient fits in `L` limbs exactly when its high half (the limbs dropped by the + // wrapping division) is zero. + let fits = full_q == exp_q.resize::(); // Both the constant-time and variable-time paths must match the reference. - for (q, r) in [ - Uint::::div_rem_wide((lo, hi), &nz), - Uint::::div_rem_wide_vartime((lo, hi), &nz), + for (q, r, quo_fits) in [ + Uint::::wrapping_div_rem_wide((lo, hi), &den), + Uint::::wrapping_div_rem_wide_vartime((lo, hi), &den), ] { - assert_eq!(q, exp_q, "div_rem_wide quotient: ({lo}, {hi}) / {den}"); - assert_eq!(r, exp_r, "div_rem_wide remainder: ({lo}, {hi}) / {den}"); + assert_eq!(q, exp_q, "quotient: ({lo}, {hi}) / {den_uint}"); + assert_eq!(r, exp_r, "remainder: ({lo}, {hi}) / {den_uint}"); + assert_eq!( + bool::from(quo_fits), + fits, + "fits: ({lo}, {hi}) / {den_uint}" + ); } - assert_eq!(Uint::::wrapping_div_wide((lo, hi), &nz), exp_q); - assert_eq!(Uint::::wrapping_div_wide_vartime((lo, hi), &nz), exp_q); + assert_eq!(Uint::::wrapping_div_wide((lo, hi), &den), exp_q); + assert_eq!(Uint::::wrapping_div_wide_vartime((lo, hi), &den), exp_q); + // `div_wide_exact` yields the quotient only when the division is exact *and* it fits. + let exact_and_fits = exact && fits; for maybe_quo in [ - Uint::::div_wide_exact((lo, hi), &nz), - Uint::::div_wide_exact_vartime((lo, hi), &nz), + Uint::::div_wide_exact((lo, hi), &den), + Uint::::div_wide_exact_vartime((lo, hi), &den), ] { - assert_eq!(bool::from(maybe_quo.is_some()), exact); - if exact { + assert_eq!(bool::from(maybe_quo.is_some()), exact_and_fits); + if exact_and_fits { assert_eq!(maybe_quo.unwrap(), exp_q); } } @@ -1171,91 +1202,117 @@ mod tests { (U128::MAX, U128::ZERO, U128::from(2u64)), // half-word divisor (two_127, U128::ONE, U128::from(3u64)), // 3 * 2^127, single-word divisor ] { - check::<{ U128::LIMBS }, { U256::LIMBS }>(lo, hi, den); + check_wide_division::<{ U128::LIMBS }, { U256::LIMBS }>(lo, hi, den.to_nz().unwrap()); } // Single-word operands (LIMBS == 1). for (lo, hi, den) in [ (U64::MAX, U64::MAX, U64::MAX), (U64::from(5u64), U64::from(9u64), U64::from(4u64)), ] { - check::<{ U64::LIMBS }, { U128::LIMBS }>(lo, hi, den); + check_wide_division::<{ U64::LIMBS }, { U128::LIMBS }>(lo, hi, den.to_nz().unwrap()); } // A couple of hand-computed values for extra confidence. - // (2^256 - 1) / (2^128 - 1) = 2^128 + 1, wrapped mod 2^128 = 1, remainder 0. - let (q, r) = U128::div_rem_wide((U128::MAX, U128::MAX), &U128::MAX.to_nz().unwrap()); + // (2^256 - 1) / (2^128 - 1) = 2^128 + 1, wrapped mod 2^128 = 1, remainder 0; the true + // quotient exceeds 2^128, so it does not fit. + let (q, r, fits) = + U128::wrapping_div_rem_wide((U128::MAX, U128::MAX), &U128::MAX.to_nz().unwrap()); assert_eq!(q, U128::ONE); assert_eq!(r, U128::ZERO); + assert!(!bool::from(fits)); - // 3 * 2^127 = 2^128 + 2^127 => hi = 1, lo = 2^127; the exact quotient is 2^127. + // 3 * 2^127 = 2^128 + 2^127 => hi = 1, lo = 2^127; the exact quotient is 2^127, which fits. let three = U128::from(3u64).to_nz().unwrap(); - let (q, r) = U128::div_rem_wide((two_127, U128::ONE), &three); + let (q, r, fits) = U128::wrapping_div_rem_wide((two_127, U128::ONE), &three); assert_eq!(q, two_127); assert_eq!(r, U128::ZERO); + assert!(bool::from(fits)); assert_eq!( U128::div_wide_exact((two_127, U128::ONE), &three).unwrap(), two_127 ); - } - /// Force a randomly generated value to be non-zero so it is a valid divisor. - #[cfg(feature = "rand_core")] - fn nz(v: Uint) -> Uint { - if bool::from(v.is_zero()) { - Uint::ONE - } else { - v - } + // A zero remainder alone is not enough for `div_wide_exact`: 2^128 / 1 is exact but the + // quotient overflows `U128`, so the result is absent. + let one = U128::ONE.to_nz().unwrap(); + assert!(bool::from( + U128::div_wide_exact((U128::ZERO, U128::ONE), &one).is_none() + )); } #[cfg(feature = "rand_core")] #[test] fn div_rem_wide_vs_concat() { + /// Random divisor with `K` significant limbs, widened to `Uint` (still non-zero). + fn narrow_nz(rng: &mut ChaCha8Rng) -> NonZero> { + NonZero::>::random_from_rng(rng) + .as_ref() + .resize::() + .to_nz() + .unwrap() + } + let mut rng = ChaCha8Rng::from_seed([9u8; 32]); for _ in 0..300 { // Single-word path (U64 operands). let lo64 = U64::random_from_rng(&mut rng); let hi64 = U64::random_from_rng(&mut rng); - check::<{ U64::LIMBS }, { U128::LIMBS }>( + check_wide_division::<{ U64::LIMBS }, { U128::LIMBS }>( lo64, hi64, - nz(U64::random_from_rng(&mut rng)), + NonZero::random_from_rng(&mut rng), ); // Multi-limb path (U128): a full-width divisor, then a single-word-value divisor // (ywords == 1, which hits the div2by1 correction). let lo = U128::random_from_rng(&mut rng); let hi = U128::random_from_rng(&mut rng); - check::<{ U128::LIMBS }, { U256::LIMBS }>(lo, hi, nz(U128::random_from_rng(&mut rng))); - let d_small = U64::random_from_rng(&mut rng).resize::<{ U128::LIMBS }>(); - check::<{ U128::LIMBS }, { U256::LIMBS }>(lo, hi, nz(d_small)); + check_wide_division::<{ U128::LIMBS }, { U256::LIMBS }>( + lo, + hi, + NonZero::random_from_rng(&mut rng), + ); + check_wide_division::<{ U128::LIMBS }, { U256::LIMBS }>( + lo, + hi, + narrow_nz::<{ U64::LIMBS }, { U128::LIMBS }>(&mut rng), + ); // Odd limb count (U192): a full-width divisor (>= 3 significant limbs, so the tail // `done` masking / vartime early-break fires) and a narrower divisor. let lo192 = U192::random_from_rng(&mut rng); let hi192 = U192::random_from_rng(&mut rng); - check::<{ U192::LIMBS }, { U384::LIMBS }>( + check_wide_division::<{ U192::LIMBS }, { U384::LIMBS }>( lo192, hi192, - nz(U192::random_from_rng(&mut rng)), + NonZero::random_from_rng(&mut rng), + ); + check_wide_division::<{ U192::LIMBS }, { U384::LIMBS }>( + lo192, + hi192, + narrow_nz::<{ U128::LIMBS }, { U192::LIMBS }>(&mut rng), ); - let d192_narrow = U128::random_from_rng(&mut rng).resize::<{ U192::LIMBS }>(); - check::<{ U192::LIMBS }, { U384::LIMBS }>(lo192, hi192, nz(d192_narrow)); // Wider operands (U256): a full-width divisor, a single-word-value divisor, and one // with exactly LIMBS-1 significant limbs (fires the vartime early-break while the // divisor is still narrower than the dividend's high half). let lo256 = U256::random_from_rng(&mut rng); let hi256 = U256::random_from_rng(&mut rng); - check::<{ U256::LIMBS }, { U512::LIMBS }>( + check_wide_division::<{ U256::LIMBS }, { U512::LIMBS }>( + lo256, + hi256, + NonZero::random_from_rng(&mut rng), + ); + check_wide_division::<{ U256::LIMBS }, { U512::LIMBS }>( + lo256, + hi256, + narrow_nz::<{ U128::LIMBS }, { U256::LIMBS }>(&mut rng), + ); + check_wide_division::<{ U256::LIMBS }, { U512::LIMBS }>( lo256, hi256, - nz(U256::random_from_rng(&mut rng)), + narrow_nz::<{ U192::LIMBS }, { U256::LIMBS }>(&mut rng), ); - let d256_narrow = U128::random_from_rng(&mut rng).resize::<{ U256::LIMBS }>(); - check::<{ U256::LIMBS }, { U512::LIMBS }>(lo256, hi256, nz(d256_narrow)); - let d256_3limb = U192::random_from_rng(&mut rng).resize::<{ U256::LIMBS }>(); - check::<{ U256::LIMBS }, { U512::LIMBS }>(lo256, hi256, nz(d256_3limb)); } } } diff --git a/src/uint/ref_type/div.rs b/src/uint/ref_type/div.rs index fd85c5a64..aada2ab70 100644 --- a/src/uint/ref_type/div.rs +++ b/src/uint/ref_type/div.rs @@ -205,7 +205,8 @@ impl UintRef { } /// Computes `x_lower_upper` / `rhs`, returning the wrapped quotient in `quo` and the - /// remainder in `rhs`. + /// remainder in `rhs`. Returns a [`Choice`] that is truthy when the quotient fit in `quo` + /// without truncation. /// /// The `x_lower_upper` tuple represents a wide (double-width) dividend `x_lo + x_hi * B`, /// where `B = 2^(x_lo.bits_precision())`. The size of `x_lower_upper.1` and of `quo` must each @@ -218,14 +219,18 @@ impl UintRef { /// # Panics /// If the divisor is zero. #[inline(always)] - pub(crate) const fn div_rem_wide( + pub(crate) const fn wrapping_div_rem_wide( x_lower_upper: (&mut Self, &mut Self), rhs: &mut Self, quo: &mut Self, - ) { + ) -> Choice { let (x_lo, x) = x_lower_upper; let y = rhs; + // The retained low half of the quotient is the whole quotient exactly when the true + // quotient is below `B`, which happens iff the high half of the dividend is `< rhs`. + let fits = UintRef::lt(x, y); + // Short circuit for single-word divisor (only reachable when the operands are one limb // wide, since `div3by2` in the main path requires a two-limb divisor). if y.nlimbs() == 1 { @@ -244,7 +249,7 @@ impl UintRef { quo.limbs[0] = x_lo.limbs[0]; y.limbs[0] = hi.shr(sh); - return; + return fits; } // Compute the size of the divisor @@ -263,10 +268,12 @@ impl UintRef { x.limbs[0] = x.limbs[0].bitor(x_lo_carry); // Perform the core division algorithm - Self::div_rem_wide_shifted((x_lo, x), x_hi, y, ywords, quo); + Self::wrapping_div_rem_wide_shifted((x_lo, x), x_hi, y, ywords, quo); // Unshift the remainder from the earlier adjustment y.shr_assign_limb(lshift); + + fits } /// Computes `x_lower_upper` / `rhs`, returning the wrapped quotient in `quo` and the @@ -277,18 +284,24 @@ impl UintRef { /// /// The `x_lower_upper` tuple represents a wide (double-width) dividend. The size of /// `x_lower_upper.1` and of `quo` must each be at least as large as `rhs`. `x_lower_upper` is - /// left in an indeterminate state. See [`UintRef::div_rem_wide`] for the wrapping semantics. + /// left in an indeterminate state. See [`UintRef::wrapping_div_rem_wide`] for the wrapping + /// semantics and the returned [`Choice`]. /// /// # Panics /// If the divisor is zero. #[inline(always)] - pub(crate) const fn div_rem_wide_vartime( + pub(crate) const fn wrapping_div_rem_wide_vartime( x_lower_upper: (&mut Self, &mut Self), rhs: &mut Self, quo: &mut Self, - ) { + ) -> Choice { let (x_lo, x) = x_lower_upper; let xsize = x.nlimbs(); + + // The retained low half of the quotient is the whole quotient exactly when the high half + // of the dividend is `< rhs` (see `wrapping_div_rem_wide`). + let fits = UintRef::lt(x, rhs); + let ysize = bitlen::to_limbs(rhs.bits_vartime()); let y = rhs.leading_mut(ysize); @@ -298,7 +311,7 @@ impl UintRef { // Empty dividend: both quotient and remainder are zero. y.fill(Limb::ZERO); quo.fill(Limb::ZERO); - return; + return fits; } (_, 1) => { // Single-word divisor: long division by one limb, most-significant limb first, @@ -326,7 +339,7 @@ impl UintRef { quo.copy_from(x_lo); y.fill(Limb::ZERO); y.limbs[0] = hi.shr(sh); - return; + return fits; } _ if ysize > xsize => { panic!("divisor too large"); @@ -349,7 +362,7 @@ impl UintRef { let reciprocal = Reciprocal::new(y.limbs[ysize - 1].to_nz().expect_copied("zero divisor")); // Perform the core division algorithm - x_hi = Self::div_rem_wide_large_shifted::( + x_hi = Self::wrapping_div_rem_wide_large_shifted::( (x_lo, x), x_hi, y, @@ -367,6 +380,8 @@ impl UintRef { // Unshift the remainder from the earlier adjustment y.shr_assign_limb_vartime(lshift); + + fits } /// Conditionally shift the limbs one position toward the most-significant end, dropping the @@ -377,11 +392,8 @@ impl UintRef { /// low `self.nlimbs()` limbs of the full-width quotient. #[inline(always)] const fn shift_in_limb(&mut self, limb: Limb, shift: Choice) { - let mut j = self.nlimbs(); - while j > 1 { - j -= 1; - self.limbs[j] = Limb::select(self.limbs[j], self.limbs[j - 1], shift); - } + // Shift the limbs up by one, inserting a zero at the bottom, then overwrite it with `limb`. + self.conditional_shl_assign_by_limbs_vartime(1, shift); if self.nlimbs() > 0 { self.limbs[0] = Limb::select(self.limbs[0], limb, shift); } @@ -397,7 +409,7 @@ impl UintRef { /// be unshifted by the caller). `x` is left in an indeterminate state. #[inline(always)] #[allow(clippy::cast_possible_truncation)] - const fn div_rem_wide_shifted( + const fn wrapping_div_rem_wide_shifted( x: (&mut Self, &mut Self), mut x_hi: Limb, y: &mut Self, @@ -412,8 +424,14 @@ impl UintRef { debug_assert!(reciprocal.shift() == 0); // Perform the core division algorithm - x_hi = - Self::div_rem_wide_large_shifted::((x_lo, x), x_hi, y, ywords, reciprocal, quo); + x_hi = Self::wrapping_div_rem_wide_large_shifted::( + (x_lo, x), + x_hi, + y, + ywords, + reciprocal, + quo, + ); // Calculate quotient and remainder for the case where the divisor is a single word. let limb_div = Choice::from_u32_eq(1, ywords); @@ -451,7 +469,7 @@ impl UintRef { /// high bit of the divisor is set, and `x_hi` holds the top bits of the dividend. #[inline(always)] #[allow(clippy::cast_possible_truncation)] - const fn div_rem_wide_large_shifted( + const fn wrapping_div_rem_wide_large_shifted( x: (&Self, &mut Self), mut x_hi: Limb, y: &Self, @@ -461,7 +479,7 @@ impl UintRef { ) -> Limb { assert!( y.nlimbs() <= x.1.nlimbs(), - "invalid input sizes for div_rem_wide_large_shifted" + "invalid input sizes for wrapping_div_rem_wide_large_shifted" ); let (x_lo, x) = x;