Add vartime GCD and Jacobi symbol functions

src/gcd.rs

@@ -0,0 +1,132 @@
+use crate::{Integer, Limb, NonZero, Word};
+
+/// Compute the greatest common divisor (GCD) of this number and another.
+/// Variable time in both arguments.
+pub fn gcd_vartime<T: Integer>(lhs: &T, rhs: &T) -> T {
+    // This we can check since it doesn't affect the return type,
+    // even though `n` will not be 0 either in the application.
+    if lhs.is_zero().into() {
+        return rhs.clone();
+    }
+
+    if rhs.is_zero().into() {
+        return lhs.clone();
+    }
+
+    // Normalize input: the resulting (a, b) are both small, a >= b, and b != 0.
+
+    let mut a = lhs.clone();
+    let mut b = rhs.clone();
+
+    let mut a_ref = &mut a;
+    let mut b_ref = &mut b;
+
+    if a_ref > b_ref {
+        (a_ref, b_ref) = (b_ref, a_ref);
+    }
+
+    // Euclidean algorithm.
+    loop {
+        let r = a_ref.rem_vartime(&NonZero::new(b_ref.clone()).unwrap());
+        if r.is_zero().into() {
+            return b_ref.clone();
+        }
+
+        (a_ref, b_ref) = (b_ref, a_ref);
+        *b_ref = r;
+    }
+}
+
+/// Compute the greatest common divisor (GCD) of this number and another.
+/// Variable time in both arguments.
+pub fn gcd_limb_vartime<T: Integer>(lhs: &T, rhs: NonZero<Limb>) -> Limb {
+    // Allowing `rhs = 0` would require us to have the return type of `T`
+    // (since `gcd(n, 0) = n`).
+
+    // This we can check since it doesn't affect the return type,
+    // even though `n` will not be 0 either in the application.
+    if lhs.is_zero().into() {
+        return rhs.0;
+    }
+
+    // Normalize input: the resulting (a, b) are both small, `a >= b`, and `b != 0`.
+    let (mut a, mut b): (Word, Word) = if lhs.bits_vartime() > Limb::BITS {
+        let n = lhs.rem_vartime(&NonZero::new(T::from(rhs.0)).unwrap());
+        (rhs.0 .0, n.as_ref()[0].0)
+    } else {
+        // In this branch `lhs` fits in a `Limb`,
+        // so we can safely take the first limb and use it.
+        let n = lhs.as_ref()[0];
+        if n > rhs.0 {
+            (n.0, rhs.0 .0)
+        } else {
+            (rhs.0 .0, n.0)
+        }
+    };
+
+    // Euclidean algorithm.
+    // Binary GCD algorithm could be used here,
+    // but the performance impact of this code is negligible.
+    loop {
+        let r = a % b;
+        if r == 0 {
+            return Limb(b);
+        }
+        (a, b) = (b, r)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::{Limb, NonZero, Word, U128};
+    use num_bigint::BigUint;
+    use num_integer::Integer;
+    use proptest::prelude::*;
+
+    use super::gcd_limb_vartime;
+
+    #[test]
+    fn corner_cases() {
+        assert_eq!(
+            gcd_limb_vartime(&U128::from(0u64), NonZero::new(Limb(5)).unwrap()),
+            Limb(5)
+        );
+        assert_eq!(
+            gcd_limb_vartime(&U128::from(1u64), NonZero::new(Limb(11 * 13 * 19)).unwrap()),
+            Limb(1)
+        );
+        assert_eq!(
+            gcd_limb_vartime(&U128::from(7u64 * 11 * 13), NonZero::new(Limb(1)).unwrap()),
+            Limb(1)
+        );
+        assert_eq!(
+            gcd_limb_vartime(
+                &U128::from(7u64 * 11 * 13),
+                NonZero::new(Limb(11 * 13 * 19)).unwrap()
+            ),
+            Limb(11 * 13)
+        );
+    }
+
+    prop_compose! {
+        fn uint()(bytes in any::<[u8; 16]>()) -> U128 {
+            U128::from_le_slice(&bytes) | U128::ONE
+        }
+    }
+
+    proptest! {
+        #[test]
+        fn fuzzy(m in any::<Word>(), n in uint()) {
+            if m == 0 {
+                return Ok(());
+            }
+
+            let m_bi = BigUint::from(m);
+            let n_bi = BigUint::from_bytes_be(n.to_be_bytes().as_ref());
+            let gcd_ref: Word = n_bi.gcd(&m_bi).try_into().unwrap();
+
+            let gcd_test = gcd_limb_vartime(&n, NonZero::new(Limb(m)).unwrap()).0;
+            assert_eq!(gcd_test, gcd_ref);
+        }
+    }
+}

src/jacobi.rs

@@ -0,0 +1,237 @@
+//! Jacobi symbol calculation.
+
+use crate::{Integer, Limb, NonZero, Odd, Word};
+
+/// Possible values of Jacobi symbol.
+#[derive(Copy, Clone, Debug, PartialEq, Eq)]
+pub enum JacobiSymbol {
+    /// `0`
+    Zero,
+    /// `1`
+    One,
+    /// `-1`
+    MinusOne,
+}
+
+impl core::ops::Neg for JacobiSymbol {
+    type Output = Self;
+    fn neg(self) -> Self {
+        match self {
+            Self::Zero => Self::Zero,
+            Self::One => Self::MinusOne,
+            Self::MinusOne => Self::One,
+        }
+    }
+}
+
+// A helper trait to generalize some functions over Word and Uint.
+trait SmallMod {
+    fn mod8(&self) -> Word;
+    fn mod4(&self) -> Word;
+}
+
+impl SmallMod for Word {
+    fn mod8(&self) -> Word {
+        self & 7
+    }
+    fn mod4(&self) -> Word {
+        self & 3
+    }
+}
+
+impl<T: Integer> SmallMod for T {
+    fn mod8(&self) -> Word {
+        self.as_ref()[0].0 & 7
+    }
+    fn mod4(&self) -> Word {
+        self.as_ref()[0].0 & 3
+    }
+}
+
+/// Transforms `(a/p)` -> `(r/p)` for odd `p`, where the resulting `r` is odd, and `a = r * 2^s`.
+/// Takes a Jacobi symbol value, and returns `r` and the new Jacobi symbol,
+/// negated if the transformation changes parity.
+///
+/// Note that the returned `r` is odd.
+fn reduce_numerator<V: SmallMod>(j: JacobiSymbol, a: Word, p: &V) -> (JacobiSymbol, Word) {
+    let p_mod_8 = p.mod8();
+    let s = a.trailing_zeros();
+    let j = if (s & 1) == 1 && (p_mod_8 == 3 || p_mod_8 == 5) {
+        -j
+    } else {
+        j
+    };
+    (j, a >> s)
+}
+
+/// Transforms `(a/p)` -> `(p/a)` for odd and coprime `a` and `p`.
+/// Takes a Jacobi symbol value, and returns the swapped pair and the new Jacobi symbol,
+/// negated if the transformation changes parity.
+fn swap<T: SmallMod, V: SmallMod>(j: JacobiSymbol, a: T, p: V) -> (JacobiSymbol, V, T) {
+    let j = if a.mod4() == 1 || p.mod4() == 1 {
+        j
+    } else {
+        -j
+    };
+    (j, p, a)
+}
+
+/// Returns the Jacobi symbol `(a/p)` given an odd `p`. Panics on even `p`.
+pub fn jacobi_symbol_vartime<T: Integer>(a: i32, p_long: &Odd<T>) -> JacobiSymbol {
+    let p_long = p_long.0.clone();
+
+    let result = JacobiSymbol::One; // Keep track of all the sign flips here.
+
+    // Deal with a negative `a` first:
+    // (-a/n) = (-1/n) * (a/n)
+    //        = (-1)^((n-1)/2) * (a/n)
+    //        = (-1 if n = 3 mod 4 else 1) * (a/n)
+    let (result, a_pos) = {
+        let result = if a < 0 && p_long.mod4() == 3 {
+            -result
+        } else {
+            result
+        };
+        (result, a.abs_diff(0))
+    };
+
+    // A degenerate case.
+    if a_pos == 1 || p_long == T::one() {
+        return result;
+    }
+
+    let a_limb = Limb::from(a_pos);
+
+    // Normalize input: at the end we want `a < p`, `p` odd, and both fitting into a `Word`.
+    let (result, a, p): (JacobiSymbol, Word, Word) = if p_long.bits_vartime() <= Limb::BITS {
+        let a = a_limb.0;
+        let p = p_long.as_ref()[0].0;
+        (result, a % p, p)
+    } else {
+        let (result, a) = reduce_numerator(result, a_limb.0, &p_long);
+        if a == 1 {
+            return result;
+        }
+        let (result, a_long, p) = swap(result, a, p_long);
+        // Can unwrap here, since `p` is swapped with `a`,
+        // and `a` would be odd after `reduce_numerator()`.
+        let a = a_long.rem_vartime(&NonZero::new(T::from(p)).unwrap());
+        (result, a.as_ref()[0].0, p)
+    };
+
+    let mut result = result;
+    let mut a = a;
+    let mut p = p;
+
+    loop {
+        if a == 0 {
+            return JacobiSymbol::Zero;
+        }
+
+        // At this point `p` is odd (either coming from outside of the `loop`,
+        // or from the previous iteration, where a previously reduced `a`
+        // was swapped into its place), so we can call this.
+        (result, a) = reduce_numerator(result, a, &p);
+
+        if a == 1 {
+            return result;
+        }
+
+        // At this point both `a` and `p` are odd: `p` was odd before,
+        // and `a` is odd after `reduce_numerator()`.
+        // Note that technically `swap()` only returns a valid `result` if `a` and `p` are coprime.
+        // But if they are not, we will return `Zero` eventually,
+        // which is not affected by any sign changes.
+        (result, a, p) = swap(result, a, p);
+
+        a %= p;
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::{Odd, U128};
+    use num_bigint::{BigInt, Sign};
+    use num_modular::ModularSymbols;
+    use proptest::prelude::*;
+
+    use super::{jacobi_symbol_vartime, JacobiSymbol};
+
+    #[test]
+    fn jacobi_symbol_neg_zero() {
+        // This does not happen during normal operation, since we return zero as soon as we get it.
+        // So just covering it for the completness' sake.
+        assert_eq!(-JacobiSymbol::Zero, JacobiSymbol::Zero);
+    }
+
+    // Reference from `num-modular` - supports long `p`, but only positive `a`.
+    fn jacobi_symbol_ref(a: i32, p: &U128) -> JacobiSymbol {
+        let a_bi = BigInt::from(a);
+        let p_bi = BigInt::from_bytes_be(Sign::Plus, p.to_be_bytes().as_ref());
+        let j = a_bi.jacobi(&p_bi);
+        if j == 1 {
+            JacobiSymbol::One
+        } else if j == -1 {
+            JacobiSymbol::MinusOne
+        } else {
+            JacobiSymbol::Zero
+        }
+    }
+
+    #[test]
+    fn small_values() {
+        // Test small values, using a reference implementation.
+        for a in -31i32..31 {
+            for p in (1u32..31).step_by(2) {
+                let p_long = U128::from(p);
+                let j_ref = jacobi_symbol_ref(a, &p_long);
+                let j = jacobi_symbol_vartime(a, &Odd::new(p_long).unwrap());
+                assert_eq!(j, j_ref);
+            }
+        }
+    }
+
+    #[test]
+    fn big_values() {
+        // a = x, p = x * y, where x and y are big primes. Should give 0.
+        let a = 2147483647i32; // 2^31 - 1, a prime
+        let p = U128::from_be_hex("000000007ffffffeffffffe28000003b"); // (2^31 - 1) * (2^64 - 59)
+        assert_eq!(
+            jacobi_symbol_vartime(a, &Odd::new(p).unwrap()),
+            JacobiSymbol::Zero
+        );
+        assert_eq!(jacobi_symbol_ref(a, &p), JacobiSymbol::Zero);
+
+        // a = x^2 mod p, should give 1.
+        let a = 659456i32; // Obtained from x = 2^70
+        let p = U128::from_be_hex("ffffffffffffffffffffffffffffff5f"); // 2^128 - 161 - not a prime
+        assert_eq!(
+            jacobi_symbol_vartime(a, &Odd::new(p).unwrap()),
+            JacobiSymbol::One
+        );
+        assert_eq!(jacobi_symbol_ref(a, &p), JacobiSymbol::One);
+
+        let a = i32::MIN; // -2^31, check that no overflow occurs
+        let p = U128::from_be_hex("000000007ffffffeffffffe28000003b"); // (2^31 - 1) * (2^64 - 59)
+        assert_eq!(
+            jacobi_symbol_vartime(a, &Odd::new(p).unwrap()),
+            JacobiSymbol::One
+        );
+        assert_eq!(jacobi_symbol_ref(a, &p), JacobiSymbol::One);
+    }
+
+    prop_compose! {
+        fn odd_uint()(bytes in any::<[u8; 16]>()) -> Odd<U128> {
+            Odd::new(U128::from_le_slice(&bytes) | U128::ONE).unwrap()
+        }
+    }
+
+    proptest! {
+        #[test]
+        fn fuzzy(a in any::<i32>(), p in odd_uint()) {
+            let j_ref = jacobi_symbol_ref(a, &p);
+            let j = jacobi_symbol_vartime(a, &p);
+            assert_eq!(j, j_ref);
+        }
+    }
+}

src/lib.rs

@@ -168,6 +168,8 @@ pub mod modular;
 mod array;
 mod checked;
 mod const_choice;
+mod gcd;
+mod jacobi;
 mod limb;
 mod non_zero;
 mod odd;
@@ -179,6 +181,8 @@ mod wrapping;
 pub use crate::{
     checked::Checked,
     const_choice::{ConstChoice, ConstCtOption},
+    gcd::{gcd_limb_vartime, gcd_vartime},
+    jacobi::{jacobi_symbol_vartime, JacobiSymbol},
     limb::{Limb, WideWord, Word},
     non_zero::NonZero,
     odd::Odd,