From c67d6ac70508157d14fbf3c257c4385632cdb725 Mon Sep 17 00:00:00 2001 From: Bogdan Opanchuk Date: Fri, 9 Aug 2024 20:17:29 -0700 Subject: [PATCH 1/2] Some simplifications of vartime division --- src/uint/boxed/div.rs | 123 ++++++++++------- src/uint/div.rs | 302 ++++++++++++++++++++++-------------------- src/uint/div_limb.rs | 24 ++-- 3 files changed, 249 insertions(+), 200 deletions(-) diff --git a/src/uint/boxed/div.rs b/src/uint/boxed/div.rs index 9d061fc9..2957b20b 100644 --- a/src/uint/boxed/div.rs +++ b/src/uint/boxed/div.rs @@ -1,9 +1,12 @@ //! [`BoxedUint`] division operations. use crate::{ - uint::{boxed, div_limb::div3by2}, + uint::{ + boxed, + div_limb::{div2by1, div3by2}, + }, BoxedUint, CheckedDiv, ConstChoice, ConstantTimeSelect, DivRemLimb, Limb, NonZero, Reciprocal, - RemLimb, Word, Wrapping, + RemLimb, Wrapping, }; use core::ops::{Div, DivAssign, Rem, RemAssign}; use subtle::CtOption; @@ -166,7 +169,7 @@ impl BoxedUint { while xi > 0 { // Divide high dividend words by the high divisor word to estimate the quotient word - let (mut quo, _) = div3by2(x_hi.0, x_lo.0, x[xi - 1].0, &reciprocal, y[size - 2].0); + let mut quo = div3by2(x_hi.0, x_lo.0, x[xi - 1].0, &reciprocal, y[size - 2].0); // This loop is a no-op once xi is smaller than the number of words in the divisor let done = ConstChoice::from_u32_lt(xi as u32, dwords - 1); @@ -206,14 +209,20 @@ impl BoxedUint { } let limb_div = ConstChoice::from_u32_eq(1, dwords); + // Calculate quotient and remainder for the case where the divisor is a single word - let (quo2, rem2) = div3by2(x_hi.0, x_lo.0, 0, &reciprocal, 0); + // 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_hi_adjusted.0, x_lo.0, &reciprocal); // Adjust the quotient for single limb division x[0] = Limb::select(x[0], Limb(quo2), limb_div); // Copy out the remainder - y[0] = Limb::select(x[0], Limb(rem2 as Word), limb_div); + y[0] = Limb::select(x[0], Limb(rem2), limb_div); i = 1; while i < size { y[i] = Limb::select(Limb::ZERO, x[i], ConstChoice::from_u32_lt(i as u32, dwords)); @@ -382,6 +391,42 @@ impl RemLimb for BoxedUint { } } +/// Computes `limbs << shift` inplace, where `0 <= shift < Limb::BITS`, returning the carry. +fn shl_limb_vartime(limbs: &mut [Limb], shift: u32) -> Limb { + if shift == 0 { + return Limb::ZERO; + } + + let lshift = shift; + let rshift = Limb::BITS - shift; + let limbs_num = limbs.len(); + + let carry = limbs[limbs_num - 1] >> rshift; + for i in (1..limbs_num).rev() { + limbs[i] = (limbs[i] << lshift) | (limbs[i - 1] >> rshift); + } + limbs[0] <<= lshift; + + carry +} + +/// Computes `limbs >> shift` inplace, where `0 <= shift < Limb::BITS`. +fn shr_limb_vartime(limbs: &mut [Limb], shift: u32) { + if shift == 0 { + return; + } + + let lshift = Limb::BITS - shift; + let rshift = shift; + + let limbs_num = limbs.len(); + + for i in 0..limbs_num - 1 { + limbs[i] = (limbs[i] >> rshift) | (limbs[i + 1] << lshift); + } + limbs[limbs_num - 1] >>= rshift; +} + /// Computes `x` / `y`, returning the quotient in `x` and the remainder in `y`. /// /// This function operates in variable-time. It will panic if the divisor is zero @@ -408,51 +453,44 @@ pub(crate) fn div_rem_vartime_in_place(x: &mut [Limb], y: &mut [Limb]) { } let lshift = y[yc - 1].leading_zeros(); - let rshift = if lshift == 0 { 0 } else { Limb::BITS - lshift }; - let mut x_hi = Limb::ZERO; - let mut carry; - - if lshift != 0 { - // 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 - carry = Limb::ZERO; - for i in 0..yc { - (y[i], carry) = (Limb((y[i].0 << lshift) | carry.0), Limb(y[i].0 >> rshift)); - } - // Shift the dividend to match - carry = Limb::ZERO; - for i in 0..xc { - (x[i], carry) = (Limb((x[i].0 << lshift) | carry.0), Limb(x[i].0 >> rshift)); - } - x_hi = carry; - } + // 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 + shl_limb_vartime(y, lshift); + + // Shift the dividend to match + let mut x_hi = shl_limb_vartime(x, lshift); let reciprocal = Reciprocal::new(y[yc - 1].to_nz().expect("zero divisor")); for xi in (yc - 1..xc).rev() { // Divide high dividend words by the high divisor word to estimate the quotient word - let (mut quo, _) = div3by2(x_hi.0, x[xi].0, x[xi - 1].0, &reciprocal, y[yc - 2].0); + let mut quo = div3by2(x_hi.0, x[xi].0, x[xi - 1].0, &reciprocal, y[yc - 2].0); // Subtract q*divisor from the dividend - carry = Limb::ZERO; - let mut borrow = Limb::ZERO; - let mut tmp; - for i in 0..yc { - (tmp, carry) = Limb::ZERO.mac(y[i], Limb(quo), carry); - (x[xi + i + 1 - yc], borrow) = x[xi + i + 1 - yc].sbb(tmp, borrow); - } - (_, borrow) = x_hi.sbb(carry, borrow); + let borrow = { + let mut carry = Limb::ZERO; + let mut borrow = Limb::ZERO; + let mut tmp; + for i in 0..yc { + (tmp, carry) = Limb::ZERO.mac(y[i], Limb(quo), carry); + (x[xi + i + 1 - yc], borrow) = x[xi + i + 1 - yc].sbb(tmp, borrow); + } + (_, borrow) = x_hi.sbb(carry, borrow); + borrow + }; // If the subtraction borrowed, then decrement q and add back the divisor // The probability of this being needed is very low, about 2/(Limb::MAX+1) - let ct_borrow = ConstChoice::from_word_mask(borrow.0); - carry = Limb::ZERO; - for i in 0..yc { - (x[xi + i + 1 - yc], carry) = - x[xi + i + 1 - yc].adc(Limb::select(Limb::ZERO, y[i], ct_borrow), carry); - } - quo = ct_borrow.select_word(quo, quo.saturating_sub(1)); + quo = { + let ct_borrow = ConstChoice::from_word_mask(borrow.0); + let mut carry = Limb::ZERO; + for i in 0..yc { + (x[xi + i + 1 - yc], carry) = + x[xi + i + 1 - yc].adc(Limb::select(Limb::ZERO, y[i], ct_borrow), carry); + } + ct_borrow.select_word(quo, quo.saturating_sub(1)) + }; // Store the quotient within dividend and set x_hi to the current highest word x_hi = x[xi]; @@ -464,12 +502,7 @@ pub(crate) fn div_rem_vartime_in_place(x: &mut [Limb], y: &mut [Limb]) { y[yc - 1] = x_hi; // Unshift the remainder from the earlier adjustment - if lshift != 0 { - carry = Limb::ZERO; - for i in (0..yc).rev() { - (y[i], carry) = (Limb((y[i].0 >> lshift) | carry.0), Limb(y[i].0 << rshift)); - } - } + shr_limb_vartime(y, lshift); // Shift the quotient to the low limbs within dividend // let x_size = xc - yc + 1; diff --git a/src/uint/div.rs b/src/uint/div.rs index 14a06199..09235390 100644 --- a/src/uint/div.rs +++ b/src/uint/div.rs @@ -1,10 +1,10 @@ //! [`Uint`] division operations. use super::div_limb::{ - div3by2, div_rem_limb_with_reciprocal, rem_limb_with_reciprocal, rem_limb_with_reciprocal_wide, - Reciprocal, + div2by1, div3by2, div_rem_limb_with_reciprocal, rem_limb_with_reciprocal, + rem_limb_with_reciprocal_wide, Reciprocal, }; -use crate::{CheckedDiv, ConstChoice, DivRemLimb, Limb, NonZero, RemLimb, Uint, Word, Wrapping}; +use crate::{CheckedDiv, ConstChoice, DivRemLimb, Limb, NonZero, RemLimb, Uint, Wrapping}; use core::ops::{Div, DivAssign, Rem, RemAssign}; use subtle::CtOption; @@ -69,7 +69,7 @@ impl Uint { while xi > 0 { // Divide high dividend words by the high divisor word to estimate the quotient word - let (mut quo, _) = div3by2(x_hi.0, x_lo.0, x[xi - 1].0, &reciprocal, y[LIMBS - 2].0); + let mut quo = div3by2(x_hi.0, x_lo.0, x[xi - 1].0, &reciprocal, y[LIMBS - 2].0); // This loop is a no-op once xi is smaller than the number of words in the divisor let done = ConstChoice::from_u32_lt(xi as u32, dwords - 1); @@ -109,14 +109,20 @@ impl Uint { } let limb_div = ConstChoice::from_u32_eq(1, dwords); + // Calculate quotient and remainder for the case where the divisor is a single word - let (quo2, rem2) = div3by2(x_hi.0, x_lo.0, 0, &reciprocal, 0); + // 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_hi_adjusted.0, x_lo.0, &reciprocal); // Adjust the quotient for single limb division x[0] = Limb::select(x[0], Limb(quo2), limb_div); // Copy out the remainder - y[0] = Limb::select(x[0], Limb(rem2 as Word), limb_div); + y[0] = Limb::select(x[0], Limb(rem2), limb_div); i = 1; while i < LIMBS { y[i] = Limb::select(Limb::ZERO, x[i], ConstChoice::from_u32_lt(i as u32, dwords)); @@ -130,6 +136,54 @@ impl Uint { ) } + /// Computes `self << shift` where `0 <= shift < Limb::BITS`, + /// returning the result and the carry. + /// + /// Note: assumes that `self` only has `limb_num` lowest non-zero limbs. + const fn shl_limb_vartime(&self, shift: u32, limbs_num: usize) -> (Self, Limb) { + if shift == 0 { + return (*self, Limb::ZERO); + } + + let mut limbs = [Limb::ZERO; LIMBS]; + + let lshift = shift; + let rshift = Limb::BITS - shift; + + let carry = self.limbs[limbs_num - 1].0 >> rshift; + let mut i = limbs_num - 1; + while i > 0 { + limbs[i] = Limb((self.limbs[i].0 << lshift) | (self.limbs[i - 1].0 >> rshift)); + i -= 1; + } + limbs[0] = Limb(self.limbs[0].0 << lshift); + + (Uint::::new(limbs), Limb(carry)) + } + + /// Computes `self >> shift` where `0 <= shift < Limb::BITS`. + /// + /// Note: assumes that `self` only has `limb_num` lowest non-zero limbs. + const fn shr_limb_vartime(&self, shift: u32, limbs_num: usize) -> Self { + if shift == 0 { + return *self; + } + + let mut limbs = [Limb::ZERO; LIMBS]; + + let lshift = Limb::BITS - shift; + let rshift = shift; + + let mut i = 0; + while i < limbs_num - 1 { + limbs[i] = Limb((self.limbs[i].0 >> rshift) | (self.limbs[i + 1].0 << lshift)); + i += 1; + } + limbs[limbs_num - 1] = Limb(self.limbs[limbs_num - 1].0 >> rshift); + + Uint::::new(limbs) + } + /// Computes `self` / `rhs`, returns the quotient (q) and the remainder (r) /// /// This is variable only with respect to `rhs`. @@ -147,81 +201,67 @@ impl Uint { let yc = ((dbits + Limb::BITS - 1) / Limb::BITS) as usize; // Short circuit for small or extra large divisors - match yc { - 1 => { - // If the divisor is a single limb, use limb division - let (q, r) = div_rem_limb_with_reciprocal( - self, - &Reciprocal::new(rhs.0.limbs[0].to_nz().expect("zero divisor")), - ); - return (q, Uint::from_word(r.0)); - } - yc if yc > LIMBS => { - // Divisor is greater than dividend. Return zero and the dividend as the - // quotient and remainder - return (Uint::ZERO, self.resize()); - } - _ => {} - }; - - let lshift = (Limb::BITS - (dbits % Limb::BITS)) % Limb::BITS; - let rshift = if lshift == 0 { 0 } else { Limb::BITS - lshift }; - let mut x = self.to_limbs(); - let mut x_hi = Limb::ZERO; - let mut xi = LIMBS - 1; - let mut y = rhs.0.to_limbs(); - let mut i; - let mut carry; + if yc == 1 { + // If the divisor is a single limb, use limb division + let (q, r) = div_rem_limb_with_reciprocal( + self, + &Reciprocal::new(rhs.0.limbs[0].to_nz().expect("zero divisor")), + ); + return (q, Uint::from_word(r.0)); + } + if yc > LIMBS { + // Divisor is greater than dividend. Return zero and the dividend as the + // quotient and remainder + return (Uint::ZERO, self.resize()); + } - if lshift != 0 { - // 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 - i = 0; - carry = Limb::ZERO; - while i < yc { - (y[i], carry) = (Limb((y[i].0 << lshift) | carry.0), Limb(y[i].0 >> rshift)); - i += 1; - } + // The shift needed to set the MSB of the highest nonzero limb of the divisor. + // 2^shift == d in the algorithm above. + let shift = (Limb::BITS - (dbits % Limb::BITS)) % Limb::BITS; - // Shift the dividend to match - i = 0; - carry = Limb::ZERO; - while i < LIMBS { - (x[i], carry) = (Limb((x[i].0 << lshift) | carry.0), Limb(x[i].0 >> rshift)); - i += 1; - } - x_hi = carry; - } + let (x, mut x_hi) = self.shl_limb_vartime(shift, LIMBS); + let mut x = x.to_limbs(); + let (y, _) = rhs.0.shl_limb_vartime(shift, yc); + let mut y = y.to_limbs(); let reciprocal = Reciprocal::new(y[yc - 1].to_nz().expect("zero divisor")); + let mut i; + + let mut xi = LIMBS - 1; + loop { // Divide high dividend words by the high divisor word to estimate the quotient word - let (mut quo, _) = div3by2(x_hi.0, x[xi].0, x[xi - 1].0, &reciprocal, y[yc - 2].0); + let mut quo = div3by2(x_hi.0, x[xi].0, x[xi - 1].0, &reciprocal, y[yc - 2].0); // Subtract q*divisor from the dividend - carry = Limb::ZERO; - let mut borrow = Limb::ZERO; - let mut tmp; - i = 0; - while i < yc { - (tmp, carry) = Limb::ZERO.mac(y[i], Limb(quo), carry); - (x[xi + i + 1 - yc], borrow) = x[xi + i + 1 - yc].sbb(tmp, borrow); - i += 1; - } - (_, borrow) = x_hi.sbb(carry, borrow); + let borrow = { + let mut carry = Limb::ZERO; + let mut borrow = Limb::ZERO; + let mut tmp; + i = 0; + while i < yc { + (tmp, carry) = Limb::ZERO.mac(y[i], Limb(quo), carry); + (x[xi + i + 1 - yc], borrow) = x[xi + i + 1 - yc].sbb(tmp, borrow); + i += 1; + } + (_, borrow) = x_hi.sbb(carry, borrow); + borrow + }; // If the subtraction borrowed, then decrement q and add back the divisor // The probability of this being needed is very low, about 2/(Limb::MAX+1) - let ct_borrow = ConstChoice::from_word_mask(borrow.0); - carry = Limb::ZERO; - i = 0; - while i < yc { - (x[xi + i + 1 - yc], carry) = - x[xi + i + 1 - yc].adc(Limb::select(Limb::ZERO, y[i], ct_borrow), carry); - i += 1; - } - quo = ct_borrow.select_word(quo, quo.saturating_sub(1)); + quo = { + let ct_borrow = ConstChoice::from_word_mask(borrow.0); + let mut carry = Limb::ZERO; + i = 0; + while i < yc { + (x[xi + i + 1 - yc], carry) = + x[xi + i + 1 - yc].adc(Limb::select(Limb::ZERO, y[i], ct_borrow), carry); + i += 1; + } + ct_borrow.select_word(quo, quo.saturating_sub(1)) + }; // Store the quotient within dividend and set x_hi to the current highest word x_hi = x[xi]; @@ -242,14 +282,7 @@ impl Uint { y[yc - 1] = x_hi; // Unshift the remainder from the earlier adjustment - if lshift != 0 { - i = yc; - carry = Limb::ZERO; - while i > 0 { - i -= 1; - (y[i], carry) = (Limb((y[i].0 >> lshift) | carry.0), Limb(y[i].0 << rshift)); - } - } + let y = Uint::new(y).shr_limb_vartime(shift, yc); // Shift the quotient to the low limbs within dividend i = 0; @@ -262,7 +295,7 @@ impl Uint { i += 1; } - (Uint::new(x), Uint::new(y)) + (Uint::new(x), y) } /// Computes `self` % `rhs`, returns the remainder. @@ -297,63 +330,63 @@ impl Uint { return Uint::from_word(r.0); } - let lshift = (Limb::BITS - (dbits % Limb::BITS)) % Limb::BITS; - let rshift = if lshift == 0 { 0 } else { Limb::BITS - lshift }; - let mut x = lower_upper.1.to_limbs(); // high limbs - let mut x_hi = Limb::ZERO; - let mut xi = LIMBS - 1; - let mut y = rhs.0.to_limbs(); - let mut extra_limbs = LIMBS; - let mut i; - let mut carry; + // The shift needed to set the MSB of the highest nonzero limb of the divisor. + // 2^shift == d in the algorithm above. + let shift = (Limb::BITS - (dbits % Limb::BITS)) % Limb::BITS; - if lshift != 0 { - // Shift divisor such that it has no leading zeros - // This ensures that the high word >= b/2, and means that div2by1 requires no extra shifts - i = 0; - carry = Limb::ZERO; - while i < yc { - (y[i], carry) = (Limb((y[i].0 << lshift) | carry.0), Limb(y[i].0 >> rshift)); - i += 1; - } + let (y, _) = rhs.0.shl_limb_vartime(shift, yc); + let y = y.to_limbs(); - // Shift the dividend to match - i = 0; - carry = Limb(lower_upper.0.limbs[LIMBS - 1].0 >> rshift); - while i < LIMBS { - (x[i], carry) = (Limb((x[i].0 << lshift) | carry.0), Limb(x[i].0 >> rshift)); - i += 1; - } - x_hi = carry; + let (x_lo, x_lo_carry) = lower_upper.0.shl_limb_vartime(shift, LIMBS); + let (x, mut x_hi) = lower_upper.1.shl_limb_vartime(shift, LIMBS); + let mut x = x.to_limbs(); + if shift > 0 { + x[0] = Limb(x[0].0 | x_lo_carry.0); } let reciprocal = Reciprocal::new(y[yc - 1].to_nz().expect("zero divisor")); + let mut xi = LIMBS - 1; + let mut extra_limbs = LIMBS; + let mut i; + + // Note that in the algorithm we only ever need to access the highest `yc` limbs + // of the dividend, and since `yc < LIMBS`, we only need to access + // the high half of the dividend. + // + // So we proceed similarly to `div_rem_vartime()` applied to the high half of the dividend, + // fetching the limbs from the lower part as we go. + loop { // Divide high dividend words by the high divisor word to estimate the quotient word - let (quo, _) = div3by2(x_hi.0, x[xi].0, x[xi - 1].0, &reciprocal, y[yc - 2].0); + let quo = div3by2(x_hi.0, x[xi].0, x[xi - 1].0, &reciprocal, y[yc - 2].0); // Subtract q*divisor from the dividend - carry = Limb::ZERO; - let mut borrow = Limb::ZERO; - let mut tmp; - i = 0; - while i < yc { - (tmp, carry) = Limb::ZERO.mac(y[i], Limb(quo), carry); - (x[xi + i + 1 - yc], borrow) = x[xi + i + 1 - yc].sbb(tmp, borrow); - i += 1; - } - (_, borrow) = x_hi.sbb(carry, borrow); + let borrow = { + let mut carry = Limb::ZERO; + let mut borrow = Limb::ZERO; + let mut tmp; + i = 0; + while i < yc { + (tmp, carry) = Limb::ZERO.mac(y[i], Limb(quo), carry); + (x[xi + i + 1 - yc], borrow) = x[xi + i + 1 - yc].sbb(tmp, borrow); + i += 1; + } + (_, borrow) = x_hi.sbb(carry, borrow); + borrow + }; // If the subtraction borrowed, then add back the divisor // The probability of this being needed is very low, about 2/(Limb::MAX+1) - let ct_borrow = ConstChoice::from_word_mask(borrow.0); - carry = Limb::ZERO; - i = 0; - while i < yc { - (x[xi + i + 1 - yc], carry) = - x[xi + i + 1 - yc].adc(Limb::select(Limb::ZERO, y[i], ct_borrow), carry); - i += 1; + { + let ct_borrow = ConstChoice::from_word_mask(borrow.0); + let mut carry = Limb::ZERO; + i = 0; + while i < yc { + (x[xi + i + 1 - yc], carry) = + x[xi + i + 1 - yc].adc(Limb::select(Limb::ZERO, y[i], ct_borrow), carry); + i += 1; + } } // Set x_hi to the current highest word @@ -367,13 +400,7 @@ impl Uint { x[i] = x[i - 1]; i -= 1; } - x[0] = lower_upper.0.limbs[extra_limbs]; - if lshift != 0 { - x[0].0 <<= lshift; - if extra_limbs > 0 { - x[0].0 |= lower_upper.0.limbs[extra_limbs - 1].0 >> rshift; - } - } + x[0] = x_lo.limbs[extra_limbs]; } else { if xi == yc - 1 { break; @@ -383,22 +410,7 @@ impl Uint { } // Unshift the remainder from the earlier adjustment - if lshift != 0 { - i = yc; - carry = Limb::ZERO; - while i > 0 { - i -= 1; - (x[i], carry) = (Limb((x[i].0 >> lshift) | carry.0), Limb(x[i].0 << rshift)); - } - } - // Clear upper limbs - i = LIMBS - 1; - while i >= yc { - x[i] = Limb::ZERO; - i -= 1; - } - - Uint::new(x) + Uint::new(x).shr_limb_vartime(shift, yc) } /// Computes `self` % 2^k. Faster than reduce since its a power of 2. diff --git a/src/uint/div_limb.rs b/src/uint/div_limb.rs index 595518b8..d9142840 100644 --- a/src/uint/div_limb.rs +++ b/src/uint/div_limb.rs @@ -145,23 +145,27 @@ pub(crate) const fn div2by1(u1: Word, u0: Word, reciprocal: &Reciprocal) -> (Wor (q1, r) } -/// Calculate the quotient and the remainder of the division of a 3-word -/// dividend `u`, with a precalculated leading-word reciprocal `v` and a second -/// divisor word `v0`. The dividend and the lower divisor word must be left-shifted -/// according to the `shift` value of the reciprocal. +/// Given two long integers `u = (..., u0, u1, u2)` and `v = (..., v0, v1)` +/// (where `u2` and `v1` are the most significant limbs), where `floor(u / v) <= Limb::MAX`, +/// calculates `q` such that `q - 1 <= floor(u / v) <= q`. +/// In place of `v1` takes its reciprocal, and assumes that `v` was already pre-shifted +/// so that v1 has its most significant bit set (that is, the reciprocal's `shift` is 0). #[inline(always)] pub(crate) const fn div3by2( u2: Word, u1: Word, u0: Word, - reciprocal: &Reciprocal, + v1_reciprocal: &Reciprocal, v0: Word, -) -> (Word, WideWord) { +) -> Word { + debug_assert!(v1_reciprocal.shift == 0); + debug_assert!(u2 <= v1_reciprocal.divisor_normalized); + // This method corresponds to Algorithm Q: // https://janmr.com/blog/2014/04/basic-multiple-precision-long-division/ - let q_maxed = ConstChoice::from_word_eq(u2, reciprocal.divisor_normalized); - let (mut quo, rem) = div2by1(q_maxed.select_word(u2, 0), u1, reciprocal); + let q_maxed = ConstChoice::from_word_eq(u2, v1_reciprocal.divisor_normalized); + let (mut quo, rem) = div2by1(q_maxed.select_word(u2, 0), u1, v1_reciprocal); // When the leading dividend word equals the leading divisor word, cap the quotient // at Word::MAX and set the remainder to the sum of the top dividend words. quo = q_maxed.select_word(quo, Word::MAX); @@ -175,11 +179,11 @@ pub(crate) const fn div3by2( let done = ConstChoice::from_word_nonzero((rem >> Word::BITS) as Word) .or(ConstChoice::from_wide_word_le(qy, rx)); quo = done.select_word(quo.saturating_sub(1), quo); - rem = done.select_wide_word(rem + (reciprocal.divisor_normalized as WideWord), rem); + rem = done.select_wide_word(rem + (v1_reciprocal.divisor_normalized as WideWord), rem); i += 1; } - (quo, rem) + quo } /// A pre-calculated reciprocal for division by a single limb. From bec0f0c3719340178639d21481cb4688366a6dfc Mon Sep 17 00:00:00 2001 From: Bogdan Opanchuk Date: Fri, 16 Aug 2024 10:54:28 -0700 Subject: [PATCH 2/2] Use wrapping subs instead of saturating subs --- src/uint/boxed/div.rs | 2 +- src/uint/div.rs | 2 +- src/uint/div_limb.rs | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/uint/boxed/div.rs b/src/uint/boxed/div.rs index 2957b20b..6663991f 100644 --- a/src/uint/boxed/div.rs +++ b/src/uint/boxed/div.rs @@ -489,7 +489,7 @@ pub(crate) fn div_rem_vartime_in_place(x: &mut [Limb], y: &mut [Limb]) { (x[xi + i + 1 - yc], carry) = x[xi + i + 1 - yc].adc(Limb::select(Limb::ZERO, y[i], ct_borrow), carry); } - ct_borrow.select_word(quo, quo.saturating_sub(1)) + ct_borrow.select_word(quo, quo.wrapping_sub(1)) }; // Store the quotient within dividend and set x_hi to the current highest word diff --git a/src/uint/div.rs b/src/uint/div.rs index 09235390..782c7bc9 100644 --- a/src/uint/div.rs +++ b/src/uint/div.rs @@ -260,7 +260,7 @@ impl Uint { x[xi + i + 1 - yc].adc(Limb::select(Limb::ZERO, y[i], ct_borrow), carry); i += 1; } - ct_borrow.select_word(quo, quo.saturating_sub(1)) + ct_borrow.select_word(quo, quo.wrapping_sub(1)) }; // Store the quotient within dividend and set x_hi to the current highest word diff --git a/src/uint/div_limb.rs b/src/uint/div_limb.rs index d9142840..775bac64 100644 --- a/src/uint/div_limb.rs +++ b/src/uint/div_limb.rs @@ -178,7 +178,7 @@ pub(crate) const fn div3by2( // If r < b and q*y[-2] > r*x[-1], then set q = q - 1 and r = r + v1 let done = ConstChoice::from_word_nonzero((rem >> Word::BITS) as Word) .or(ConstChoice::from_wide_word_le(qy, rx)); - quo = done.select_word(quo.saturating_sub(1), quo); + quo = done.select_word(quo.wrapping_sub(1), quo); rem = done.select_wide_word(rem + (v1_reciprocal.divisor_normalized as WideWord), rem); i += 1; }