zenith_float_num/
modular.rs

1//! Modular exponentiation, inverses, Miller–Rabin, and Pollard–Brent on [`ExactInt`].
2
3use crate::ExactInt;
4use core::cmp::Ordering;
5
6/// Maximum `f` evaluations in one [`pollard_rho`] run before `None`.
7pub const POLLARD_RHO_ITER_MAX: usize = 1_048_576;
8
9/// Brent inner-loop product batch.
10const POLLARD_BRENT_M: usize = 128;
11
12fn two() -> ExactInt {
13    ExactInt::from_u64(2)
14}
15
16fn abs_int(a: &ExactInt) -> ExactInt {
17    if a.is_negative() {
18        a.neg()
19    } else {
20        a.clone()
21    }
22}
23
24fn is_even(a: &ExactInt) -> bool {
25    a.low_word() & 1 == 0
26}
27
28fn rem_nonneg(a: &ExactInt, m: &ExactInt) -> Option<ExactInt> {
29    let mabs = abs_int(m);
30    if mabs.is_zero() {
31        return None;
32    }
33    let (_, r) = a.div_rem(&mabs)?;
34    if r.is_negative() {
35        Some(r.add(&mabs))
36    } else {
37        Some(r)
38    }
39}
40
41fn mul_mod(a: &ExactInt, b: &ExactInt, m: &ExactInt) -> Option<ExactInt> {
42    rem_nonneg(&a.mul(b), m)
43}
44
45fn add_mod(a: &ExactInt, b: &ExactInt, m: &ExactInt) -> Option<ExactInt> {
46    rem_nonneg(&a.add(b), m)
47}
48
49/// `base^exp mod modulus`. `exp < 0` or a zero modulus is `None`.
50/// `|modulus| = 1` is `0`.
51pub fn mod_pow(base: &ExactInt, exp: &ExactInt, modulus: &ExactInt) -> Option<ExactInt> {
52    if exp.is_negative() {
53        return None;
54    }
55    let m = abs_int(modulus);
56    if m.is_zero() {
57        return None;
58    }
59    if m.is_one() {
60        return Some(ExactInt::zero());
61    }
62    let mut acc = ExactInt::one();
63    let mut b = rem_nonneg(base, &m)?;
64    let mut e = exp.clone();
65    let two = two();
66    while !e.is_zero() {
67        if !is_even(&e) {
68            acc = mul_mod(&acc, &b, &m)?;
69        }
70        e = e.div_rem(&two)?.0;
71        if !e.is_zero() {
72            b = mul_mod(&b, &b, &m)?;
73        }
74    }
75    Some(acc)
76}
77
78fn egcd(a: &ExactInt, b: &ExactInt) -> (ExactInt, ExactInt, ExactInt) {
79    if b.is_zero() {
80        return (a.clone(), ExactInt::one(), ExactInt::zero());
81    }
82    let (q, r) = a.div_rem(b).unwrap_or((ExactInt::zero(), ExactInt::zero()));
83    let (g, x, y) = egcd(b, &r);
84    // x1 = y, y1 = x - q y
85    (g, y.clone(), x.sub(&q.mul(&y)))
86}
87
88/// Modular inverse of `a` modulo `modulus`. `None` if not invertible.
89pub fn mod_inv(a: &ExactInt, modulus: &ExactInt) -> Option<ExactInt> {
90    let m = abs_int(modulus);
91    if m.is_zero() || m.is_one() {
92        return None;
93    }
94    let aa = rem_nonneg(a, &m)?;
95    if aa.is_zero() {
96        return None;
97    }
98    let (g, x, _) = egcd(&aa, &m);
99    if !g.is_one() {
100        return None;
101    }
102    rem_nonneg(&x, &m)
103}
104
105/// Miller–Rabin on `|n|` with the given bases. `n < 2` is `false`.
106///
107/// Deterministic for `n < 3·10^{18}` when the witness list is a complete
108/// Jaeschke set; this function only uses the callers' bases.
109pub fn miller_rabin(n: &ExactInt, witnesses: &[ExactInt]) -> bool {
110    let n = abs_int(n);
111    if n.cmp(&two()) == Ordering::Less {
112        return false;
113    }
114    if n == two() || n == ExactInt::from_u64(3) {
115        return true;
116    }
117    if is_even(&n) {
118        return false;
119    }
120    let one = ExactInt::one();
121    let n_minus = n.sub(&one);
122    let mut d = n_minus.clone();
123    let mut s = 0usize;
124    let two = two();
125    while !d.is_zero() && is_even(&d) {
126        d = d
127            .div_rem(&two)
128            .map(|(q, _)| q)
129            .unwrap_or_else(ExactInt::zero);
130        s += 1;
131    }
132    if s == 0 {
133        return false;
134    }
135    'wit: for a in witnesses {
136        let a = rem_nonneg(a, &n).unwrap_or_else(ExactInt::zero);
137        if a.is_zero() || a.is_one() {
138            continue;
139        }
140        let mut x = match mod_pow(&a, &d, &n) {
141            Some(v) => v,
142            None => return false,
143        };
144        if x.is_one() || x == n_minus {
145            continue;
146        }
147        for _ in 1..s {
148            x = match mul_mod(&x, &x, &n) {
149                Some(v) => v,
150                None => return false,
151            };
152            if x == n_minus {
153                continue 'wit;
154            }
155            if x.is_one() {
156                return false;
157            }
158        }
159        return false;
160    }
161    !witnesses.is_empty()
162}
163
164fn pollard_f(x: &ExactInt, c: &ExactInt, n: &ExactInt) -> Option<ExactInt> {
165    add_mod(&mul_mod(x, x, n)?, c, n)
166}
167
168fn brent_once(n: &ExactInt, c: &ExactInt, max_f: usize) -> Option<ExactInt> {
169    let mut y = ExactInt::zero();
170    let mut g = ExactInt::one();
171    let mut r = 1usize;
172    let mut q = ExactInt::one();
173    let mut f_used = 0usize;
174    let mut x = ExactInt::zero();
175    let mut ys = ExactInt::zero();
176    while g.is_one() {
177        x = y.clone();
178        for _ in 0..r {
179            y = pollard_f(&y, c, n)?;
180            f_used += 1;
181            if f_used >= max_f {
182                return None;
183            }
184        }
185        let mut k = 0usize;
186        while k < r && g.is_one() {
187            ys = y.clone();
188            let steps = POLLARD_BRENT_M.min(r - k);
189            for _ in 0..steps {
190                y = pollard_f(&y, c, n)?;
191                f_used += 1;
192                if f_used >= max_f {
193                    return None;
194                }
195                let diff = abs_int(&x.sub(&y));
196                q = mul_mod(&q, &diff, n)?;
197            }
198            g = q.gcd(n);
199            k += steps;
200        }
201        if r > usize::MAX / 2 {
202            return None;
203        }
204        r = r.saturating_mul(2);
205    }
206    if g == *n {
207        loop {
208            ys = pollard_f(&ys, c, n)?;
209            f_used += 1;
210            if f_used >= max_f {
211                return None;
212            }
213            g = abs_int(&x.sub(&ys)).gcd(n);
214            if !g.is_one() {
215                break;
216            }
217        }
218    }
219    if g.is_one() || g == *n {
220        None
221    } else {
222        Some(g)
223    }
224}
225
226/// Brent Pollard ρ. A proper factor of `|n|`, or `None` if prime / cap hit.
227pub fn pollard_rho(n: &ExactInt) -> Option<ExactInt> {
228    let n = abs_int(n);
229    if n.cmp(&two()) != Ordering::Greater {
230        return None;
231    }
232    if is_even(&n) {
233        return Some(two());
234    }
235    for c in 1u64..=32 {
236        if let Some(f) = brent_once(&n, &ExactInt::from_u64(c), POLLARD_RHO_ITER_MAX) {
237            if !f.is_one() && f != n {
238                return Some(f);
239            }
240        }
241    }
242    None
243}
244
245#[cfg(test)]
246mod tests {
247    use super::*;
248
249    #[test]
250    fn modular_pow_inv_miller_rho() {
251        let two = ExactInt::from_u64(2);
252        let exp = ExactInt::from_u64(100);
253        let m = ExactInt::from_u64(1_000_000_007);
254        let got = mod_pow(&two, &exp, &m).expect("mod_pow");
255        assert_eq!(got, ExactInt::from_u64(976_371_285));
256
257        let inv = mod_inv(&ExactInt::from_u64(3), &ExactInt::from_u64(7)).expect("inv");
258        assert_eq!(inv, ExactInt::from_u64(5));
259        assert!(mod_inv(&ExactInt::from_u64(2), &ExactInt::from_u64(4)).is_none());
260
261        let mersenne = two.pow(31).sub(&ExactInt::one());
262        let wits: Vec<ExactInt> = [2u64, 3, 5, 7]
263            .into_iter()
264            .map(ExactInt::from_u64)
265            .collect();
266        assert!(miller_rabin(&mersenne, &wits));
267        assert!(!miller_rabin(&ExactInt::from_u64(9), &wits));
268
269        let n = ExactInt::from_u64(8051);
270        let f = pollard_rho(&n).expect("rho");
271        let other = n.div_rem(&f).expect("div").0;
272        let a = ExactInt::from_u64(83);
273        let b = ExactInt::from_u64(97);
274        assert!(
275            (f == a && other == b) || (f == b && other == a),
276            "factor {f:?}"
277        );
278        assert!(pollard_rho(&ExactInt::from_u64(31)).is_none());
279        assert!(mod_pow(&two, &ExactInt::from_i64(-1), &m).is_none());
280    }
281}
zenith_float_num/modular.rs

zenith_float_num/
modular.rs
