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            141481,
445            467280,
446            1543321,
447            5097243,
448            16835050,
449            55602393,
450            183642229,
451            606529080,
452            2003229469,
453            6616217487,
454            21851881930,
455            72171863277,
456            238367471761,
457            787274278560,
458            2600190307441,
459        ];
460        let m = random::<u16>();
461        for n in 2..p3qm1.len() {
462            let (uk, _) = LucasUtils::lucasm(3, -1, m as u64, n as u64);
463            assert_eq!(uk, p3qm1[n] % (m as u64));
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[n] % (m as u64)));
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 = random::<u8>() as u16;
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={}, q={}, m={}, n={}",
502                p,
503                q,
504                m,
505                n
506            );
507        }
508    }
509
510    #[test]
511    fn lucas_prp_test() {
512        // Some known cases
513        assert!(19u8.is_lprp(Some(3), Some(-1)));
514        assert!(5u8.is_lprp(None, None));
515        assert!(7u8.is_lprp(None, None));
516        assert!(!9u8.is_lprp(None, None));
517        assert!(!5719u16.is_lprp(None, None));
518        assert!(!1239u16.is_eslprp(None));
519
520        // least lucas pseudo primes for Q=-1 and Jacobi(D/n) = -1 (from Wikipedia)
521        let plimit: [u16; 5] = [323, 35, 119, 9, 9];
522        for (i, l) in plimit.iter().cloned().enumerate() {
523            let p = i + 1;
524            assert!(l.is_lprp(Some(p), Some(-1)));
525
526            // test four random numbers under limit
527            for _ in 0..10 {
528                let n = random::<u16>() % l;
529                if n <= 3 || n.is_sprp(2) {
530                    continue;
531                } // skip real primes
532                let d = (p * p + 4) as u16;
533                if n.is_odd() && d.jacobi(&n) != -1 {
534                    continue;
535                }
536                assert!(
537                    !n.is_lprp(Some(p), Some(-1)),
538                    "lucas prp test on {} with p = {}",
539                    n,
540                    p
541                );
542            }
543        }
544
545        // least strong lucas pseudoprimes for Q=-1 and Jacobi(D/n) = -1 (from Wikipedia)
546        let plimit: [u16; 3] = [4181, 169, 119];
547        for (i, l) in plimit.iter().cloned().enumerate() {
548            let p = i + 1;
549            assert!(l.is_slprp(Some(p), Some(-1)));
550
551            // test random numbers under limit
552            for _ in 0..10 {
553                let n = random::<u16>() % l;
554                if n <= 3 || (n.is_sprp(2) && n.is_sprp(3)) {
555                    continue;
556                } // skip real primes
557                let d = (p * p + 4) as u16;
558                if n.is_odd() && d.jacobi(&n) != -1 {
559                    continue;
560                }
561                assert!(
562                    !n.is_slprp(Some(p), Some(-1)),
563                    "strong lucas prp test on {} with p = {}",
564                    n,
565                    p
566                );
567            }
568        }
569
570        // lucas pseudoprimes (OEIS:A217120) under 10000
571        let lpsp: [u16; 9] = [323, 377, 1159, 1829, 3827, 5459, 5777, 9071, 9179];
572        for psp in lpsp {
573            assert!(
574                psp.is_lprp(None, None),
575                "lucas prp test on pseudoprime {}",
576                psp
577            );
578        }
579        for _ in 0..50 {
580            let n = random::<u16>() % 10000;
581            if n <= 3 || (n.is_sprp(2) && n.is_sprp(3)) {
582                continue;
583            } // skip real primes
584            if lpsp.iter().find(|&x| x == &n).is_some() {
585                continue;
586            } // skip pseudoprimes
587            assert!(!n.is_lprp(None, None), "lucas prp test on {}", n);
588        }
589
590        // strong lucas pseudoprimes (OEIS:A217255) under 10000
591        let slpsp: [u16; 2] = [5459, 5777];
592        for psp in slpsp {
593            assert!(
594                psp.is_slprp(None, None),
595                "strong lucas prp test on pseudoprime {}",
596                psp
597            );
598        }
599        for _ in 0..50 {
600            let n = random::<u16>() % 10000;
601            if n <= 3 || (n.is_sprp(2) && n.is_sprp(3)) {
602                continue;
603            } // skip real primes
604            if slpsp.iter().find(|&x| x == &n).is_some() {
605                continue;
606            } // skip pseudoprimes
607            assert!(!n.is_slprp(None, None), "strong lucas prp test on {}", n);
608        }
609
610        // extra strong lucas pseudoprimes (OEIS:A217719) under 10000
611        let eslpsp: [u16; 3] = [989, 3239, 5777];
612        for psp in eslpsp {
613            assert!(
614                psp.is_eslprp(None),
615                "extra strong lucas prp test on pseudoprime {}",
616                psp
617            );
618        }
619        for _ in 0..50 {
620            let n = random::<u16>() % 10000;
621            if n <= 3 || (n.is_sprp(2) && n.is_sprp(3)) {
622                continue;
623            } // skip real primes
624            if eslpsp.iter().find(|&x| x == &n).is_some() {
625                continue;
626            } // skip pseudoprimes
627            assert!(!n.is_eslprp(None), "extra strong lucas prp test on {}", n);
628        }
629    }
630}
num_prime/primality.rs

num_prime/
primality.rs
