crypto_primes/hazmat/
lucas.rs

1//! Lucas primality test.
2use core::num::NonZero;
3use crypto_bigint::{Integer, Limb, Monty, Odd, Square, Word};
4
5use super::{
6    gcd::gcd_vartime,
7    jacobi::{jacobi_symbol_vartime, JacobiSymbol},
8    Primality,
9};
10
11/// The maximum number of attempts to find `D` such that `(D/n) == -1`.
12// This is widely believed to be impossible.
13// So if we exceed it, we will panic reporting the value of `n`.
14const MAX_ATTEMPTS: usize = 10_000;
15
16/// The number of attempts to find `D` such that `(D/n) == -1`
17/// before checking that `n` is a square (in which case such `D` does not exist).
18// This check is relatively expensive compared to calculating the Jacobi symbol
19// (~30x for 1024-bit numbers, ~100x for 2048 bit).
20// On the other hand, if `n` is a non-square we expect to find a `D`
21// in just a few attempts on average (an estimate for the Selfridge method
22// can be found in [^Baillie1980], section 7; for the brute force method
23// it seems to be about the same).
24const ATTEMPTS_BEFORE_SQRT: usize = 30;
25
26/// A method for selecting the base `(P, Q)` for the Lucas primality test.
27pub trait LucasBase {
28    /// Given an odd integer, returns `Ok((P, abs(Q), is_negative(Q)))` on success,
29    /// or `Err(Primality)` if the primality for the given integer was discovered
30    /// during the search for a base.
31    fn generate<T: Integer>(&self, n: &Odd<T>) -> Result<(Word, Word, bool), Primality>;
32}
33
34/// "Method A" for selecting the base given in Baillie & Wagstaff[^Baillie1980],
35/// attributed to Selfridge.
36///
37/// Try `D = 1 - 4Q = 5, -7, 9, -11, 13, ...` until `Jacobi(D, n) = -1`.
38/// Return `P = 1, Q = (1 - D) / 4)`.
39///
40/// [^Baillie1980]:
41///   R. Baillie, S. S. Wagstaff, "Lucas pseudoprimes",
42///   Math. Comp. 35 1391-1417 (1980),
43///   DOI: [10.2307/2006406](https://dx.doi.org/10.2307/2006406),
44///   <http://mpqs.free.fr/LucasPseudoprimes.pdf>
45#[derive(Copy, Clone, Debug, PartialEq, Eq)]
46pub struct SelfridgeBase;
47
48impl LucasBase for SelfridgeBase {
49    fn generate<T: Integer>(&self, n: &Odd<T>) -> Result<(Word, Word, bool), Primality> {
50        let mut abs_d = 5;
51        let mut d_is_negative = false;
52        let n_is_small = n.bits_vartime() < Word::BITS; // if true, `n` fits into one `Word`
53        let small_n = n.as_ref().as_ref()[0].0;
54        let mut attempts = 0;
55        loop {
56            if attempts >= MAX_ATTEMPTS {
57                panic!("internal error: cannot find (D/n) = -1 for {:?}", n)
58            }
59
60            if attempts >= ATTEMPTS_BEFORE_SQRT {
61                let sqrt_n = n.sqrt_vartime();
62                if &sqrt_n.wrapping_mul(&sqrt_n) == n.as_ref() {
63                    return Err(Primality::Composite);
64                }
65            }
66
67            let j = jacobi_symbol_vartime(abs_d, d_is_negative, n);
68
69            if j == JacobiSymbol::MinusOne {
70                break;
71            }
72            if j == JacobiSymbol::Zero {
73                // Modification of Method A by Baillie, in an example to OEIS:A217120
74                // (https://oeis.org/A217120/a217120_1.txt):
75                // If `d == (+,-)n`, (e.g., `n` = 5 or 11) try the next `d` instead of quitting;
76                // this small modification of Selfridge's method A
77                // enables 5 and 11 to be classified as Lucas probable primes.
78                // Otherwise GCD(D, n) > 1, and therefore n is not prime.
79                if !(n_is_small && small_n == abs_d) {
80                    return Err(Primality::Composite);
81                }
82            }
83
84            attempts += 1;
85            d_is_negative = !d_is_negative;
86            abs_d += 2;
87        }
88
89        // Calculate `q = (1 - d) / 4`.
90        // No remainder from division by 4, by construction of `d`.
91        let (abs_q, q_is_negative) = if d_is_negative {
92            ((abs_d + 1) / 4, false)
93        } else {
94            ((abs_d - 1) / 4, true)
95        };
96
97        Ok((1, abs_q, q_is_negative))
98    }
99}
100
101/// "Method A*" for selecting the base given in Baillie & Wagstaff[^Baillie1980].
102///
103/// Same as [`SelfridgeBase`], but returns `(P = 5, Q = 5)` if the Selfridge base set `Q = -1`.
104///
105/// [^Baillie1980]:
106///   R. Baillie, S. S. Wagstaff, "Lucas pseudoprimes",
107///   Math. Comp. 35 1391-1417 (1980),
108///   DOI: [10.2307/2006406](https://dx.doi.org/10.2307/2006406),
109///   <http://mpqs.free.fr/LucasPseudoprimes.pdf>
110#[derive(Copy, Clone, Debug, PartialEq, Eq)]
111pub struct AStarBase;
112
113impl LucasBase for AStarBase {
114    fn generate<T: Integer>(&self, n: &Odd<T>) -> Result<(Word, Word, bool), Primality> {
115        SelfridgeBase.generate(n).map(|(p, abs_q, q_is_negative)| {
116            if abs_q == 1 && q_is_negative {
117                (5, 5, false)
118            } else {
119                (p, abs_q, q_is_negative)
120            }
121        })
122    }
123}
124
125/// "Method C" for selecting the base given by Baillie[^Baillie].
126///
127/// Try `P = 3, 4, 5, ...` until `Jacobi(D, n) = -1`, where `D = P^2 - 4Q`.
128/// Returns the found `P`, and `Q = 1`.
129///
130/// [^Baillie]:
131///   R. Baillie, Mathematica code for extra strong Lucas pseudoprimes,
132///   <https://oeis.org/A217719/a217719.txt>
133#[derive(Copy, Clone, Debug, PartialEq, Eq)]
134pub struct BruteForceBase;
135
136impl LucasBase for BruteForceBase {
137    fn generate<T: Integer>(&self, n: &Odd<T>) -> Result<(Word, Word, bool), Primality> {
138        let mut p = 3;
139        let mut attempts = 0;
140
141        loop {
142            if attempts >= MAX_ATTEMPTS {
143                panic!("internal error: cannot find (D/n) = -1 for {:?}", n)
144            }
145
146            if attempts >= ATTEMPTS_BEFORE_SQRT {
147                let sqrt_n = n.sqrt_vartime();
148                if &sqrt_n.wrapping_mul(&sqrt_n) == n.as_ref() {
149                    return Err(Primality::Composite);
150                }
151            }
152
153            // Can unwrap here since `p` is always small (see the condition above).
154            let j = jacobi_symbol_vartime(p * p - 4, false, n);
155
156            if j == JacobiSymbol::MinusOne {
157                break;
158            }
159            if j == JacobiSymbol::Zero {
160                // D = P^2 - 4 = (P - 2)(P + 2).
161                // If (D/n) == 0 then D shares a prime factor with n.
162                // Since the loop proceeds in increasing P and starts with P - 2 == 1,
163                // the shared prime factor must be P + 2.
164                // If P + 2 == n, then n is prime; otherwise P + 2 is a proper factor of n.
165                let primality = if n.as_ref() == &T::from_limb_like(Limb::from(p + 2), n.as_ref()) {
166                    Primality::Prime
167                } else {
168                    Primality::Composite
169                };
170                return Err(primality);
171            }
172
173            attempts += 1;
174            p += 1;
175        }
176
177        Ok((p, 1, false))
178    }
179}
180
181/// For the given odd `n`, finds `s` and odd `d` such that `n + 1 == 2^s * d`.
182fn decompose<T: Integer>(n: &Odd<T>) -> (u32, Odd<T>) {
183    // Need to be careful here since `n + 1` can overflow.
184    // Instead of adding 1 and counting trailing 0s, we count trailing ones on the original `n`.
185
186    let one = T::one_like(n);
187    let s = n.trailing_ones_vartime();
188    let d = if s < n.bits_precision() {
189        // The shift won't overflow because of the check above.
190        // The addition won't overflow since the original `n` was odd,
191        // so we right-shifted at least once.
192        n.as_ref()
193            .overflowing_shr_vartime(s)
194            .expect("shift should be within range by construction")
195            .checked_add(&one)
196            .expect("addition should not overflow by construction")
197    } else {
198        one
199    };
200
201    (s, Odd::new(d).expect("`d` should be odd by construction"))
202}
203
204/// The checks to perform in the Lucas test.
205///
206/// Given the Lucas sequence built from some base `(P, Q)` (see [`LucasBase`])
207/// up to the elements `V(d)`, `U(d)`, where `d * 2^s == n - (D/n)`, `d` odd, and `D = P^2 - 4Q`,
208/// the checks are defined as follows:
209#[derive(Copy, Clone, Debug, PartialEq, Eq)]
210pub enum LucasCheck {
211    /// Introduced by Baillie & Wagstaff[^Baillie1980].
212    /// If either of the following is true:
213    /// - any of `V(d*2^r) == 0` for `0 <= r < s`,
214    /// - `U(d) == 0`,
215    ///
216    /// report the number as prime.
217    ///
218    /// If the base is [`SelfridgeBase`], known false positives constitute OEIS:A217255[^A217255].
219    ///
220    /// [^Baillie1980]:
221    ///   R. Baillie, S. S. Wagstaff, "Lucas pseudoprimes",
222    ///   Math. Comp. 35 1391-1417 (1980),
223    ///   DOI: [10.2307/2006406](https://dx.doi.org/10.2307/2006406),
224    ///   <http://mpqs.free.fr/LucasPseudoprimes.pdf>
225    ///
226    /// [^A217255]: <https://oeis.org/A217255>
227    Strong,
228
229    /// A [`LucasCheck::ExtraStrong`] without checking for `U(d)`.
230    /// That is, if either of the following is true:
231    /// - any of `V(d*2^r) == 0` for `0 <= r < s`,
232    /// - `V(d) == ±2`,
233    ///
234    /// report the number as prime.
235    ///
236    /// Note: the second condition is only checked if `Q == 1`,
237    /// otherwise it is considered to be true.
238    ///
239    /// If the base is [`BruteForceBase`], some known false positives
240    /// are listed by Jacobsen[^Jacobsen].
241    ///
242    /// Note: this option is intended for testing against known pseudoprimes;
243    /// do not use unless you know what you are doing.
244    ///
245    /// [^Jacobsen]:
246    ///   D. Jacobsen, "Pseudoprime Statistics, Tables, and Data",
247    ///   <http://ntheory.org/pseudoprimes.html>
248    AlmostExtraStrong,
249
250    /// Introduced by Mo[^Mo1993], and also described by Grantham[^Grantham2001].
251    /// If either of the following is true:
252    /// - any of `V(d*2^r) == 0` for `0 <= r < s`,
253    /// - `U(d) == 0` and `V(d) == ±2`,
254    ///
255    /// report the number as prime.
256    ///
257    /// Note that this check only differs from [`LucasCheck::Strong`] if `Q == 1`.
258    ///
259    /// If the base is [`BruteForceBase`], known false positives constitute OEIS:A217719[^A217719].
260    ///
261    /// [^Mo1993]:
262    ///   Zhaiyu Mo, "Diophantine equations, Lucas sequences and pseudoprimes",
263    ///   graduate thesis, University of Calgary, Calgary, AB (1993)
264    ///   DOI: [10.11575/PRISM/10820](https://dx.doi.org/10.11575/PRISM/10820)
265    ///
266    /// [^Grantham2001]:
267    ///   J. Grantham, "Frobenius pseudoprimes",
268    ///   Math. Comp. 70 873-891 (2001),
269    ///   DOI: [10.1090/S0025-5718-00-01197-2](https://dx.doi.org/10.1090/S0025-5718-00-01197-2)
270    ///
271    /// [^A217719]: <https://oeis.org/A217719>
272    ExtraStrong,
273
274    /// Introduced by Baillie et al[^Baillie2021].
275    /// If `V(n+1) == 2 Q`, report the number as prime.
276    ///
277    /// If the base is [`AStarBase`], some known false positives
278    /// are provided by Baillie et al[^Baillie2021].
279    ///
280    /// [^Baillie2021]: R. Baillie, A. Fiori, S. S. Wagstaff,
281    ///   "Strengthening the Baillie-PSW primality test",
282    ///   Math. Comp. 90 1931-1955 (2021),
283    ///   DOI: [10.1090/mcom/3616](https://doi.org/10.1090/mcom/3616)
284    LucasV,
285}
286
287/// Performs the primality test based on Lucas sequence.
288/// See [`LucasCheck`] for possible checks, and the implementors of [`LucasBase`]
289/// for the corresponding bases.
290pub fn lucas_test<T: Integer>(candidate: Odd<T>, base: impl LucasBase, check: LucasCheck) -> Primality {
291    // The comments in this function use references in `LucasCheck`, plus this one:
292    //
293    // [^Crandall2005]:
294    //   R. Crandall, C. Pomerance, "Prime numbers: a computational perspective",
295    //   2nd ed., Springer (2005) (ISBN: 0-387-25282-7, 978-0387-25282-7)
296
297    // A word-to-big integer conversion helper
298    let to_integer = |x: Word| T::from_limb_like(Limb::from(x), candidate.as_ref());
299
300    // Find the base for the Lucas sequence.
301    let (p, abs_q, q_is_negative) = match base.generate(&candidate) {
302        Ok(pq) => pq,
303        Err(primality) => return primality,
304    };
305
306    // Discriminant `d = p^2 - 4q`
307    let (abs_d, d_is_negative) = if q_is_negative {
308        (p * p + 4 * abs_q, false)
309    } else {
310        let t1 = p * p;
311        let t2 = 4 * abs_q;
312        if t2 > t1 {
313            (t2 - t1, true)
314        } else {
315            (t1 - t2, false)
316        }
317    };
318
319    // If either is true, it allows us to optimize certain parts of the calculations.
320    let p_is_one = p == 1;
321    let q_is_one = abs_q == 1 && !q_is_negative;
322
323    // See the references for the specific checks in the docstrings for [`LucasCheck`].
324
325    // All of the definitions require gcd(n, 2QD) == 1.
326    // We know gcd(n, D) = 1 by construction of the base (D is chosen such that (D/n) != 0).
327    // We know gcd(n, 2) = 1 (we checked for it earlier).
328    // In practice, gcd(n, Q) = 1 is always true, because the Lucas test is preceded by a sieve,
329    // and since `Q` is always small, division by it would have been already checked.
330    // But in order to avoid an implicit assumption that a sieve has been run,
331    // we check that gcd(n, Q) = 1 anyway - again, since `Q` is small,
332    // it does not noticeably affect the performance.
333    if abs_q != 1
334        && gcd_vartime(
335            candidate.as_ref(),
336            NonZero::new(abs_q).expect("q is not zero by construction"),
337        ) != 1
338        && candidate.as_ref() > &to_integer(abs_q)
339    {
340        return Primality::Composite;
341    }
342
343    // Find `d` and `s`, such that `d` is odd and `d * 2^s = n - (D/n)`.
344    // Since `(D/n) == -1` by construction, we're looking for `d * 2^s = n + 1`.
345    let (s, d) = decompose(&candidate);
346
347    // Some constants in Montgomery form
348    // TODO(dp): Cloning `candidate` should not be necessary. It's only needed because the
349    // `to_integer` closure captures it, but all we really need is the `BITS`
350    let params = <T as Integer>::Monty::new_params_vartime(candidate.clone());
351
352    let zero = <T as Integer>::Monty::zero(params.clone());
353    let one = <T as Integer>::Monty::one(params.clone());
354    let two = one.clone() + &one;
355    let minus_two = -two.clone();
356
357    // Convert Q to Montgomery form
358
359    let q = if q_is_one {
360        one.clone()
361    } else {
362        let abs_q = <T as Integer>::Monty::new(to_integer(abs_q), params.clone());
363        if q_is_negative {
364            -abs_q
365        } else {
366            abs_q
367        }
368    };
369
370    // Convert P to Montgomery form
371
372    let p = if p_is_one {
373        one.clone()
374    } else {
375        <T as Integer>::Monty::new(to_integer(p), params.clone())
376    };
377
378    // Compute d-th element of Lucas sequence (U_d(P, Q), V_d(P, Q)), where:
379    //
380    // V_0 = 2
381    // U_0 = 0
382    //
383    // U_{2k} = U_k V_k
384    // V_{2k} = V_k^2 - 2 Q^k
385    //
386    // U_{k+1} = (P U_k + V_k) / 2
387    // V_{k+1} = (D U_k + P V_k) / 2
388    //
389    // (The propagation method is due to [^Baillie2021], Eqs. 13, 14, 16, 17)
390    // We can therefore start with k=0 and build up to k=d in log2(d) steps.
391
392    // Starting with k = 0
393    let mut vk = two.clone(); // keeps V_k
394    let mut uk = <T as Integer>::Monty::zero(params.clone()); // keeps U_k
395    let mut qk = one.clone(); // keeps Q^k
396
397    // D in Montgomery representation - note that it can be negative.
398    let abs_d = <T as Integer>::Monty::new(to_integer(abs_d), params);
399    let d_m = if d_is_negative { -abs_d } else { abs_d };
400
401    for i in (0..d.bits_vartime()).rev() {
402        // k' = k * 2
403
404        let u_2k = uk * &vk;
405        let v_2k = vk.square() - &(qk.clone() + &qk);
406        let q_2k = qk.square();
407
408        uk = u_2k;
409        vk = v_2k;
410        qk = q_2k;
411
412        if d.bit_vartime(i) {
413            // k' = k + 1
414
415            let (p_uk, p_vk) = if p_is_one {
416                (uk.clone(), vk.clone())
417            } else {
418                (p.clone() * &uk, p.clone() * &vk)
419            };
420
421            let u_k1 = (p_uk + &vk).div_by_2();
422            let v_k1 = (d_m.clone() * &uk + &p_vk).div_by_2();
423            let q_k1 = qk * &q;
424
425            uk = u_k1;
426            vk = v_k1;
427            qk = q_k1;
428        }
429    }
430
431    // Now k=d, so vk = V_d and uk = U_d.
432
433    // Check for the first sufficient condition in various strong checks.
434
435    if check == LucasCheck::Strong && uk == zero {
436        // Strong check: `U_d == 0 mod n`.
437        return Primality::ProbablyPrime;
438    } else if check == LucasCheck::ExtraStrong || check == LucasCheck::AlmostExtraStrong {
439        // Extra strong check (from [^Mo1993]): `V_d == ±2 mod n` and `U_d == 0 mod n`.
440        //
441        // Note that the first identity only applies if `Q = 1`, since it is a consequence
442        // of a property of Lucas series: `V_k^2 - 4 Q^k = D U_k^2 mod n`.
443        // If `Q = 1` we can easily decompose the left side of the equation
444        // leading to the check above.
445        //
446        // If `Q != 1` we just consider it passed (we don't have a corresponding
447        // pseudoprime list anyway).
448
449        let vk_equals_two = !q_is_one || (vk == two || vk == minus_two);
450
451        if check == LucasCheck::ExtraStrong && uk == zero && vk_equals_two {
452            return Primality::ProbablyPrime;
453        }
454
455        // "Almost extra strong" check skips the `U_d` check.
456        // Since we have `U_d` anyway, it does not improve performance,
457        // so it is only here for testing purposes, since we have a corresponding pseudoprime list.
458        if check == LucasCheck::AlmostExtraStrong && vk_equals_two {
459            return Primality::ProbablyPrime;
460        }
461    }
462
463    // Second sufficient condition requires further propagating `V_k` up to `V_{n+1}`.
464
465    // Check if V_{2^t d} == 0 mod n for some 0 <= t < s.
466    // (unless we're in Lucas-V mode, then we just propagate V_k)
467
468    if check != LucasCheck::LucasV && vk == zero {
469        return Primality::ProbablyPrime;
470    }
471
472    for _ in 1..s {
473        // Optimization: V_k = ±2 is a fixed point for V_k' = V_k^2 - 2 Q^k with Q = 1,
474        // so if V_k = ±2, we can stop: we will never find a future V_k == 0.
475        if check != LucasCheck::LucasV && q_is_one && (vk == two || vk == minus_two) {
476            return Primality::Composite;
477        }
478
479        // k' = 2k
480        // V_{k'} = V_k^2 - 2 Q^k
481        vk = vk.square() - &qk - &qk;
482
483        if check != LucasCheck::LucasV && vk == zero {
484            return Primality::ProbablyPrime;
485        }
486
487        if !q_is_one {
488            qk = qk.square();
489        }
490    }
491
492    if check == LucasCheck::LucasV {
493        // At this point vk = V_{d * 2^(s-1)}.
494        // Double the index again:
495        vk = vk.square() - &qk - &qk; // now vk = V_{d * 2^s} = V_{n+1}
496
497        // Lucas-V check[^Baillie2021]: if V_{n+1} == 2 Q, report `n` as prime.
498        if vk == q.clone() + &q {
499            return Primality::ProbablyPrime;
500        }
501    }
502
503    Primality::Composite
504}
505
506#[cfg(test)]
507mod tests {
508
509    use alloc::format;
510
511    use crypto_bigint::{Integer, Odd, Uint, Word, U128, U64};
512
513    #[cfg(feature = "tests-exhaustive")]
514    use num_prime::nt_funcs::is_prime64;
515
516    use super::{decompose, lucas_test, AStarBase, BruteForceBase, LucasBase, LucasCheck, SelfridgeBase};
517    use crate::hazmat::{primes, pseudoprimes, Primality};
518
519    #[test]
520    fn bases_derived_traits() {
521        assert_eq!(format!("{SelfridgeBase:?}"), "SelfridgeBase");
522        assert_eq!(SelfridgeBase.clone(), SelfridgeBase);
523        assert_eq!(format!("{AStarBase:?}"), "AStarBase");
524        assert_eq!(AStarBase.clone(), AStarBase);
525        assert_eq!(format!("{BruteForceBase:?}"), "BruteForceBase");
526        assert_eq!(BruteForceBase.clone(), BruteForceBase);
527
528        assert_eq!(format!("{:?}", LucasCheck::Strong), "Strong");
529        assert_eq!(LucasCheck::Strong.clone(), LucasCheck::Strong);
530    }
531
532    #[test]
533    fn base_for_square() {
534        // We can't find a base with Jacobi symbol = -1 for a square,
535        // check that it is handled properly.
536        let num = Odd::new(U64::from(131u32).square()).unwrap();
537        assert_eq!(SelfridgeBase.generate(&num), Err(Primality::Composite));
538        assert_eq!(AStarBase.generate(&num), Err(Primality::Composite));
539        assert_eq!(BruteForceBase.generate(&num), Err(Primality::Composite));
540    }
541
542    #[test]
543    fn base_early_quit() {
544        // 5 is flagged as prime at the base generation stage
545        assert_eq!(
546            BruteForceBase.generate(&Odd::new(U64::from(5u32)).unwrap()),
547            Err(Primality::Prime)
548        )
549    }
550
551    #[test]
552    fn gcd_check() {
553        // Test that `gcd(2QD, n) == 1` is checked after the base is found.
554        // All the bases already produce D such that `gcd(D, n) == 1`.
555        // We need a special test base generator to test the situation where
556        // `gcd(n, Q) != 1` and `n > Q`, because it just doesn't seem to happen normally
557        // for "production" bases.
558
559        struct TestBase;
560
561        impl LucasBase for TestBase {
562            fn generate<T: Integer>(&self, _n: &Odd<T>) -> Result<(Word, Word, bool), Primality> {
563                Ok((5, 5, false))
564            }
565        }
566
567        assert_eq!(
568            lucas_test(Odd::new(U64::from(15u32)).unwrap(), TestBase, LucasCheck::Strong),
569            Primality::Composite
570        );
571    }
572
573    #[test]
574    fn decomposition() {
575        assert_eq!(
576            decompose(&Odd::new(U128::MAX).unwrap()),
577            (128, Odd::new(U128::ONE).unwrap())
578        );
579        assert_eq!(
580            decompose(&Odd::new(U128::ONE).unwrap()),
581            (1, Odd::new(U128::ONE).unwrap())
582        );
583        assert_eq!(
584            decompose(&Odd::new(U128::from(7766015u32)).unwrap()),
585            (15, Odd::new(U128::from(237u32)).unwrap())
586        );
587    }
588
589    fn is_slpsp(num: u32) -> bool {
590        pseudoprimes::STRONG_LUCAS.iter().any(|x| *x == num)
591    }
592
593    fn is_aeslpsp(num: u32) -> bool {
594        pseudoprimes::ALMOST_EXTRA_STRONG_LUCAS.iter().any(|x| *x == num)
595    }
596
597    fn is_eslpsp(num: u32) -> bool {
598        pseudoprimes::EXTRA_STRONG_LUCAS.iter().any(|x| *x == num)
599    }
600
601    fn is_vpsp(num: u32) -> bool {
602        pseudoprimes::LUCAS_V.iter().any(|x| *x == num)
603    }
604
605    fn test_composites_selfridge(numbers: &[u32], expected_result: bool) {
606        for num in numbers.iter() {
607            let false_positive = is_slpsp(*num);
608            let actual_expected_result = if false_positive { true } else { expected_result };
609
610            // Test both single-limb and multi-limb, just in case.
611            assert_eq!(
612                lucas_test(
613                    Odd::new(Uint::<1>::from(*num)).unwrap(),
614                    SelfridgeBase,
615                    LucasCheck::Strong
616                )
617                .is_probably_prime(),
618                actual_expected_result
619            );
620            assert_eq!(
621                lucas_test(
622                    Odd::new(Uint::<2>::from(*num)).unwrap(),
623                    SelfridgeBase,
624                    LucasCheck::Strong
625                )
626                .is_probably_prime(),
627                actual_expected_result
628            );
629        }
630    }
631
632    fn test_composites_a_star(numbers: &[u32], expected_result: bool) {
633        for num in numbers.iter() {
634            let false_positive = is_vpsp(*num);
635            let actual_expected_result = if false_positive { true } else { expected_result };
636
637            // Test both single-limb and multi-limb, just in case.
638            assert_eq!(
639                lucas_test(Odd::new(Uint::<1>::from(*num)).unwrap(), AStarBase, LucasCheck::LucasV).is_probably_prime(),
640                actual_expected_result
641            );
642            assert_eq!(
643                lucas_test(Odd::new(Uint::<2>::from(*num)).unwrap(), AStarBase, LucasCheck::LucasV).is_probably_prime(),
644                actual_expected_result
645            );
646        }
647    }
648
649    fn test_composites_brute_force(numbers: &[u32], almost_extra: bool, expected_result: bool) {
650        for num in numbers.iter() {
651            let false_positive = if almost_extra {
652                is_aeslpsp(*num)
653            } else {
654                is_eslpsp(*num)
655            };
656            let actual_expected_result = if false_positive { true } else { expected_result };
657            let check = if almost_extra {
658                LucasCheck::AlmostExtraStrong
659            } else {
660                LucasCheck::ExtraStrong
661            };
662
663            // Test both single-limb and multi-limb, just in case.
664            assert_eq!(
665                lucas_test(Odd::new(Uint::<1>::from(*num)).unwrap(), BruteForceBase, check).is_probably_prime(),
666                actual_expected_result,
667                "Brute force base, n = {num}, almost_extra = {almost_extra}",
668            );
669            assert_eq!(
670                lucas_test(Odd::new(Uint::<2>::from(*num)).unwrap(), BruteForceBase, check).is_probably_prime(),
671                actual_expected_result
672            );
673        }
674    }
675
676    #[test]
677    fn strong_fibonacci_pseudoprimes() {
678        // Can't use `test_composites()` since `STRONG_FIBONACCI` is `U64`.
679        // Good thing we don't need to test for intersection
680        // with `EXTRA_STRONG_LUCAS` or `STRONG_LUCAS` - there's none.
681        for num in pseudoprimes::STRONG_FIBONACCI.iter() {
682            assert!(!lucas_test(Odd::new(*num).unwrap(), SelfridgeBase, LucasCheck::Strong).is_probably_prime());
683            assert!(!lucas_test(Odd::new(*num).unwrap(), BruteForceBase, LucasCheck::ExtraStrong).is_probably_prime());
684        }
685    }
686
687    #[test]
688    fn fibonacci_pseudoprimes() {
689        test_composites_selfridge(pseudoprimes::FIBONACCI, false);
690        test_composites_a_star(pseudoprimes::FIBONACCI, false);
691        test_composites_brute_force(pseudoprimes::FIBONACCI, false, false);
692        test_composites_brute_force(pseudoprimes::FIBONACCI, true, false);
693    }
694
695    #[test]
696    fn bruckman_lucas_pseudoprimes() {
697        test_composites_selfridge(pseudoprimes::BRUCKMAN_LUCAS, false);
698        test_composites_a_star(pseudoprimes::BRUCKMAN_LUCAS, false);
699        test_composites_brute_force(pseudoprimes::BRUCKMAN_LUCAS, false, false);
700        test_composites_brute_force(pseudoprimes::BRUCKMAN_LUCAS, true, false);
701    }
702
703    #[test]
704    fn almost_extra_strong_lucas_pseudoprimes() {
705        test_composites_selfridge(pseudoprimes::ALMOST_EXTRA_STRONG_LUCAS, false);
706        test_composites_a_star(pseudoprimes::ALMOST_EXTRA_STRONG_LUCAS, false);
707
708        // Check for the difference between the almost extra strong and extra strong tests.
709        test_composites_brute_force(pseudoprimes::ALMOST_EXTRA_STRONG_LUCAS, false, false);
710        test_composites_brute_force(pseudoprimes::ALMOST_EXTRA_STRONG_LUCAS, true, true);
711    }
712
713    #[test]
714    fn extra_strong_lucas_pseudoprimes() {
715        test_composites_selfridge(pseudoprimes::EXTRA_STRONG_LUCAS, false);
716        test_composites_a_star(pseudoprimes::EXTRA_STRONG_LUCAS, false);
717
718        // These are the known false positives for the extra strong test
719        // with brute force base selection.
720        test_composites_brute_force(pseudoprimes::EXTRA_STRONG_LUCAS, false, true);
721    }
722
723    #[test]
724    fn lucas_pseudoprimes() {
725        test_composites_selfridge(pseudoprimes::LUCAS, false);
726        test_composites_a_star(pseudoprimes::LUCAS, false);
727        test_composites_brute_force(pseudoprimes::LUCAS, false, false);
728        test_composites_brute_force(pseudoprimes::LUCAS, true, false);
729    }
730
731    #[test]
732    fn strong_lucas_pseudoprimes() {
733        // These are the known false positives for the strong test
734        // with Selfridge base selection.
735        test_composites_selfridge(pseudoprimes::STRONG_LUCAS, true);
736
737        test_composites_a_star(pseudoprimes::STRONG_LUCAS, false);
738        test_composites_brute_force(pseudoprimes::STRONG_LUCAS, false, false);
739        test_composites_brute_force(pseudoprimes::STRONG_LUCAS, true, false);
740    }
741
742    #[test]
743    fn strong_pseudoprimes_base_2() {
744        // Cross-test against the pseudoprimes that circumvent the MR test base 2.
745        // We expect the Lucas test to correctly classify them as composites.
746        test_composites_selfridge(pseudoprimes::STRONG_BASE_2, false);
747        test_composites_a_star(pseudoprimes::STRONG_BASE_2, false);
748        test_composites_brute_force(pseudoprimes::STRONG_BASE_2, false, false);
749        test_composites_brute_force(pseudoprimes::STRONG_BASE_2, true, false);
750    }
751
752    #[test]
753    fn large_carmichael_number() {
754        let p = Odd::new(pseudoprimes::LARGE_CARMICHAEL_NUMBER).unwrap();
755        assert!(!lucas_test(p, SelfridgeBase, LucasCheck::Strong).is_probably_prime());
756        assert!(!lucas_test(p, AStarBase, LucasCheck::LucasV).is_probably_prime());
757        assert!(!lucas_test(p, BruteForceBase, LucasCheck::AlmostExtraStrong).is_probably_prime());
758        assert!(!lucas_test(p, BruteForceBase, LucasCheck::ExtraStrong).is_probably_prime());
759    }
760
761    fn test_large_primes<const L: usize>(nums: &[Uint<L>]) {
762        for num in nums {
763            let num = Odd::new(*num).unwrap();
764            assert!(lucas_test(num, SelfridgeBase, LucasCheck::Strong).is_probably_prime());
765            assert!(lucas_test(num, AStarBase, LucasCheck::LucasV).is_probably_prime());
766            assert!(lucas_test(num, BruteForceBase, LucasCheck::AlmostExtraStrong).is_probably_prime());
767            assert!(lucas_test(num, BruteForceBase, LucasCheck::ExtraStrong).is_probably_prime());
768        }
769    }
770
771    #[test]
772    fn large_primes() {
773        test_large_primes(primes::PRIMES_128);
774        test_large_primes(primes::PRIMES_256);
775        test_large_primes(primes::PRIMES_384);
776        test_large_primes(primes::PRIMES_512);
777        test_large_primes(primes::PRIMES_1024);
778    }
779
780    #[test]
781    fn test_lucas_v_pseudoprimes() {
782        for num in pseudoprimes::LARGE_LUCAS_V {
783            let num = Odd::new(*num).unwrap();
784            // These are false positives for Lucas-V test
785            assert!(lucas_test(num, AStarBase, LucasCheck::LucasV).is_probably_prime());
786
787            // These tests should work correctly
788            assert!(!lucas_test(num, SelfridgeBase, LucasCheck::Strong).is_probably_prime());
789            assert!(!lucas_test(num, BruteForceBase, LucasCheck::AlmostExtraStrong).is_probably_prime());
790            assert!(!lucas_test(num, BruteForceBase, LucasCheck::ExtraStrong).is_probably_prime());
791        }
792    }
793
794    #[test]
795    fn corner_cases() {
796        // By convention, 1 is composite. That's what `num-prime` returns.
797        let res = lucas_test(
798            Odd::new(U64::ONE).unwrap(),
799            BruteForceBase,
800            LucasCheck::AlmostExtraStrong,
801        );
802        assert_eq!(res, Primality::Composite);
803    }
804
805    #[cfg(feature = "tests-exhaustive")]
806    #[test]
807    fn exhaustive() {
808        // Test all the odd numbers up to the limit where we know the false positives,
809        // and compare the results with the reference.
810        for num in (3..pseudoprimes::EXHAUSTIVE_TEST_LIMIT).step_by(2) {
811            let res_ref = is_prime64(num.into());
812
813            let eslpsp = is_eslpsp(num);
814            let aeslpsp = is_aeslpsp(num);
815            let slpsp = is_slpsp(num);
816            let vpsp = is_vpsp(num);
817
818            let odd_num = Odd::new(Uint::<1>::from(num)).unwrap();
819
820            let res = lucas_test(odd_num, BruteForceBase, LucasCheck::AlmostExtraStrong).is_probably_prime();
821            let expected = aeslpsp || res_ref;
822            assert_eq!(
823                res, expected,
824                "Brute force base, almost extra strong: n={num}, expected={expected}, actual={res}",
825            );
826
827            let res = lucas_test(odd_num, BruteForceBase, LucasCheck::ExtraStrong).is_probably_prime();
828            let expected = eslpsp || res_ref;
829            assert_eq!(
830                res, expected,
831                "Brute force base: n={num}, expected={expected}, actual={res}",
832            );
833
834            let res = lucas_test(odd_num, SelfridgeBase, LucasCheck::Strong).is_probably_prime();
835            let expected = slpsp || res_ref;
836            assert_eq!(
837                res, expected,
838                "Selfridge base: n={num}, expected={expected}, actual={res}",
839            );
840
841            let res = lucas_test(odd_num, AStarBase, LucasCheck::LucasV).is_probably_prime();
842            let expected = vpsp || res_ref;
843
844            assert_eq!(
845                res, expected,
846                "A* base, Lucas-V: n={num}, expected={expected}, actual={res}",
847            );
848        }
849    }
850}
crypto_primes/hazmat/lucas.rs

crypto_primes/hazmat/
lucas.rs
