Implemented div_limb and div_rem_limb_with_reciprocal

Author: xuganyu96

src/uint/boxed.rs

@@ -8,6 +8,7 @@ mod bits;
 mod cmp;
 mod ct;
 mod div;
+mod div_limb;
 pub(crate) mod encoding;
 mod inv_mod;
 mod mul;
@@ -51,6 +52,20 @@ pub struct BoxedUint {
 }
 
 impl BoxedUint {
+    /// Convenience function to instantiate BoxedUint from Limb(s)
+    pub fn from_limbs(limbs: Vec<Limb>) -> Self {
+        return Self {
+            limbs: limbs.into(),
+        };
+    }
+
+    /// Instantiate BoxedUint by cloning from a slice of Limb(s)
+    pub fn copy_from_limb_slice(limbs: &[Limb]) -> Self {
+        return Self {
+            limbs: limbs.to_vec().into(),
+        };
+    }
+
     fn limbs_for_precision(at_least_bits_precision: u32) -> usize {
         ((at_least_bits_precision + Limb::BITS - 1) / Limb::BITS) as usize
     }

src/uint/boxed/div.rs

@@ -1,9 +1,11 @@
 //! [`BoxedUint`] division operations.
 
-use crate::{BoxedUint, CheckedDiv, Limb, NonZero, Wrapping, Reciprocal};
+use crate::{BoxedUint, CheckedDiv, Limb, NonZero, Reciprocal, Wrapping};
 use core::ops::{Div, DivAssign, Rem, RemAssign};
 use subtle::{Choice, ConstantTimeEq, ConstantTimeLess, CtOption};
 
+use super::div_limb::div_rem_limb_with_reciprocal;
+
 impl BoxedUint {
     /// Computes self / rhs, returns the quotient, remainder.
     pub fn div_rem(&self, rhs: &NonZero<Self>) -> (Self, Self) {
@@ -24,11 +26,17 @@ impl BoxedUint {
         self.div_rem_vartime_unchecked(rhs.as_ref())
     }
 
-    /// Computes `self` / `rhs` using a pre-made reciprocal
-    /// returns the quotient (q) and remainder (r)
+    /// Computes `self` / `rhs` using a pre-made reciprocal,
+    /// returns the quotient (q) and remainder (r).
+    #[inline(always)]
     pub fn div_rem_limb_with_reciprocal(&self, reciprocal: &Reciprocal) -> (Self, Limb) {
-        todo!("Need to implement shl_limb");
-        // No need to implement div2by1
+        div_rem_limb_with_reciprocal(self, reciprocal)
+    }
+
+    /// Computes `self` / `rhs`, returns the quotient (q) and remainder (r).
+    #[inline(always)]
+    pub fn div_rem_limb(&self, rhs: NonZero<Limb>) -> (Self, Limb) {
+        div_rem_limb_with_reciprocal(self, &Reciprocal::new(rhs))
     }
 
     /// Computes self % rhs, returns the remainder.

src/uint/boxed/div_limb.rs

@@ -0,0 +1,90 @@
+//! Implementation of constant-time division via reciprocal precomputation, as described in
+//! "Improved Division by Invariant Integers" by Niels Möller and Torbjorn Granlund
+//! (DOI: 10.1109/TC.2010.143, <https://gmplib.org/~tege/division-paper.pdf>).
+use crate::{BoxedUint, ConstChoice, Limb, Reciprocal, WideWord, Word};
+
+/// Multiplies `x` and `y`, returning the most significant
+/// and the least significant words as `(hi, lo)`.
+#[inline(always)]
+const fn mulhilo(x: Word, y: Word) -> (Word, Word) {
+    let res = (x as WideWord) * (y as WideWord);
+    ((res >> Word::BITS) as Word, res as Word)
+}
+
+/// Adds wide numbers represented by pairs of (most significant word, least significant word)
+/// and returns the result in the same format `(hi, lo)`.
+#[inline(always)]
+const fn addhilo(x_hi: Word, x_lo: Word, y_hi: Word, y_lo: Word) -> (Word, Word) {
+    let res = (((x_hi as WideWord) << Word::BITS) | (x_lo as WideWord))
+        + (((y_hi as WideWord) << Word::BITS) | (y_lo as WideWord));
+    ((res >> Word::BITS) as Word, res as Word)
+}
+
+/// Calculate the quotient and the remainder of the division of a wide word
+/// (supplied as high and low words) by `d`, with a precalculated reciprocal `v`.
+#[inline(always)]
+const fn div2by1(u1: Word, u0: Word, reciprocal: &Reciprocal) -> (Word, Word) {
+    let d = reciprocal.divisor_normalized;
+
+    debug_assert!(d >= (1 << (Word::BITS - 1)));
+    debug_assert!(u1 < d);
+
+    let (q1, q0) = mulhilo(reciprocal.reciprocal, u1);
+    let (q1, q0) = addhilo(q1, q0, u1, u0);
+    let q1 = q1.wrapping_add(1);
+    let r = u0.wrapping_sub(q1.wrapping_mul(d));
+
+    let r_gt_q0 = ConstChoice::from_word_lt(q0, r);
+    let q1 = r_gt_q0.select_word(q1, q1.wrapping_sub(1));
+    let r = r_gt_q0.select_word(r, r.wrapping_add(d));
+
+    // If this was a normal `if`, we wouldn't need wrapping ops, because there would be no overflow.
+    // But since we calculate both results either way, we have to wrap.
+    // Added an assert to still check the lack of overflow in debug mode.
+    debug_assert!(r < d || q1 < Word::MAX);
+    let r_ge_d = ConstChoice::from_word_le(d, r);
+    let q1 = r_ge_d.select_word(q1, q1.wrapping_add(1));
+    let r = r_ge_d.select_word(r, r.wrapping_sub(d));
+
+    (q1, r)
+}
+
+/// Divides `u` by the divisor encoded in the `reciprocal`, and returns
+/// the quotient and the remainder.
+#[inline(always)]
+#[allow(dead_code, unused_variables)]
+pub(crate) fn div_rem_limb_with_reciprocal(
+    u: &BoxedUint,
+    reciprocal: &Reciprocal,
+) -> (BoxedUint, Limb) {
+    let nlimbs = u.as_limbs().len();
+    let (u_shifted, u_hi) = u.shl_limb(reciprocal.shift);
+    let mut r = u_hi.0;
+    let mut q = vec![Limb::ZERO; nlimbs];
+
+    let mut j = nlimbs;
+    while j > 0 {
+        j -= 1;
+        let (qj, rj) = div2by1(r, u_shifted.as_limbs()[j].0, reciprocal);
+        q[j] = Limb(qj);
+        r = rj;
+    }
+    (BoxedUint { limbs: q.into() }, Limb(r >> reciprocal.shift))
+}
+
+#[cfg(test)]
+mod tests {
+    use super::{div2by1, Reciprocal};
+    use crate::{Limb, NonZero, Word};
+    #[test]
+    fn div2by1_overflow() {
+        // A regression test for a situation when in div2by1() an operation (`q1 + 1`)
+        // that is protected from overflowing by a condition in the original paper (`r >= d`)
+        // still overflows because we're calculating the results for both branches.
+        let r = Reciprocal::new(NonZero::new(Limb(Word::MAX - 1)).unwrap());
+        assert_eq!(
+            div2by1(Word::MAX - 2, Word::MAX - 63, &r),
+            (Word::MAX, Word::MAX - 65)
+        );
+    }
+}

src/uint/boxed/shl.rs

@@ -1,6 +1,6 @@
 //! [`BoxedUint`] bitwise left shift operations.
 
-use crate::{BoxedUint, ConstChoice, Limb, Zero};
+use crate::{BoxedUint, ConstChoice, Limb, Word, Zero};
 use core::ops::{Shl, ShlAssign};
 use subtle::{Choice, ConstantTimeLess};
 
@@ -71,6 +71,36 @@ impl BoxedUint {
         Some(())
     }
 
+    /// Computes `self << shift` where `0 <= shift < Limb::BITS`,
+    /// returning the result and the carry.
+    #[inline(always)]
+    pub(crate) fn shl_limb(&self, shift: u32) -> (Self, Limb) {
+        let nlimbs = self.limbs.len();
+        let mut limbs = vec![Limb::ZERO; nlimbs];
+
+        let nz = ConstChoice::from_u32_nonzero(shift);
+        let lshift = shift;
+        let rshift = nz.if_true_u32(Limb::BITS - shift);
+        let carry = nz.if_true_word(self.limbs[nlimbs - 1].0.wrapping_shr(Word::BITS - shift));
+
+        limbs[0] = Limb(self.limbs[0].0 << lshift);
+        let mut i = 1;
+        while i < nlimbs {
+            let mut limb = self.limbs[i].0 << lshift;
+            let hi = self.limbs[i - 1].0 >> rshift;
+            limb |= nz.if_true_word(hi);
+            limbs[i] = Limb(limb);
+            i += 1
+        }
+
+        (
+            BoxedUint {
+                limbs: limbs.into(),
+            },
+            Limb(carry),
+        )
+    }
+
     /// Computes `self << shift`.
     /// Returns `None` if `shift >= self.bits_precision()`.
     ///

src/uint/div_limb.rs

@@ -160,9 +160,9 @@ const fn div2by1(u1: Word, u0: Word, reciprocal: &Reciprocal) -> (Word, Word) {
 /// A pre-calculated reciprocal for division by a single limb.
 #[derive(Copy, Clone, Debug, PartialEq, Eq)]
 pub struct Reciprocal {
-    divisor_normalized: Word,
-    shift: u32,
-    reciprocal: Word,
+    pub(crate) divisor_normalized: Word,
+    pub(crate) shift: u32,
+    pub(crate) reciprocal: Word,
 }
 
 impl Reciprocal {