Add Jacobi symbol calculation

benches/const_monty.rs

@@ -119,6 +119,22 @@ fn bench_montgomery_ops<M: Measurement>(group: &mut BenchmarkGroup<'_, M>) {
         )
     });
 
+    group.bench_function("jacobi_symbol", |b| {
+        b.iter_batched(
+            || ConstMontyForm::random(&mut rng),
+            |a| a.jacobi_symbol(),
+            BatchSize::SmallInput,
+        )
+    });
+
+    group.bench_function("jacobi_symbol_vartime", |b| {
+        b.iter_batched(
+            || ConstMontyForm::random(&mut rng),
+            |a| a.jacobi_symbol_vartime(),
+            BatchSize::SmallInput,
+        )
+    });
+
     group.bench_function("lincomb, U256*U256+U256*U256", |b| {
         b.iter_batched(
             || ConstMontyForm::random(&mut rng),

benches/uint.rs

@@ -468,6 +468,46 @@ fn bench_gcd(c: &mut Criterion) {
     group.finish();
 }
 
+fn mod_symbols_bench<const LIMBS: usize>(g: &mut BenchmarkGroup<WallTime>, _x: Uint<LIMBS>) {
+    let mut rng = make_rng();
+
+    g.bench_function(BenchmarkId::new("jacobi_symbol", LIMBS), |b| {
+        b.iter_batched(
+            || {
+                (
+                    OddUint::<LIMBS>::random(&mut rng),
+                    Uint::<LIMBS>::random(&mut rng),
+                )
+            },
+            |(f, g)| black_box(g.jacobi_symbol(&f)),
+            BatchSize::SmallInput,
+        )
+    });
+
+    g.bench_function(BenchmarkId::new("jacobi_symbol_vartime", LIMBS), |b| {
+        b.iter_batched(
+            || {
+                (
+                    OddUint::<LIMBS>::random(&mut rng),
+                    Uint::<LIMBS>::random(&mut rng),
+                )
+            },
+            |(f, g)| black_box(g.jacobi_symbol_vartime(&f)),
+            BatchSize::SmallInput,
+        )
+    });
+}
+
+fn bench_mod_symbols(c: &mut Criterion) {
+    let mut group = c.benchmark_group("modular symbols");
+
+    mod_symbols_bench(&mut group, Uint::<1>::ZERO);
+    mod_symbols_bench(&mut group, Uint::<4>::ZERO);
+    mod_symbols_bench(&mut group, Uint::<64>::ZERO);
+
+    group.finish();
+}
+
 fn bench_shl(c: &mut Criterion) {
     let mut group = c.benchmark_group("left shift");
 
@@ -631,6 +671,7 @@ criterion_group!(
     bench_mul,
     bench_division,
     bench_gcd,
+    bench_mod_symbols,
     bench_shl,
     bench_shr,
     bench_invert_mod,

src/const_choice.rs

@@ -157,6 +157,12 @@ impl ConstChoice {
         Self::from_u32_lsb(bit)
     }
 
+    /// Returns the truthy value if `x == y`, and the falsy value otherwise.
+    #[inline]
+    pub(crate) const fn from_i64_eq(x: i64, y: i64) -> Self {
+        Self::from_word_nonzero(x as Word ^ y as Word).not()
+    }
+
     #[inline]
     pub(crate) const fn not(&self) -> Self {
         Self(!self.0)

src/int.rs

@@ -23,6 +23,7 @@ mod encoding;
 mod from;
 mod gcd;
 mod invert_mod;
+mod mod_symbol;
 mod mul;
 mod mul_uint;
 mod neg;

src/int/mod_symbol.rs

@@ -0,0 +1,112 @@
+//! Support for computing modular symbols.
+
+use crate::{ConstChoice, Int, JacobiSymbol, Odd, Uint};
+
+impl<const LIMBS: usize> Int<LIMBS> {
+    /// Compute the Jacobi symbol `(self|rhs)`.
+    ///
+    /// For prime `rhs`, this corresponds to the Legendre symbol and
+    /// indicates whether `self` is quadratic residue modulo `rhs`.
+    pub const fn jacobi_symbol(&self, rhs: &Odd<Uint<LIMBS>>) -> JacobiSymbol {
+        let (abs, sign) = self.abs_sign();
+        let jacobi = abs.jacobi_symbol(rhs) as i64;
+        // (-self|rhs) = -(self|rhs) iff rhs = 3 mod 4
+        let swap = sign.and(ConstChoice::from_word_eq(rhs.as_ref().limbs[0].0 & 3, 3));
+        JacobiSymbol::from_i8(swap.select_i64(jacobi, -jacobi) as i8)
+    }
+
+    /// Compute the Jacobi symbol `(self|rhs)`.
+    ///
+    /// For prime `rhs`, this corresponds to the Legendre symbol and
+    /// indicates whether `self` is quadratic residue modulo `rhs`.
+    ///
+    /// This method executes in variable-time for the value of `self`.
+    pub const fn jacobi_symbol_vartime(&self, rhs: &Odd<Uint<LIMBS>>) -> JacobiSymbol {
+        let (abs, sign) = self.abs_sign();
+        let jacobi = abs.jacobi_symbol_vartime(rhs);
+        JacobiSymbol::from_i8(
+            if sign.is_true_vartime() && rhs.as_ref().limbs[0].0 & 3 == 3 {
+                -(jacobi as i8)
+            } else {
+                jacobi as i8
+            },
+        )
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::{I256, JacobiSymbol, U256};
+
+    #[test]
+    fn jacobi_quad_residue() {
+        // Two semiprimes with no common factors, and
+        // f is quadratic residue modulo g
+        let f = I256::from(59i32 * 67);
+        let g = U256::from(61u32 * 71).to_odd().unwrap();
+        let res = f.jacobi_symbol(&g);
+        let res_vartime = f.jacobi_symbol_vartime(&g);
+        assert_eq!(res, JacobiSymbol::One);
+        assert_eq!(res, res_vartime);
+    }
+
+    #[test]
+    fn jacobi_non_quad_residue() {
+        // f and g have no common factors, but
+        // f is not quadratic residue modulo g
+        let f = I256::from(59i32 * 67 + 2);
+        let g = U256::from(61u32 * 71).to_odd().unwrap();
+        let res = f.jacobi_symbol(&g);
+        let res_vartime = f.jacobi_symbol_vartime(&g);
+        assert_eq!(res, JacobiSymbol::MinusOne);
+        assert_eq!(res, res_vartime);
+    }
+
+    #[test]
+    fn jacobi_non_coprime() {
+        let f = I256::from(4391633i32);
+        let g = U256::from(2022161u32).to_odd().unwrap();
+        let res = f.jacobi_symbol(&g);
+        let res_vartime = f.jacobi_symbol_vartime(&g);
+        assert_eq!(res, JacobiSymbol::Zero);
+        assert_eq!(res, res_vartime);
+    }
+
+    #[test]
+    fn jacobi_zero() {
+        assert_eq!(
+            I256::ZERO.jacobi_symbol(&U256::ONE.to_odd().unwrap()),
+            JacobiSymbol::One
+        );
+    }
+
+    #[test]
+    fn jacobi_neg_one() {
+        let f = I256::ONE;
+        assert_eq!(
+            f.jacobi_symbol(&U256::ONE.to_odd().unwrap()),
+            JacobiSymbol::One
+        );
+        assert_eq!(
+            f.jacobi_symbol(&U256::from(3u8).to_odd().unwrap()),
+            JacobiSymbol::One
+        );
+    }
+
+    #[test]
+    fn jacobi_int() {
+        // Two semiprimes with no common factors, and
+        // f is quadratic residue modulo g
+        let f = I256::from(59i32 * 67);
+        let g = U256::from(61u32 * 71).to_odd().unwrap();
+        let res = f.jacobi_symbol(&g);
+        let res_vartime = f.jacobi_symbol_vartime(&g);
+        assert_eq!(res, JacobiSymbol::One);
+        assert_eq!(res, res_vartime);
+
+        let res = f.checked_neg().unwrap().jacobi_symbol(&g);
+        let res_vartime = f.checked_neg().unwrap().jacobi_symbol_vartime(&g);
+        assert_eq!(res, JacobiSymbol::MinusOne);
+        assert_eq!(res, res_vartime);
+    }
+}

src/modular/bingcd/compact.rs

@@ -53,8 +53,8 @@ impl<const LIMBS: usize> Uint<LIMBS> {
             .bitand(&mask)
     }
 
-    /// Compact `self` to a form containing the concatenation of its bit ranges `[0, K-1)`
-    /// and `[n-K-1, n)`.
+    /// Compact `self` to a form containing the concatenation of its bit ranges `[0, K)`
+    /// and `[n-K, n)`.
     ///
     /// Assumes `K ≤ Uint::<SUMMARY_LIMBS>::BITS`, `n ≤ Self::BITS` and `n ≥ 2K`.
     #[inline(always)]
@@ -67,13 +67,13 @@ impl<const LIMBS: usize> Uint<LIMBS> {
         debug_assert!(n >= 2 * K);
 
         // safe to vartime; this function is vartime in length only, which is a public constant
-        let hi = self.section_vartime_length(n - K - 1, K + 1);
+        let hi = self.section_vartime_length(n - K, K);
         // safe to vartime; this function is vartime in idx and length only, which are both public
         // constants
-        let lo = self.section_vartime(0, K - 1);
+        let lo = self.section_vartime(0, K);
         // safe to vartime; shl_vartime is variable in the value of shift only. Since this shift
         // is a public constant, the constant time property of this algorithm is not impacted.
-        hi.shl_vartime(K - 1).bitxor(&lo)
+        hi.shl_vartime(K).bitxor(&lo)
     }
 
     /// Vartime equivalent of [`Self::compact`].
@@ -87,13 +87,13 @@ impl<const LIMBS: usize> Uint<LIMBS> {
         debug_assert!(n >= 2 * K);
 
         // safe to vartime; this function is vartime in length only, which is a public constant
-        let hi = self.section_vartime(n - K - 1, K + 1);
+        let hi = self.section_vartime(n - K, K);
         // safe to vartime; this function is vartime in idx and length only, which are both public
         // constants
-        let lo = self.section_vartime(0, K - 1);
+        let lo = self.section_vartime(0, K);
         // safe to vartime; shl_vartime is variable in the value of shift only. Since this shift
         // is a public constant, the constant time property of this algorithm is not impacted.
-        hi.shl_vartime(K - 1).bitxor(&lo)
+        hi.shl_vartime(K).bitxor(&lo)
     }
 }
 
@@ -105,7 +105,7 @@ mod tests {
     fn test_compact() {
         let val =
             U256::from_be_hex("CFCF1535CEBE19BBF289933AB8645189397450A32BFEC57579FB7EB14E27D101");
-        let target = U128::from_be_hex("BBF289933AB86451F9FB7EB14E27D101");
+        let target = U128::from_be_hex("BBF289933AB8645179FB7EB14E27D101");
 
         let compact = val.compact::<64, { U128::LIMBS }>(200);
         assert_eq!(compact, target);