|
| 1 | +//! Jacobi symbol calculation. |
| 2 | +
|
| 3 | +use crate::{Integer, Limb, NonZero, Odd, Word}; |
| 4 | + |
| 5 | +/// Possible values of Jacobi symbol. |
| 6 | +#[derive(Copy, Clone, Debug, PartialEq, Eq)] |
| 7 | +pub enum JacobiSymbol { |
| 8 | + /// `0` |
| 9 | + Zero, |
| 10 | + /// `1` |
| 11 | + One, |
| 12 | + /// `-1` |
| 13 | + MinusOne, |
| 14 | +} |
| 15 | + |
| 16 | +impl core::ops::Neg for JacobiSymbol { |
| 17 | + type Output = Self; |
| 18 | + fn neg(self) -> Self { |
| 19 | + match self { |
| 20 | + Self::Zero => Self::Zero, |
| 21 | + Self::One => Self::MinusOne, |
| 22 | + Self::MinusOne => Self::One, |
| 23 | + } |
| 24 | + } |
| 25 | +} |
| 26 | + |
| 27 | +// A helper trait to generalize some functions over Word and Uint. |
| 28 | +trait SmallMod { |
| 29 | + fn mod8(&self) -> Word; |
| 30 | + fn mod4(&self) -> Word; |
| 31 | +} |
| 32 | + |
| 33 | +impl SmallMod for Word { |
| 34 | + fn mod8(&self) -> Word { |
| 35 | + self & 7 |
| 36 | + } |
| 37 | + fn mod4(&self) -> Word { |
| 38 | + self & 3 |
| 39 | + } |
| 40 | +} |
| 41 | + |
| 42 | +impl<T: Integer> SmallMod for T { |
| 43 | + fn mod8(&self) -> Word { |
| 44 | + self.as_ref()[0].0 & 7 |
| 45 | + } |
| 46 | + fn mod4(&self) -> Word { |
| 47 | + self.as_ref()[0].0 & 3 |
| 48 | + } |
| 49 | +} |
| 50 | + |
| 51 | +/// Transforms `(a/p)` -> `(r/p)` for odd `p`, where the resulting `r` is odd, and `a = r * 2^s`. |
| 52 | +/// Takes a Jacobi symbol value, and returns `r` and the new Jacobi symbol, |
| 53 | +/// negated if the transformation changes parity. |
| 54 | +/// |
| 55 | +/// Note that the returned `r` is odd. |
| 56 | +fn reduce_numerator<V: SmallMod>(j: JacobiSymbol, a: Word, p: &V) -> (JacobiSymbol, Word) { |
| 57 | + let p_mod_8 = p.mod8(); |
| 58 | + let s = a.trailing_zeros(); |
| 59 | + let j = if (s & 1) == 1 && (p_mod_8 == 3 || p_mod_8 == 5) { |
| 60 | + -j |
| 61 | + } else { |
| 62 | + j |
| 63 | + }; |
| 64 | + (j, a >> s) |
| 65 | +} |
| 66 | + |
| 67 | +/// Transforms `(a/p)` -> `(p/a)` for odd and coprime `a` and `p`. |
| 68 | +/// Takes a Jacobi symbol value, and returns the swapped pair and the new Jacobi symbol, |
| 69 | +/// negated if the transformation changes parity. |
| 70 | +fn swap<T: SmallMod, V: SmallMod>(j: JacobiSymbol, a: T, p: V) -> (JacobiSymbol, V, T) { |
| 71 | + let j = if a.mod4() == 1 || p.mod4() == 1 { |
| 72 | + j |
| 73 | + } else { |
| 74 | + -j |
| 75 | + }; |
| 76 | + (j, p, a) |
| 77 | +} |
| 78 | + |
| 79 | +/// Returns the Jacobi symbol `(a/p)` given an odd `p`. Panics on even `p`. |
| 80 | +pub fn jacobi_symbol_vartime<T: Integer>(a: i32, p_long: &Odd<T>) -> JacobiSymbol { |
| 81 | + let p_long = p_long.0.clone(); |
| 82 | + |
| 83 | + let result = JacobiSymbol::One; // Keep track of all the sign flips here. |
| 84 | + |
| 85 | + // Deal with a negative `a` first: |
| 86 | + // (-a/n) = (-1/n) * (a/n) |
| 87 | + // = (-1)^((n-1)/2) * (a/n) |
| 88 | + // = (-1 if n = 3 mod 4 else 1) * (a/n) |
| 89 | + let (result, a_pos) = { |
| 90 | + let result = if a < 0 && p_long.mod4() == 3 { |
| 91 | + -result |
| 92 | + } else { |
| 93 | + result |
| 94 | + }; |
| 95 | + (result, a.abs_diff(0)) |
| 96 | + }; |
| 97 | + |
| 98 | + // A degenerate case. |
| 99 | + if a_pos == 1 || p_long == T::one() { |
| 100 | + return result; |
| 101 | + } |
| 102 | + |
| 103 | + let a_limb = Limb::from(a_pos); |
| 104 | + |
| 105 | + // Normalize input: at the end we want `a < p`, `p` odd, and both fitting into a `Word`. |
| 106 | + let (result, a, p): (JacobiSymbol, Word, Word) = if p_long.bits_vartime() <= Limb::BITS { |
| 107 | + let a = a_limb.0; |
| 108 | + let p = p_long.as_ref()[0].0; |
| 109 | + (result, a % p, p) |
| 110 | + } else { |
| 111 | + let (result, a) = reduce_numerator(result, a_limb.0, &p_long); |
| 112 | + if a == 1 { |
| 113 | + return result; |
| 114 | + } |
| 115 | + let (result, a_long, p) = swap(result, a, p_long); |
| 116 | + // Can unwrap here, since `p` is swapped with `a`, |
| 117 | + // and `a` would be odd after `reduce_numerator()`. |
| 118 | + let a = a_long.rem_vartime(&NonZero::new(T::from(p)).unwrap()); |
| 119 | + (result, a.as_ref()[0].0, p) |
| 120 | + }; |
| 121 | + |
| 122 | + let mut result = result; |
| 123 | + let mut a = a; |
| 124 | + let mut p = p; |
| 125 | + |
| 126 | + loop { |
| 127 | + if a == 0 { |
| 128 | + return JacobiSymbol::Zero; |
| 129 | + } |
| 130 | + |
| 131 | + // At this point `p` is odd (either coming from outside of the `loop`, |
| 132 | + // or from the previous iteration, where a previously reduced `a` |
| 133 | + // was swapped into its place), so we can call this. |
| 134 | + (result, a) = reduce_numerator(result, a, &p); |
| 135 | + |
| 136 | + if a == 1 { |
| 137 | + return result; |
| 138 | + } |
| 139 | + |
| 140 | + // At this point both `a` and `p` are odd: `p` was odd before, |
| 141 | + // and `a` is odd after `reduce_numerator()`. |
| 142 | + // Note that technically `swap()` only returns a valid `result` if `a` and `p` are coprime. |
| 143 | + // But if they are not, we will return `Zero` eventually, |
| 144 | + // which is not affected by any sign changes. |
| 145 | + (result, a, p) = swap(result, a, p); |
| 146 | + |
| 147 | + a %= p; |
| 148 | + } |
| 149 | +} |
| 150 | + |
| 151 | +#[cfg(test)] |
| 152 | +mod tests { |
| 153 | + use crate::{Odd, U128}; |
| 154 | + use num_bigint::{BigInt, Sign}; |
| 155 | + use num_modular::ModularSymbols; |
| 156 | + use proptest::prelude::*; |
| 157 | + |
| 158 | + use super::{jacobi_symbol_vartime, JacobiSymbol}; |
| 159 | + |
| 160 | + #[test] |
| 161 | + fn jacobi_symbol_neg_zero() { |
| 162 | + // This does not happen during normal operation, since we return zero as soon as we get it. |
| 163 | + // So just covering it for the completness' sake. |
| 164 | + assert_eq!(-JacobiSymbol::Zero, JacobiSymbol::Zero); |
| 165 | + } |
| 166 | + |
| 167 | + // Reference from `num-modular` - supports long `p`, but only positive `a`. |
| 168 | + fn jacobi_symbol_ref(a: i32, p: &U128) -> JacobiSymbol { |
| 169 | + let a_bi = BigInt::from(a); |
| 170 | + let p_bi = BigInt::from_bytes_be(Sign::Plus, p.to_be_bytes().as_ref()); |
| 171 | + let j = a_bi.jacobi(&p_bi); |
| 172 | + if j == 1 { |
| 173 | + JacobiSymbol::One |
| 174 | + } else if j == -1 { |
| 175 | + JacobiSymbol::MinusOne |
| 176 | + } else { |
| 177 | + JacobiSymbol::Zero |
| 178 | + } |
| 179 | + } |
| 180 | + |
| 181 | + #[test] |
| 182 | + fn small_values() { |
| 183 | + // Test small values, using a reference implementation. |
| 184 | + for a in -31i32..31 { |
| 185 | + for p in (1u32..31).step_by(2) { |
| 186 | + let p_long = U128::from(p); |
| 187 | + let j_ref = jacobi_symbol_ref(a, &p_long); |
| 188 | + let j = jacobi_symbol_vartime(a, &Odd::new(p_long).unwrap()); |
| 189 | + assert_eq!(j, j_ref); |
| 190 | + } |
| 191 | + } |
| 192 | + } |
| 193 | + |
| 194 | + #[test] |
| 195 | + fn big_values() { |
| 196 | + // a = x, p = x * y, where x and y are big primes. Should give 0. |
| 197 | + let a = 2147483647i32; // 2^31 - 1, a prime |
| 198 | + let p = U128::from_be_hex("000000007ffffffeffffffe28000003b"); // (2^31 - 1) * (2^64 - 59) |
| 199 | + assert_eq!( |
| 200 | + jacobi_symbol_vartime(a, &Odd::new(p).unwrap()), |
| 201 | + JacobiSymbol::Zero |
| 202 | + ); |
| 203 | + assert_eq!(jacobi_symbol_ref(a, &p), JacobiSymbol::Zero); |
| 204 | + |
| 205 | + // a = x^2 mod p, should give 1. |
| 206 | + let a = 659456i32; // Obtained from x = 2^70 |
| 207 | + let p = U128::from_be_hex("ffffffffffffffffffffffffffffff5f"); // 2^128 - 161 - not a prime |
| 208 | + assert_eq!( |
| 209 | + jacobi_symbol_vartime(a, &Odd::new(p).unwrap()), |
| 210 | + JacobiSymbol::One |
| 211 | + ); |
| 212 | + assert_eq!(jacobi_symbol_ref(a, &p), JacobiSymbol::One); |
| 213 | + |
| 214 | + let a = i32::MIN; // -2^31, check that no overflow occurs |
| 215 | + let p = U128::from_be_hex("000000007ffffffeffffffe28000003b"); // (2^31 - 1) * (2^64 - 59) |
| 216 | + assert_eq!( |
| 217 | + jacobi_symbol_vartime(a, &Odd::new(p).unwrap()), |
| 218 | + JacobiSymbol::One |
| 219 | + ); |
| 220 | + assert_eq!(jacobi_symbol_ref(a, &p), JacobiSymbol::One); |
| 221 | + } |
| 222 | + |
| 223 | + prop_compose! { |
| 224 | + fn odd_uint()(bytes in any::<[u8; 16]>()) -> Odd<U128> { |
| 225 | + Odd::new(U128::from_le_slice(&bytes) | U128::ONE).unwrap() |
| 226 | + } |
| 227 | + } |
| 228 | + |
| 229 | + proptest! { |
| 230 | + #[test] |
| 231 | + fn fuzzy(a in any::<i32>(), p in odd_uint()) { |
| 232 | + let j_ref = jacobi_symbol_ref(a, &p); |
| 233 | + let j = jacobi_symbol_vartime(a, &p); |
| 234 | + assert_eq!(j, j_ref); |
| 235 | + } |
| 236 | + } |
| 237 | +} |
0 commit comments