Implement modular arithmetic; MSRV 1.61 (#130)

Adds the following functionality:

- Shifting left and right by 1 and returning the overflow as a `Choice`
- `Neg` for `Wrapping<Uint<LIMBS>>`
- Modular multiplicative inverse
- Computing parameters for the Montgomery form (MF)
- Going to and from MF
- Modular addition (in MF)
- Modular multiplication (in MF)
- Modular exponentiation (in MF)

Author: jellevos

.github/workflows/crypto-bigint.yml

@@ -19,7 +19,7 @@ jobs:
     strategy:
       matrix:
         rust:
-          - 1.57.0 # MSRV
+          - 1.61.0 # MSRV
           - stable
         target:
           - thumbv7em-none-eabi
@@ -49,15 +49,15 @@ jobs:
         include:
           # 32-bit Linux
           - target: i686-unknown-linux-gnu
-            rust: 1.57.0 # MSRV
+            rust: 1.61.0 # MSRV
             deps: sudo apt update && sudo apt install gcc-multilib
           - target: i686-unknown-linux-gnu
             rust: stable
             deps: sudo apt update && sudo apt install gcc-multilib
 
           # 64-bit Linux
           - target: x86_64-unknown-linux-gnu
-            rust: 1.57.0 # MSRV
+            rust: 1.61.0 # MSRV
           - target: x86_64-unknown-linux-gnu
             rust: stable
     steps:
@@ -122,7 +122,7 @@ jobs:
       - uses: actions/checkout@v3
       - uses: actions-rs/toolchain@v1
         with:
-          toolchain: 1.57.0
+          toolchain: 1.61.0
           components: clippy
           override: true
           profile: minimal

Cargo.toml

@@ -14,7 +14,7 @@ keywords = ["arbitrary", "crypto", "bignum", "integer", "precision"]
 readme = "README.md"
 resolver = "2"
 edition = "2021"
-rust-version = "1.57"
+rust-version = "1.61"
 
 [dependencies]
 subtle = { version = "2.4", default-features = false }

src/lib.rs

@@ -47,6 +47,15 @@
 //! pub const MODULUS_SHR1: U256 = MODULUS.shr_vartime(1);
 //! ```
 //!
+//! Note that large constant computations may accidentally trigger a the `const_eval_limit` of the compiler.
+//! The current way to deal with this problem is to either simplify this computation,
+//! or increase the compiler's limit (currently a nightly feature).
+//! One can completely remove the compiler's limit using:
+//! ```ignore
+//! #![feature(const_eval_limit)]
+//! #![const_eval_limit = "0"]
+//! ```
+//!
 //! ### Trait-based usage
 //!
 //! The [`UInt`] type itself does not implement the standard arithmetic traits
@@ -100,6 +109,10 @@
 //! assert_eq!(b, U256::ZERO);
 //! ```
 //!
+//! It also supports modular arithmetic over constant moduli using `Residue`.
+//! That includes modular exponentiation and multiplicative inverses.
+//! These features are described in the [`modular`] module.
+//!
 //! ### Random number generation
 //!
 //! When the `rand_core` or `rand` features of this crate are enabled, it's

src/uint.rs

@@ -25,6 +25,7 @@ mod from;
 mod inv_mod;
 mod mul;
 mod mul_mod;
+mod neg;
 mod neg_mod;
 mod resize;
 mod shl;
@@ -33,6 +34,9 @@ mod sqrt;
 mod sub;
 mod sub_mod;
 
+/// Implements modular arithmetic for constant moduli.
+pub mod modular;
+
 #[cfg(feature = "generic-array")]
 mod array;
 

src/uint/add.rs

@@ -1,6 +1,6 @@
 //! [`UInt`] addition operations.
 
-use crate::{Checked, CheckedAdd, Limb, UInt, Wrapping, Zero};
+use crate::{Checked, CheckedAdd, Limb, UInt, Word, Wrapping, Zero};
 use core::ops::{Add, AddAssign};
 use subtle::CtOption;
 
@@ -36,6 +36,14 @@ impl<const LIMBS: usize> UInt<LIMBS> {
     pub const fn wrapping_add(&self, rhs: &Self) -> Self {
         self.adc(rhs, Limb::ZERO).0
     }
+
+    /// Perform wrapping addition, returning the overflow bit as a `Word` that is either 0...0 or 1...1.
+    pub(crate) const fn conditional_wrapping_add(&self, rhs: &Self, choice: Word) -> (Self, Word) {
+        let actual_rhs = UInt::ct_select(UInt::ZERO, *rhs, choice);
+        let (sum, carry) = self.adc(&actual_rhs, Limb::ZERO);
+
+        (sum, carry.0.wrapping_mul(Word::MAX))
+    }
 }
 
 impl<const LIMBS: usize> CheckedAdd<&UInt<LIMBS>> for UInt<LIMBS> {

src/uint/bit_and.rs

@@ -46,6 +46,7 @@ impl<const LIMBS: usize> BitAnd for UInt<LIMBS> {
 impl<const LIMBS: usize> BitAnd<&UInt<LIMBS>> for UInt<LIMBS> {
     type Output = UInt<LIMBS>;
 
+    #[allow(clippy::needless_borrow)]
     fn bitand(self, rhs: &UInt<LIMBS>) -> UInt<LIMBS> {
         (&self).bitand(rhs)
     }

src/uint/bit_not.rs

@@ -23,6 +23,7 @@ impl<const LIMBS: usize> UInt<LIMBS> {
 impl<const LIMBS: usize> Not for UInt<LIMBS> {
     type Output = Self;
 
+    #[allow(clippy::needless_borrow)]
     fn not(self) -> <Self as Not>::Output {
         (&self).not()
     }

src/uint/bit_or.rs

@@ -46,6 +46,7 @@ impl<const LIMBS: usize> BitOr for UInt<LIMBS> {
 impl<const LIMBS: usize> BitOr<&UInt<LIMBS>> for UInt<LIMBS> {
     type Output = UInt<LIMBS>;
 
+    #[allow(clippy::needless_borrow)]
     fn bitor(self, rhs: &UInt<LIMBS>) -> UInt<LIMBS> {
         (&self).bitor(rhs)
     }

src/uint/bit_xor.rs

@@ -46,6 +46,7 @@ impl<const LIMBS: usize> BitXor for UInt<LIMBS> {
 impl<const LIMBS: usize> BitXor<&UInt<LIMBS>> for UInt<LIMBS> {
     type Output = UInt<LIMBS>;
 
+    #[allow(clippy::needless_borrow)]
     fn bitxor(self, rhs: &UInt<LIMBS>) -> UInt<LIMBS> {
         (&self).bitxor(rhs)
     }

src/uint/cmp.rs

@@ -24,6 +24,14 @@ impl<const LIMBS: usize> UInt<LIMBS> {
         UInt { limbs }
     }
 
+    #[inline]
+    pub(crate) const fn ct_swap(a: UInt<LIMBS>, b: UInt<LIMBS>, c: Word) -> (Self, Self) {
+        let new_a = Self::ct_select(a, b, c);
+        let new_b = Self::ct_select(b, a, c);
+
+        (new_a, new_b)
+    }
+
     /// Returns all 1's if `self`!=0 or 0 if `self`==0.
     ///
     /// Const-friendly: we can't yet use `subtle` in `const fn` contexts.
@@ -38,6 +46,10 @@ impl<const LIMBS: usize> UInt<LIMBS> {
         Limb::is_nonzero(Limb(b))
     }
 
+    pub(crate) const fn ct_is_odd(&self) -> Word {
+        (self.limbs[0].0 & 1).wrapping_mul(Word::MAX)
+    }
+
     /// Returns -1 if self < rhs
     ///          0 if self == rhs
     ///          1 if self > rhs

src/uint/div.rs

@@ -68,6 +68,43 @@ impl<const LIMBS: usize> UInt<LIMBS> {
         (rem, (is_some & 1) as u8)
     }
 
+    /// Computes `self` % `rhs`, returns the remainder and
+    /// and 1 for is_some or 0 for is_none. The results can be wrapped in [`CtOption`].
+    /// NOTE: Use only if you need to access const fn. Otherwise use `reduce`
+    /// This is variable only with respect to `rhs`.
+    ///
+    /// When used with a fixed `rhs`, this function is constant-time with respect
+    /// to `self`.
+    #[allow(dead_code)]
+    pub(crate) const fn ct_reduce_wide(lower_upper: (Self, Self), rhs: &Self) -> (Self, u8) {
+        let mb = rhs.bits_vartime();
+
+        // The number of bits to consider is two sets of limbs * BIT_SIZE - mb (modulus bitcount)
+        let mut bd = (2 * LIMBS * Limb::BIT_SIZE) - mb;
+
+        // The wide integer to reduce, split into two halves
+        let (mut lower, mut upper) = lower_upper;
+
+        // Factor of the modulus, split into two halves
+        let mut c = Self::shl_vartime_wide((*rhs, UInt::ZERO), bd);
+
+        loop {
+            let (lower_sub, borrow) = lower.sbb(&c.0, Limb::ZERO);
+            let (upper_sub, borrow) = upper.sbb(&c.1, borrow);
+
+            lower = Self::ct_select(lower_sub, lower, borrow.0);
+            upper = Self::ct_select(upper_sub, upper, borrow.0);
+            if bd == 0 {
+                break;
+            }
+            bd -= 1;
+            c = Self::shr_vartime_wide(c, 1);
+        }
+
+        let is_some = Limb(mb as Word).is_nonzero();
+        (lower, (is_some & 1) as u8)
+    }
+
     /// Computes `self` % 2^k. Faster than reduce since its a power of 2.
     /// Limited to 2^16-1 since UInt doesn't support higher.
     pub const fn reduce2k(&self, k: usize) -> Self {
@@ -466,6 +503,19 @@ mod tests {
         assert_eq!(r, U256::from(3u8));
     }
 
+    #[test]
+    fn reduce_tests_wide_zero_padded() {
+        let (r, is_some) = U256::ct_reduce_wide((U256::from(10u8), U256::ZERO), &U256::from(2u8));
+        assert_eq!(is_some, 1);
+        assert_eq!(r, U256::ZERO);
+        let (r, is_some) = U256::ct_reduce_wide((U256::from(10u8), U256::ZERO), &U256::from(3u8));
+        assert_eq!(is_some, 1);
+        assert_eq!(r, U256::ONE);
+        let (r, is_some) = U256::ct_reduce_wide((U256::from(10u8), U256::ZERO), &U256::from(7u8));
+        assert_eq!(is_some, 1);
+        assert_eq!(r, U256::from(3u8));
+    }
+
     #[test]
     fn reduce_max() {
         let mut a = U256::ZERO;

src/uint/inv_mod.rs

@@ -1,5 +1,7 @@
+use subtle::{Choice, CtOption};
+
 use super::UInt;
-use crate::Limb;
+use crate::{Limb, Word};
 
 impl<const LIMBS: usize> UInt<LIMBS> {
     /// Computes 1/`self` mod 2^k as specified in Algorithm 4 from
@@ -25,11 +27,72 @@ impl<const LIMBS: usize> UInt<LIMBS> {
         }
         x
     }
+
+    /// Computes the multiplicative inverse of `self` mod `modulus`. In other words `self^-1 mod modulus`. Returns `(inverse, 1...1)` if an inverse exists, otherwise `(undefined, 0...0)`. The algorithm is the same as in GMP 6.2.1's `mpn_sec_invert`.
+    pub const fn inv_odd_mod(self, modulus: UInt<LIMBS>) -> (Self, Word) {
+        debug_assert!(modulus.ct_is_odd() == Word::MAX);
+
+        let mut a = self;
+
+        let mut u = UInt::ONE;
+        let mut v = UInt::ZERO;
+
+        let mut b = modulus;
+
+        // TODO: This can be lower if `self` is known to be small.
+        let bit_size = 2 * LIMBS * 64;
+
+        let mut m1hp = modulus;
+        let (m1hp_new, carry) = m1hp.shr_1();
+        debug_assert!(carry == Word::MAX);
+        m1hp = m1hp_new.wrapping_add(&UInt::ONE);
+
+        let mut i = 0;
+        while i < bit_size {
+            debug_assert!(b.ct_is_odd() == Word::MAX);
+
+            let self_odd = a.ct_is_odd();
+
+            // Set `self -= b` if `self` is odd.
+            let (new_a, swap) = a.conditional_wrapping_sub(&b, self_odd);
+            // Set `b += self` if `swap` is true.
+            b = UInt::ct_select(b, b.wrapping_add(&new_a), swap);
+            // Negate `self` if `swap` is true.
+            a = new_a.conditional_wrapping_neg(swap);
+
+            let (new_u, new_v) = UInt::ct_swap(u, v, swap);
+            let (new_u, cy) = new_u.conditional_wrapping_sub(&new_v, self_odd);
+            let (new_u, cyy) = new_u.conditional_wrapping_add(&modulus, cy);
+            debug_assert!(cy == cyy);
+
+            let (new_a, overflow) = a.shr_1();
+            debug_assert!(overflow == 0);
+            let (new_u, cy) = new_u.shr_1();
+            let (new_u, cy) = new_u.conditional_wrapping_add(&m1hp, cy);
+            debug_assert!(cy == 0);
+
+            a = new_a;
+            u = new_u;
+            v = new_v;
+
+            i += 1;
+        }
+
+        debug_assert!(a.ct_cmp(&UInt::ZERO) == 0);
+
+        (v, b.ct_not_eq(&UInt::ONE) ^ Word::MAX)
+    }
+
+    /// Computes the multiplicative inverse of `self` mod `modulus`. In other words `self^-1 mod modulus`. Returns `None` if the inverse does not exist. The algorithm is the same as in GMP 6.2.1's `mpn_sec_invert`.
+    pub fn inv_odd_mod_option(self, modulus: UInt<LIMBS>) -> CtOption<Self> {
+        let (inverse, exists) = self.inv_odd_mod(modulus);
+        CtOption::new(inverse, Choice::from((exists == Word::MAX) as u8))
+    }
 }
 
 #[cfg(test)]
 mod tests {
-    use crate::U256;
+    use crate::{U1024, U256, U64};
 
     #[test]
     fn inv_mod2k() {
@@ -59,4 +122,35 @@ mod tests {
         let a = v.inv_mod2k(256);
         assert_eq!(e, a);
     }
+
+    #[test]
+    fn test_invert() {
+        let a = U1024::from_be_hex("000225E99153B467A5B451979A3F451DAEF3BF8D6C6521D2FA24BBB17F29544E347A412B065B75A351EA9719E2430D2477B11CC9CF9C1AD6EDEE26CB15F463F8BCC72EF87EA30288E95A48AA792226CEC959DCB0672D8F9D80A54CBBEA85CAD8382EC224DEB2F5784E62D0CC2F81C2E6AD14EBABE646D6764B30C32B87688985");
+        let m = U1024::from_be_hex("D509E7854ABDC81921F669F1DC6F61359523F3949803E58ED4EA8BC16483DC6F37BFE27A9AC9EEA2969B357ABC5C0EE214BE16A7D4C58FC620D5B5A20AFF001AD198D3155E5799DC4EA76652D64983A7E130B5EACEBAC768D28D589C36EC749C558D0B64E37CD0775C0D0104AE7D98BA23C815185DD43CD8B16292FD94156767");
+
+        let res = a.inv_odd_mod_option(m);
+
+        let expected = U1024::from_be_hex("B03623284B0EBABCABD5C5881893320281460C0A8E7BF4BFDCFFCBCCBF436A55D364235C8171E46C7D21AAD0680676E57274A8FDA6D12768EF961CACDD2DAE5788D93DA5EB8EDC391EE3726CDCF4613C539F7D23E8702200CB31B5ED5B06E5CA3E520968399B4017BF98A864FABA2B647EFC4998B56774D4F2CB026BC024A336");
+        assert_eq!(res.unwrap(), expected);
+    }
+
+    #[test]
+    fn test_invert_small() {
+        let a = U64::from(3u64);
+        let m = U64::from(13u64);
+
+        let res = a.inv_odd_mod_option(m);
+
+        assert_eq!(U64::from(9u64), res.unwrap());
+    }
+
+    #[test]
+    fn test_no_inverse_small() {
+        let a = U64::from(14u64);
+        let m = U64::from(49u64);
+
+        let res = a.inv_odd_mod_option(m);
+
+        assert!(res.is_none().unwrap_u8() == 1);
+    }
 }

src/uint/modular/add.rs

@@ -0,0 +1,14 @@
+use crate::UInt;
+
+pub trait AddResidue {
+    /// Computes the (reduced) sum of two residues.
+    fn add(&self, rhs: &Self) -> Self;
+}
+
+pub(crate) const fn add_montgomery_form<const LIMBS: usize>(
+    a: &UInt<LIMBS>,
+    b: &UInt<LIMBS>,
+    modulus: &UInt<LIMBS>,
+) -> UInt<LIMBS> {
+    a.add_mod(b, modulus)
+}