Merge pull request #20 from fjarri/lucas-sequence

Lucas sequence optimization

Author: fjarri

.github/workflows/crypto-primes.yml

@@ -111,7 +111,7 @@ jobs:
           components: clippy
           override: true
           profile: minimal
-      - run: cargo clippy --all --all-features --tests -- -D warnings
+      - run: cargo clippy --all --all-features --tests --benches -- -D warnings
 
   rustfmt:
     runs-on: ubuntu-latest

CHANGELOG.md

@@ -4,6 +4,21 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 
+## [Unreleased]
+
+### Added
+
+- An alternative propagation method for Lucas sequences improving the performance of `lucas_test()`. ([#20])
+
+
+### Fixed
+
+- Some mistakes in the description of Lucas checks (the logic itself was fine). ([#20])
+
+
+[#20]: https://github.com/nucypher/rust-umbral/pull/20
+
+
 ## [0.2.0] - 2023-03-06
 
 ### Changed

Cargo.toml

@@ -9,7 +9,7 @@ readme = "README.md"
 categories = ["cryptography", "no-std"]
 
 [dependencies]
-crypto-bigint = { version = "0.5", default-features = false, features = ["rand_core"] }
+crypto-bigint = { version = "0.5.2", default-features = false, features = ["rand_core"] }
 rand_core = { version = "0.6", default-features = false }
 getrandom = { version = "0.2", optional = true, default-features = false, features = ["js"] }
 openssl = { version = "0.10", optional = true, features = ["vendored"] }

benches/bench.rs

@@ -1,25 +1,26 @@
-use rand_core::OsRng;
-
 use criterion::{
     criterion_group, criterion_main, measurement::Measurement, BenchmarkGroup, Criterion,
 };
 use crypto_bigint::{U1024, U128};
+use rand_chacha::ChaCha8Rng;
+use rand_core::{OsRng, SeedableRng};
 
 use crypto_primes::{
     hazmat::{
-        lucas_test, random_odd_uint, BruteForceBase, LucasCheck, MillerRabin, SelfridgeBase, Sieve,
+        lucas_test, random_odd_uint, AStarBase, BruteForceBase, LucasCheck, MillerRabin,
+        SelfridgeBase, Sieve,
     },
     safe_prime_with_rng,
 };
 
-fn bench_sieve<'a, M: Measurement>(group: &mut BenchmarkGroup<'a, M>) {
+fn bench_sieve<M: Measurement>(group: &mut BenchmarkGroup<'_, M>) {
     let start: U1024 = random_odd_uint(&mut OsRng, 1024);
     group.bench_function("(U1024) Sieve, 1000 samples", |b| {
         b.iter(|| Sieve::new(&start, 1024, false).take(1000).for_each(drop))
     });
 }
 
-fn bench_miller_rabin<'a, M: Measurement>(group: &mut BenchmarkGroup<'a, M>) {
+fn bench_miller_rabin<M: Measurement>(group: &mut BenchmarkGroup<'_, M>) {
     let start: U1024 = random_odd_uint(&mut OsRng, 1024);
     group.bench_function("(U1024) Miller-Rabin creation", |b| {
         b.iter(|| {
@@ -40,15 +41,28 @@ fn bench_miller_rabin<'a, M: Measurement>(group: &mut BenchmarkGroup<'a, M>) {
     );
 }
 
-fn bench_lucas<'a, M: Measurement>(group: &mut BenchmarkGroup<'a, M>) {
-    let start: U1024 = random_odd_uint(&mut OsRng, 1024);
+fn bench_lucas<M: Measurement>(group: &mut BenchmarkGroup<'_, M>) {
+    let mut rng = ChaCha8Rng::from_seed(*b"01234567890123456789012345678901");
+
+    let start: U1024 = random_odd_uint(&mut rng, 1024);
     let mut sieve = Sieve::new(&start, 1024, false);
-    group.bench_function("(U1024) Sieve + Lucas test (Selfridge base)", |b| {
+    group.bench_function(
+        "(U1024) Sieve + Lucas test (Selfridge base, strong check)",
+        |b| {
+            b.iter(|| {
+                lucas_test(&sieve.next().unwrap(), SelfridgeBase, LucasCheck::Strong);
+            })
+        },
+    );
+
+    let mut sieve = Sieve::new(&start, 1024, false);
+    group.bench_function("(U1024) Sieve + Lucas test (A* base, Lucas-V check)", |b| {
         b.iter(|| {
-            lucas_test(&sieve.next().unwrap(), SelfridgeBase, LucasCheck::Strong);
+            lucas_test(&sieve.next().unwrap(), AStarBase, LucasCheck::LucasV);
         })
     });
 
+    let mut sieve = Sieve::new(&start, 1024, false);
     group.bench_function(
         "(U1024) Sieve + Lucas test (brute force base, almost extra strong)",
         |b| {
@@ -62,6 +76,7 @@ fn bench_lucas<'a, M: Measurement>(group: &mut BenchmarkGroup<'a, M>) {
         },
     );
 
+    let mut sieve = Sieve::new(&start, 1024, false);
     group.bench_function(
         "(U1024) Sieve + Lucas test (brute force base, extra strong)",
         |b| {
@@ -102,7 +117,7 @@ fn bench_primality_tests(c: &mut Criterion) {
     group.finish();
 }
 
-fn bench_presets<'a, M: Measurement>(group: &mut BenchmarkGroup<'a, M>) {
+fn bench_presets<M: Measurement>(group: &mut BenchmarkGroup<'_, M>) {
     group.bench_function("(U128) Random safe prime", |b| {
         b.iter(|| {
             let p: U128 = safe_prime_with_rng(&mut OsRng, 128);

src/hazmat/lucas.rs

@@ -1,7 +1,7 @@
 //! Lucas primality test.
 use crypto_bigint::{
     modular::runtime_mod::{DynResidue, DynResidueParams},
-    Integer, Invert, Limb, Uint, Word,
+    Integer, Limb, Uint, Word,
 };
 
 use super::{
@@ -197,7 +197,10 @@ fn decompose<const L: usize>(n: &Uint<L>) -> (u32, Uint<L>) {
 #[derive(Copy, Clone, Debug, PartialEq, Eq)]
 pub enum LucasCheck {
     /// Introduced by Baillie & Wagstaff[^Baillie1980].
-    /// If any of `V(d*2^r) == 0` for `0 <= r < s`, and `U(d) == 0`, report the number as prime.
+    /// If either of the following is true:
+    /// - any of `V(d*2^r) == 0` for `0 <= r < s`,
+    /// - `U(d) == 0`,
+    /// report the number as prime.
     ///
     /// If the base is [`SelfridgeBase`], known false positives constitute OEIS:A217255[^A217255].
     ///
@@ -210,8 +213,14 @@ pub enum LucasCheck {
     /// [^A217255]: <https://oeis.org/A217255>
     Strong,
 
-    /// If any of `V(d*2^r) == 0` for `0 <= r < s`, and `V(d) == ±2` report the number as prime.
-    /// The second condition is only checked if `Q == 1`, otherwise it is considered to be true.
+    /// A [`LucasCheck::ExtraStrong`] without checking for `U(d)`.
+    /// That is, if either of the following is true:
+    /// - any of `V(d*2^r) == 0` for `0 <= r < s`,
+    /// - `V(d) == ±2`,
+    /// report the number as prime.
+    ///
+    /// Note: the second condition is only checked if `Q == 1`,
+    /// otherwise it is considered to be true.
     ///
     /// If the base is [`BruteForceBase`], some known false positives
     /// are listed by Jacobsen[^Jacobsen].
@@ -225,7 +234,9 @@ pub enum LucasCheck {
     AlmostExtraStrong,
 
     /// Introduced by Mo[^Mo1993], and also described by Grantham[^Grantham2001].
-    /// If [`LucasCheck::Strong`] check passes, and `V(d) == ±2`,
+    /// If either of the following is true:
+    /// - any of `V(d*2^r) == 0` for `0 <= r < s`,
+    /// - `U(d) == 0` and `V(d) == ±2`,
     /// report the number as prime.
     ///
     /// Note that this check only differs from [`LucasCheck::Strong`] if `Q == 1`.
@@ -335,102 +346,89 @@ pub fn lucas_test<const L: usize>(
         DynResidue::<L>::new(&Uint::<L>::from(p), params)
     };
 
-    // Compute d-th element of Lucas sequence V_d(P, Q), where:
-    //
-    //  V_0 = 2
-    //  V_1 = P
-    //  V_k = P V_{k-1} - Q V_{k-2}.
+    // Compute d-th element of Lucas sequence (U_d(P, Q), V_d(P, Q)), where:
     //
-    // In general V(k) = α^k + β^k, where α and β are roots of x^2 - Px + Q.
-    // [^Crandall2005], eq. (3.14) observe that for 0 <= j <= k,
+    // V_0 = 2
+    // U_0 = 1
     //
-    //  V_{j+k} = V_j V_k - Q^j * V_{k-j}.
+    // U_{2k} = U_k V_k
+    // V_{2k} = V_k^2 - 2 Q^k
     //
-    // So in particular, to quickly double the subscript:
-    //
-    //  V_{2k} = V_k^2 - 2 * Q^k
-    //  V_{2k+1} = V_k V_{k+1} - Q^k
+    // U_{k+1} = (P U_k + V_k) / 2
+    // V_{k+1} = (D U_k + P V_k) / 2
     //
+    // (The propagation method is due to [^Baillie2021], Eqs. 13, 14, 16, 17)
     // We can therefore start with k=0 and build up to k=d in log2(d) steps.
 
+    // Starting with k = 0
     let mut vk = two; // keeps V_k
-    let mut vk1 = p; // keeps V_{k+1}
+    let mut uk = DynResidue::<L>::zero(params); // keeps U_k
     let mut qk = one; // keeps Q^k
-    let mut qk_times_p = if p_is_one { one } else { p }; // keeps P Q^{k}
+
+    // D in Montgomery representation - note that it can be negative.
+    let abs_d = DynResidue::<L>::new(&Uint::<L>::from(discriminant.abs_diff(0)), params);
+    let d_m = if discriminant < 0 { -abs_d } else { abs_d };
 
     for i in (0..d.bits_vartime()).rev() {
-        if d.bit_vartime(i) {
-            // k' = 2k+1
+        // k' = k * 2
 
-            // V_k' = V_{2k+1} = V_k V_{k+1} - P Q^k
-            vk = vk * vk1 - qk_times_p;
+        let u_2k = uk * vk;
+        let v_2k = vk.square() - (qk + qk);
+        let q_2k = qk.square();
 
-            // V_{k'+1} = V_{2k+2} = V_{k+1}^2 - 2 Q^{k+1}
-            let qk1 = qk * q; // Q^{k+1}
-            let two_qk1 = if q_is_one { two } else { qk1 + qk1 }; // 2 Q^{k+1}
-            vk1 = vk1.square() - two_qk1;
-            qk *= qk1;
-        } else {
-            // k' = 2k
+        uk = u_2k;
+        vk = v_2k;
+        qk = q_2k;
 
-            // V_{k'+1} = V_{2k+1} = V_k V_{k+1} - P Q^k
-            vk1 = vk * vk1 - qk_times_p;
+        if d.bit_vartime(i) {
+            // k' = k + 1
 
-            // V_k' = V_{2k} = V_k^2 - 2 Q^k
-            let two_qk = if q_is_one { two } else { qk + qk }; // 2 Q^k
-            vk = vk.square() - two_qk;
-            qk = qk.square();
-        }
+            let (p_uk, p_vk) = if p_is_one { (uk, vk) } else { (p * uk, p * vk) };
 
-        if p_is_one {
-            qk_times_p = qk;
-        } else {
-            qk_times_p = qk * p;
+            let u_k1 = (p_uk + vk).div_by_2();
+            let v_k1 = (d_m * uk + p_vk).div_by_2();
+            let q_k1 = qk * q;
+
+            uk = u_k1;
+            vk = v_k1;
+            qk = q_k1;
         }
     }
 
-    // Now k=d, so vk = V_d, vk_1 = V_{d+1}.
+    // Now k=d, so vk = V_d and uk = U_d.
 
-    // Extra strong check (from [^Mo1993]): `V_d == ±2 mod n`.
-    // Do it first since it is cheap.
-    //
-    // Note that it only applies if Q = 1, since it is a consequence
-    // of a property of Lucas series: V_k^2 - 4 Q^k = D U_k^2 mod n.
-    // If Q = 1 we can easily decompose the left side of the equation leading to the check above.
-    let vk_equals_two = if q_is_one {
-        vk == two || vk == minus_two
-    } else {
-        true
-    };
+    // Check for the first sufficient condition in various strong checks.
 
-    if vk_equals_two {
-        // Strong check:`U_d == 0 mod n`.
-        // As suggested by [^Jacobsen], apply Eq. (3.13) from [^Crandall2005]:
+    if check == LucasCheck::Strong && uk == zero {
+        // Strong check: `U_d == 0 mod n`.
+        return Primality::ProbablyPrime;
+    } else if check == LucasCheck::ExtraStrong || check == LucasCheck::AlmostExtraStrong {
+        // Extra strong check (from [^Mo1993]): `V_d == ±2 mod n` and `U_d == 0 mod n`.
         //
-        //  U_k = D^{-1} (2 V_{k+1} - P V_k)
+        // Note that the first identity only applies if `Q = 1`, since it is a consequence
+        // of a property of Lucas series: `V_k^2 - 4 Q^k = D U_k^2 mod n`.
+        // If `Q = 1` we can easily decompose the left side of the equation
+        // leading to the check above.
         //
-        // Some implementations just test for 2 V_{k+1} == P V_{k},
-        // but we don't have any reference pseudoprime lists for this, so we are not doing it.
-        if check == LucasCheck::Strong || check == LucasCheck::ExtraStrong {
-            let abs_d = DynResidue::<L>::new(&Uint::<L>::from(discriminant.abs_diff(0)), params);
-            let d_m = if discriminant < 0 { -abs_d } else { abs_d };
-            // `d` is guaranteed non-zero by construction, so we can safely unwrap
-            let inv_d = <DynResidue<L> as Invert>::invert(&d_m).unwrap();
-
-            let vk_times_p = if p_is_one { vk } else { vk * p };
-            let uk = inv_d * (vk1 + vk1 - vk_times_p);
-
-            if uk == zero {
-                return Primality::ProbablyPrime;
-            }
-        } else {
-            // This is "almost extra strong check": we only checked for `V_d` earlier.
-            if check == LucasCheck::AlmostExtraStrong {
-                return Primality::ProbablyPrime;
-            }
+        // If `Q != 1` we just consider it passed (we don't have a corresponding
+        // pseudoprime list anyway).
+
+        let vk_equals_two = !q_is_one || (vk == two || vk == minus_two);
+
+        if check == LucasCheck::ExtraStrong && uk == zero && vk_equals_two {
+            return Primality::ProbablyPrime;
+        }
+
+        // "Almost extra strong" check skips the `U_d` check.
+        // Since we have `U_d` anyway, it does not improve performance,
+        // so it is only here for testing purposes, since we have a corresponding pseudoprime list.
+        if check == LucasCheck::AlmostExtraStrong && vk_equals_two {
+            return Primality::ProbablyPrime;
         }
     }
 
+    // Second sufficient condition requires further propagating `V_k` up to `V_{n+1}`.
+
     // Check if V_{2^t d} == 0 mod n for some 0 <= t < s.
     // (unless we're in Lucas-V mode, then we just propagate V_k)
 
@@ -446,7 +444,7 @@ pub fn lucas_test<const L: usize>(
         }
 
         // k' = 2k
-        // V(k') = V(2k) = V(k)² - 2 * Q^k
+        // V_{k'} = V_k^2 - 2 Q^k
         vk = vk * vk - qk - qk;
 
         if check != LucasCheck::LucasV && vk == zero {