num_prime/
primality.rs

1use crate::traits::{BitTest, ExactRoots, PrimalityUtils};
2use either::Either;
3use num_integer::{Integer, Roots};
4use num_modular::{ModularCoreOps, ModularRefOps, ModularUnaryOps};
5use num_traits::{FromPrimitive, NumRef, RefNum, ToPrimitive};
6
7/// Utilities for the Lucas pseudoprime test
8pub trait LucasUtils {
9    /// Find Lucas sequence to n with modulo m, i.e. $U_{n+1}(P,Q)$ mod m and $V_{n+1}(P,Q)$
10    fn lucasm(p: usize, q: isize, m: Self, n: Self) -> (Self, Self)
11    where
12        Self: Sized;
13
14    /// Find proper parameters P and Q of Lucas pseudoprime test for n, such that Jacobi(D|n) is -1,
15    /// using Selfridge's method (referred as method A by Baillie).
16    fn pq_selfridge(n: &Self) -> (usize, isize);
17
18    /// Find proper parameters P of extra strong Lucas pseudoprime test for n, such that Jacobi(D|n) is -1,
19    /// using brute-force searching (referred as method C by Baillie)
20    fn p_bruteforce(n: &Self) -> usize;
21}
22
23impl<T> LucasUtils for T
24where
25    T: Integer
26        + FromPrimitive
27        + ToPrimitive
28        + NumRef
29        + BitTest
30        + ExactRoots
31        + Clone
32        + ModularRefOps,
33    for<'r> &'r T:
34        RefNum<T> + ModularUnaryOps<&'r T, Output = T> + ModularCoreOps<&'r T, &'r T, Output = T>,
35{
36    fn lucasm(p: usize, q: isize, m: Self, n: Self) -> (Self, Self) {
37        // Reference: <https://en.wikipedia.org/wiki/Lucas_sequence>
38        // and Peter Hackman, "Elementary Number Theory", section "L.XVII Scalar" <http://hackmat.se/kurser/TATM54/booktot.pdf>
39
40        let p = T::from_usize(p).unwrap() % &m;
41        let q = if q >= 0 {
42            T::from_isize(q).unwrap() % &m
43        } else {
44            T::from_isize(-q).unwrap().negm(&m)
45        };
46
47        let mut uk = T::zero() % &m; // U(k), mod m for montgomery form
48        let mut uk1 = T::one() % &m; // U(k+1)
49
50        for i in (0..n.bits()).rev() {
51            if n.bit(i) {
52                // k' = 2k+1
53                // U(k'+1) = U(2k+2) = PU(k+1)² - 2*QU(k+1)U(k)
54                let uk1sq = (&uk1).sqm(&m);
55                let t1 = (&p).mulm(&uk1sq, &m);
56                let t2 = (&q).mulm(&uk1, &m).mulm(&uk, &m).dblm(&m);
57                let new_uk1 = t1.subm(&t2, &m);
58                // U(k') = U(2k+1) = U(k+1)² - QU(k)²
59                let t1 = uk1sq;
60                let t2 = (&q).mulm(&uk.sqm(&m), &m);
61                let new_uk = t1.subm(&t2, &m);
62                uk1 = new_uk1;
63                uk = new_uk;
64            } else {
65                // k' = 2k
66                // U(k'+1) = U(2k+1) = U(k+1)² - QU(k)²
67                let uksq = (&uk).sqm(&m);
68                let t1 = (&uk1).sqm(&m);
69                let t2 = (&uksq).mulm(&q, &m);
70                let new_uk1 = t1.subm(&t2, &m);
71                // U(k') = U(2k) = 2U(k+1)U(k) - PU(k)²
72                let t1 = uk1.mulm(&uk, &m).dblm(&m);
73                let t2 = uksq.mulm(&p, &m);
74                let new_uk = t1.subm(&t2, &m);
75                uk1 = new_uk1;
76                uk = new_uk;
77            }
78        }
79
80        // V(k) = 2U(k+1) - PU(k)
81        let vk = (&uk1).dblm(&m).subm(&p.mulm(&uk, &m), &m);
82        (uk, vk)
83    }
84
85    fn pq_selfridge(n: &Self) -> (usize, isize) {
86        let mut d = T::from_u8(5).unwrap();
87        let mut neg = false;
88        loop {
89            // check if n is a square number after several trials
90            if d == T::from_u8(13).unwrap() && (*n).is_square() {
91                break (0, 0);
92            }
93
94            let sd = if neg { (&d).negm(n) } else { d.clone() };
95            let j = sd.jacobi(n);
96            if j == 0 && &d != n {
97                break (0, 0);
98            } // modification from Baillie, see https://oeis.org/A217120/a217120_1.txt
99            if j == -1 {
100                let d = if neg {
101                    -d.to_isize().unwrap()
102                } else {
103                    d.to_isize().unwrap()
104                };
105                break (1, (1 - d) / 4);
106            }
107
108            d = d + T::from_u8(2).unwrap();
109            neg = !neg;
110        }
111    }
112
113    fn p_bruteforce(n: &Self) -> usize {
114        let mut p: usize = 3;
115        loop {
116            // check if n is a square number after several trials
117            if p == 10 && (*n).is_square() {
118                break 0;
119            }
120
121            let d = T::from_usize(p * p - 4).unwrap();
122            let j = d.jacobi(n);
123            if j == 0 && &d != n {
124                break 0;
125            }
126            if j == -1 {
127                break p;
128            }
129            p += 1;
130        }
131    }
132}
133
134impl<T: Integer + NumRef + Clone + FromPrimitive + LucasUtils + BitTest + ModularRefOps>
135    PrimalityUtils for T
136where
137    for<'r> &'r T:
138        RefNum<T> + std::ops::Shr<usize, Output = T> + ModularUnaryOps<&'r T, Output = T>,
139{
140    #[inline]
141    fn is_prp(&self, base: Self) -> bool {
142        if self < &Self::one() {
143            return false;
144        }
145        let tm1 = self - Self::one();
146        base.powm(&tm1, self).is_one()
147    }
148
149    #[inline]
150    fn is_sprp(&self, base: Self) -> bool {
151        self.test_sprp(base).either(|v| v, |_| false)
152    }
153
154    fn test_sprp(&self, base: T) -> Either<bool, T> {
155        if self < &Self::one() {
156            return Either::Left(false);
157        }
158
159        // find 2^shift*u + 1 = n
160        let tm1 = self - T::one();
161        let shift = tm1.trailing_zeros();
162        let u = &tm1 >> shift;
163
164        // prevent reduction if the input is in montgomery form
165        let m1 = T::one() % self;
166        let mm1 = (&m1).negm(self);
167
168        let mut x = base.powm(&u, self);
169        if x == m1 || x == mm1 {
170            return Either::Left(true);
171        }
172
173        for _ in 0..shift {
174            let y = (&x).sqm(self);
175            if y == m1 {
176                return Either::Right(self.gcd(&(x - T::one())));
177            }
178            if y == mm1 {
179                return Either::Left(true);
180            }
181            x = y;
182        }
183
184        Either::Left(x == m1)
185    }
186
187    fn is_lprp(&self, p: Option<usize>, q: Option<isize>) -> bool {
188        if self < &Self::one() {
189            return false;
190        }
191        if self.is_even() {
192            return false;
193        }
194
195        let (p, q) = match (p, q) {
196            (Some(sp), Some(sq)) => (sp, sq),
197            (_, _) => {
198                let (sp, sq) = LucasUtils::pq_selfridge(self);
199                if sp == 0 {
200                    return false;
201                }; // is a perfect power
202                (sp, sq)
203            }
204        };
205
206        let d = (p * p) as isize - 4 * q;
207        let d = if d > 0 {
208            Self::from_isize(d).unwrap()
209        } else {
210            Self::from_isize(-d).unwrap().negm(self)
211        };
212        let delta = match d.jacobi(self) {
213            0 => self.clone(),
214            -1 => self + Self::one(),
215            1 => self - Self::one(),
216            _ => unreachable!(),
217        };
218        let (u, _) = LucasUtils::lucasm(p, q, self.clone(), delta);
219        u.is_zero()
220    }
221
222    fn is_slprp(&self, p: Option<usize>, q: Option<isize>) -> bool {
223        if self < &Self::one() {
224            return false;
225        }
226        if self.is_even() {
227            return false;
228        }
229
230        let (p, q) = match (p, q) {
231            (Some(sp), Some(sq)) => (sp, sq),
232            (_, _) => {
233                let (sp, sq) = LucasUtils::pq_selfridge(self);
234                if sp == 0 {
235                    return false;
236                }; // is a perfect power
237                (sp, sq)
238            }
239        };
240
241        let d = (p * p) as isize - 4 * q;
242        let d = if d > 0 {
243            Self::from_isize(d).unwrap()
244        } else {
245            Self::from_isize(-d).unwrap().negm(self)
246        };
247        let delta = match d.jacobi(self) {
248            0 => self.clone(),
249            -1 => self + Self::one(),
250            1 => self - Self::one(),
251            _ => unreachable!(),
252        };
253
254        let shift = delta.trailing_zeros();
255        let base = &delta >> shift;
256
257        let (ud, vd) = LucasUtils::lucasm(p, q, self.clone(), base.clone());
258        if ud.is_zero() || vd.is_zero() {
259            return true;
260        }
261
262        // do first iteration to remove the sign on Q
263        if shift == 0 {
264            return false;
265        }
266        let mut qk = Self::from_isize(q.abs()).unwrap().powm(&base, self);
267        // V(2k) = V(k)² - 2Q^k
268        let mut vd = if q >= 0 {
269            vd.sqm(self).subm(&(&qk).dblm(self), self)
270        } else {
271            vd.sqm(self).addm(&(&qk).dblm(self), self)
272        };
273        if vd.is_zero() {
274            return true;
275        }
276
277        for _ in 1..shift {
278            // V(2k) = V(k)² - 2Q^k
279            qk = qk.sqm(self);
280            vd = vd.sqm(self).subm(&(&qk).dblm(self), self);
281            if vd.is_zero() {
282                return true;
283            }
284        }
285        false
286    }
287
288    fn is_eslprp(&self, p: Option<usize>) -> bool {
289        if self < &Self::one() {
290            return false;
291        }
292        if self.is_even() {
293            return false;
294        }
295
296        let p = match p {
297            Some(sp) => sp,
298            None => {
299                let sp = LucasUtils::p_bruteforce(self);
300                if sp == 0 {
301                    return false;
302                }; // is a perfect power
303                sp
304            }
305        };
306
307        let d = (p * p) as isize - 4;
308        let d = if d > 0 {
309            Self::from_isize(d).unwrap()
310        } else {
311            Self::from_isize(-d).unwrap().negm(self)
312        };
313        let delta = match d.jacobi(self) {
314            0 => self.clone(),
315            -1 => self + Self::one(),
316            1 => self - Self::one(),
317            _ => unreachable!(),
318        };
319
320        let shift = delta.trailing_zeros();
321        let base = &delta >> shift;
322
323        let (ud, mut vd) = LucasUtils::lucasm(p, 1, self.clone(), base.clone());
324        let two = Self::from_u8(2).unwrap();
325        // U(d) = 0 or V(d) = ±2
326        if ud.is_zero() && (vd == two || vd == self - &two) {
327            return true;
328        }
329        if vd.is_zero() {
330            return true;
331        }
332
333        for _ in 1..(shift - 1) {
334            // V(2k) = V(k)² - 2
335            vd = vd.sqm(self).subm(&two, self);
336            if vd.is_zero() {
337                return true;
338            }
339        }
340        false
341    }
342}
343
344/// A dummy trait for integer type. All types that implements this and [`PrimalityRefBase`]
345/// will be supported by most functions in `num-primes`
346pub trait PrimalityBase:
347    Integer
348    + Roots
349    + NumRef
350    + Clone
351    + FromPrimitive
352    + ToPrimitive
353    + ExactRoots
354    + BitTest
355    + ModularRefOps
356{
357}
358impl<
359        T: Integer
360            + Roots
361            + NumRef
362            + Clone
363            + FromPrimitive
364            + ToPrimitive
365            + ExactRoots
366            + BitTest
367            + ModularRefOps,
368    > PrimalityBase for T
369{
370}
371
372/// A dummy trait for integer reference type. All types that implements this and [`PrimalityBase`]
373/// will be supported by most functions in `num-primes`
374pub trait PrimalityRefBase<Base>:
375    RefNum<Base>
376    + std::ops::Shr<usize, Output = Base>
377    + for<'r> ModularUnaryOps<&'r Base, Output = Base>
378    + for<'r> ModularCoreOps<&'r Base, &'r Base, Output = Base>
379{
380}
381impl<T, Base> PrimalityRefBase<Base> for T where
382    T: RefNum<Base>
383        + std::ops::Shr<usize, Output = Base>
384        + for<'r> ModularUnaryOps<&'r Base, Output = Base>
385        + for<'r> ModularCoreOps<&'r Base, &'r Base, Output = Base>
386{
387}
388
389#[cfg(test)]
390mod tests {
391    use super::*;
392    use crate::mint::SmallMint;
393    use num_modular::{ModularAbs, ModularSymbols};
394    use rand::random;
395
396    #[cfg(feature = "num-bigint")]
397    use num_bigint::BigUint;
398
399    #[test]
400    fn fermat_prp_test() {
401        // 341 is the smallest pseudoprime for base 2
402        assert!(341u16.is_prp(2));
403        assert!(!340u16.is_prp(2));
404        assert!(!105u16.is_prp(2));
405
406        // Carmichael number 561 = 3*11*17 is fermat pseudoprime for any base coprime to 561
407        for p in [2, 5, 7, 13, 19] {
408            assert!(561u32.is_prp(p));
409        }
410    }
411
412    #[test]
413    fn sprp_test() {
414        // strong pseudoprimes of base 2 (OEIS:A001262) under 10000
415        let spsp: [u16; 5] = [2047, 3277, 4033, 4681, 8321];
416        for psp in spsp {
417            assert!(psp.is_sprp(2));
418            assert!(SmallMint::from(psp).is_sprp(2.into())); // test Mint execution
419        }
420
421        // test cofactor return
422        assert_eq!(341u16.test_sprp(2), Either::Right(31));
423        assert_eq!(
424            SmallMint::from(341u16).test_sprp(2.into()),
425            Either::Right(31.into())
426        );
427    }
428
429    #[test]
430    fn lucas_mod_test() {
431        // OEIS:A006190
432        let p3qm1: [u64; 26] = [
433            0,
434            1,
435            3,
436            10,
437            33,
438            109,
439            360,
440            1189,
441            3927,
442            12970,
443            42837,
444            141_481,
445            467_280,
446            1_543_321,
447            5_097_243,
448            16_835_050,
449            55_602_393,
450            183_642_229,
451            606_529_080,
452            2_003_229_469,
453            6_616_217_487,
454            21_851_881_930,
455            72_171_863_277,
456            238_367_471_761,
457            787_274_278_560,
458            2_600_190_307_441,
459        ];
460        let m = random::<u16>();
461        for (n, &p3qm1_val) in p3qm1.iter().enumerate().skip(2) {
462            let (uk, _) = LucasUtils::lucasm(3, -1, u64::from(m), n as u64);
463            assert_eq!(uk, p3qm1_val % u64::from(m));
464
465            #[cfg(feature = "num-bigint")]
466            {
467                let (uk, _) = LucasUtils::lucasm(3, -1, BigUint::from(m), BigUint::from(n));
468                assert_eq!(uk, BigUint::from(p3qm1_val % u64::from(m)));
469            }
470        }
471
472        fn lucasm_naive(p: usize, q: isize, m: u16, n: u16) -> (u16, u16) {
473            if n == 0 {
474                return (0, 2);
475            }
476
477            let m = m as usize;
478            let q = q.absm(&m);
479            let (mut um1, mut u) = (0, 1); // U_{n-1}, U_{n}
480            let (mut vm1, mut v) = (2, p % m); // V_{n-1}, V_{n}
481
482            for _ in 1..n {
483                let new_u = p.mulm(u, &m).subm(q.mulm(um1, &m), &m);
484                um1 = u;
485                u = new_u;
486
487                let new_v = p.mulm(v, &m).subm(q.mulm(vm1, &m), &m);
488                vm1 = v;
489                v = new_v;
490            }
491            (u as u16, v as u16)
492        }
493        for _ in 0..10 {
494            let n = u16::from(random::<u8>());
495            let m = random::<u16>();
496            let p = random::<u16>() as usize;
497            let q = random::<i16>() as isize;
498            assert_eq!(
499                LucasUtils::lucasm(p, q, m, n),
500                lucasm_naive(p, q, m, n),
501                "failed with Lucas settings: p={p}, q={q}, m={m}, n={n}"
502            );
503        }
504    }
505
506    #[test]
507    fn lucas_prp_test() {
508        // Some known cases
509        assert!(19u8.is_lprp(Some(3), Some(-1)));
510        assert!(5u8.is_lprp(None, None));
511        assert!(7u8.is_lprp(None, None));
512        assert!(!9u8.is_lprp(None, None));
513        assert!(!5719u16.is_lprp(None, None));
514        assert!(!1239u16.is_eslprp(None));
515
516        // least lucas pseudo primes for Q=-1 and Jacobi(D/n) = -1 (from Wikipedia)
517        let plimit: [u16; 5] = [323, 35, 119, 9, 9];
518        for (i, l) in plimit.iter().copied().enumerate() {
519            let p = i + 1;
520            assert!(l.is_lprp(Some(p), Some(-1)));
521
522            // test four random numbers under limit
523            for _ in 0..10 {
524                let n = random::<u16>() % l;
525                if n <= 3 || n.is_sprp(2) {
526                    continue;
527                } // skip real primes
528                let d = (p * p + 4) as u16;
529                if n.is_odd() && d.jacobi(&n) != -1 {
530                    continue;
531                }
532                assert!(
533                    !n.is_lprp(Some(p), Some(-1)),
534                    "lucas prp test on {} with p = {}",
535                    n,
536                    p
537                );
538            }
539        }
540
541        // least strong lucas pseudoprimes for Q=-1 and Jacobi(D/n) = -1 (from Wikipedia)
542        let plimit: [u16; 3] = [4181, 169, 119];
543        for (i, l) in plimit.iter().copied().enumerate() {
544            let p = i + 1;
545            assert!(l.is_slprp(Some(p), Some(-1)));
546
547            // test random numbers under limit
548            for _ in 0..10 {
549                let n = random::<u16>() % l;
550                if n <= 3 || (n.is_sprp(2) && n.is_sprp(3)) {
551                    continue;
552                } // skip real primes
553                let d = (p * p + 4) as u16;
554                if n.is_odd() && d.jacobi(&n) != -1 {
555                    continue;
556                }
557                assert!(
558                    !n.is_slprp(Some(p), Some(-1)),
559                    "strong lucas prp test on {} with p = {}",
560                    n,
561                    p
562                );
563            }
564        }
565
566        // lucas pseudoprimes (OEIS:A217120) under 10000
567        let lpsp: [u16; 9] = [323, 377, 1159, 1829, 3827, 5459, 5777, 9071, 9179];
568        for psp in lpsp {
569            assert!(
570                psp.is_lprp(None, None),
571                "lucas prp test on pseudoprime {}",
572                psp
573            );
574        }
575        for _ in 0..50 {
576            let n = random::<u16>() % 10000;
577            if n <= 3 || (n.is_sprp(2) && n.is_sprp(3)) {
578                continue;
579            } // skip real primes
580            if lpsp.iter().any(|x| x == &n) {
581                continue;
582            } // skip pseudoprimes
583            assert!(!n.is_lprp(None, None), "lucas prp test on {}", n);
584        }
585
586        // strong lucas pseudoprimes (OEIS:A217255) under 10000
587        let slpsp: [u16; 2] = [5459, 5777];
588        for psp in slpsp {
589            assert!(
590                psp.is_slprp(None, None),
591                "strong lucas prp test on pseudoprime {}",
592                psp
593            );
594        }
595        for _ in 0..50 {
596            let n = random::<u16>() % 10000;
597            if n <= 3 || (n.is_sprp(2) && n.is_sprp(3)) {
598                continue;
599            } // skip real primes
600            if slpsp.iter().any(|x| x == &n) {
601                continue;
602            } // skip pseudoprimes
603            assert!(!n.is_slprp(None, None), "strong lucas prp test on {}", n);
604        }
605
606        // extra strong lucas pseudoprimes (OEIS:A217719) under 10000
607        let eslpsp: [u16; 3] = [989, 3239, 5777];
608        for psp in eslpsp {
609            assert!(
610                psp.is_eslprp(None),
611                "extra strong lucas prp test on pseudoprime {}",
612                psp
613            );
614        }
615        for _ in 0..50 {
616            let n = random::<u16>() % 10000;
617            if n <= 3 || (n.is_sprp(2) && n.is_sprp(3)) {
618                continue;
619            } // skip real primes
620            if eslpsp.iter().any(|x| x == &n) {
621                continue;
622            } // skip pseudoprimes
623            assert!(!n.is_eslprp(None), "extra strong lucas prp test on {}", n);
624        }
625    }
626}
num_prime/primality.rs

num_prime/
primality.rs
