Draft impl of JP06 pseudo-uniform sieve

Cargo.toml

@@ -18,8 +18,10 @@ openssl = { version = "0.10.39", optional = true, features = ["vendored"] }
 rug = { version = "1.26", default-features = false, features = ["integer"], optional = true }
 glass_pumpkin = { version = "1", optional = true }
 rayon = { version = "1", optional = true }
+tracing = { version = "0.1", default-features = false }
 
 [dev-dependencies]
+rand = "0.8"
 # need `crypto-bigint` with `alloc` to test `BoxedUint`
 crypto-bigint = { version = "0.6.0-rc.6", default-features = false, features = ["alloc"] }
 rand_chacha = "0.3"
@@ -30,6 +32,8 @@ num-integer = "0.1"
 proptest = "1"
 num-prime = "0.4.3"
 num_cpus = "1.16"
+test-log = { version = "0.2.16", default-features = false, features = ["trace", "color"] }
+statrs = "0.17.1"
 
 [features]
 default = ["default-rng"]
@@ -49,3 +53,10 @@ rustdoc-args = ["--cfg", "docsrs"]
 [[bench]]
 name = "bench"
 harness = false
+
+[profile.bench]
+lto = true
+codegen-units = 1
+
+[[example]]
+name = "distribution"

benches/bench.rs

@@ -1,7 +1,7 @@
 use core::num::NonZero;
 
 use criterion::{criterion_group, criterion_main, BatchSize, Criterion};
-use crypto_bigint::{nlimbs, BoxedUint, Integer, Odd, RandomBits, Uint, U1024, U128, U256};
+use crypto_bigint::{nlimbs, BoxedUint, Integer, Odd, RandomBits, Uint, U1024, U128, U2048, U256, U512};
 use rand_chacha::ChaCha8Rng;
 use rand_core::{CryptoRngCore, OsRng, SeedableRng};
 
@@ -39,6 +39,35 @@ fn make_presieved_num<const L: usize>(rng: &mut impl CryptoRngCore) -> Odd<Uint<
     Odd::new(sieve.next().unwrap()).unwrap()
 }
 
+fn bench_uniform_sieve(c: &mut Criterion) {
+    use crypto_primes::hazmat::UniformSieve;
+    let mut group = c.benchmark_group("Uniform sieve");
+
+    let mut rng = make_rng();
+    group.bench_function("(U128) Random prime", |b| {
+        b.iter(|| U128::generate_prime_with_rng(&mut rng));
+    });
+
+    let mut rng = make_rng();
+    group.bench_function("(U256) Random prime", |b| {
+        b.iter(|| U256::generate_prime_with_rng(&mut rng));
+    });
+
+    let mut rng = make_rng();
+    group.bench_function("(U512) Random prime", |b| {
+        b.iter(|| U256::generate_prime_with_rng(&mut rng));
+    });
+
+    let mut rng = make_rng();
+    group.bench_function("(U1024) Random prime", |b| {
+        b.iter(|| U1024::generate_prime_with_rng(&mut rng));
+    });
+
+    let mut rng = make_rng();
+    group.bench_function("(U2048) Random prime", |b| {
+        b.iter(|| U2048::generate_prime_with_rng(&mut rng));
+    });
+}
 fn bench_sieve(c: &mut Criterion) {
     let mut group = c.benchmark_group("Sieve");
 
@@ -227,11 +256,26 @@ fn bench_presets(c: &mut Criterion) {
         b.iter(|| generate_prime_with_rng::<U128>(&mut rng, 128))
     });
 
+    let mut rng = make_rng();
+    group.bench_function("(U256) Random prime", |b| {
+        b.iter(|| generate_prime_with_rng::<U256>(&mut rng, 128))
+    });
+
+    let mut rng = make_rng();
+    group.bench_function("(U512) Random prime", |b| {
+        b.iter(|| generate_prime_with_rng::<U512>(&mut rng, 128))
+    });
+
     let mut rng = make_rng();
     group.bench_function("(U1024) Random prime", |b| {
         b.iter(|| generate_prime_with_rng::<U1024>(&mut rng, 1024))
     });
 
+    let mut rng = make_rng();
+    group.bench_function("(U2048) Random prime", |b| {
+        b.iter(|| generate_prime_with_rng::<U2048>(&mut rng, 1024))
+    });
+
     let mut rng = make_rng();
     group.bench_function("(U128) Random safe prime", |b| {
         b.iter(|| generate_safe_prime_with_rng::<U128>(&mut rng, 128))
@@ -523,6 +567,7 @@ fn bench_glass_pumpkin(_c: &mut Criterion) {}
 
 criterion_group!(
     benches,
+    bench_uniform_sieve,
     bench_sieve,
     bench_miller_rabin,
     bench_lucas,

examples/distribution.rs

@@ -0,0 +1,90 @@
+use crypto_bigint::{FixedInteger, NonZero, RandomBits, RandomMod};
+use crypto_primes::hazmat::UniformSieve;
+use statrs::distribution::{ChiSquared, ContinuousCDF};
+
+/// Chi-squared statistical test for some collection of primes.
+
+fn main() {
+    let sample_count = 5000;
+    let num_intervals = 20;
+    type T = crypto_bigint::U512;
+
+    let end = T::MAX;
+    let (counts, elapsed) = distribution::<T>(sample_count, num_intervals, end);
+    let expected_count = sample_count as f64 / num_intervals as f64;
+
+    // Chi-Square test for uniformity
+    let (chi_square_stat, p_value) = chi_square_test(&counts, expected_count);
+
+    println!("Prime Counts in Each Interval: {:?}", counts);
+    println!("Collected {sample_count} primes in {elapsed:?}");
+    println!("Expected Count per Interval: {:.2}", expected_count);
+    println!("Chi-Square Statistic: {:.2}", chi_square_stat);
+    println!("P-Value: {:.4}", p_value);
+
+    if p_value < 0.05 {
+        println!("The prime distribution significantly deviates from uniformity.");
+    } else {
+        println!("The prime distribution is close to uniform.");
+    }
+}
+
+// Divide interval and count primes in each subinterval. Returns stats and the time elapsed.
+fn distribution<T>(num_primes: usize, intervals: usize, end: T) -> (Vec<u64>, std::time::Duration)
+where
+    T: FixedInteger + RandomBits + RandomMod,
+    T: UniformSieve<T>,
+{
+    let now = std::time::Instant::now();
+    let primes: Vec<T> = collect_primes(num_primes);
+    let elapsed = now.elapsed();
+    println!("primes={:?}", &primes[0..5]);
+    let start = T::ZERO;
+    assert!(end > start);
+    let interval_width = (end - start).div(NonZero::new(T::from(intervals as u64)).unwrap());
+    let interval_width_nz = NonZero::new(interval_width).unwrap();
+    println!("intervals_width={interval_width:?}");
+    // Initialize counters for each interval
+    let mut counts = vec![0; intervals];
+
+    // Count primes in each interval
+    for &prime in &primes {
+        let interval = (prime - start).div(interval_width_nz).as_ref()[0].0 as usize;
+        if interval < intervals {
+            counts[interval] += 1;
+        } else {
+            println!("interval={interval} is outside the intervals count {intervals}");
+        }
+    }
+
+    (counts, elapsed)
+}
+
+// Collect primes
+fn collect_primes<T>(num_primes: usize) -> Vec<T>
+where
+    T: FixedInteger + RandomBits + RandomMod,
+    T: UniformSieve<T>,
+{
+    (0..num_primes).fold(Vec::with_capacity(num_primes), |mut acc, _| {
+        acc.push(T::generate_prime());
+        acc
+    })
+}
+
+// Calculate Chi-Square statistic and p-value
+fn chi_square_test(counts: &[u64], expected: f64) -> (f64, f64) {
+    let chi_square_stat: f64 = counts
+        .iter()
+        .map(|&count| {
+            let diff = count as f64 - expected;
+            diff * diff / expected
+        })
+        .sum();
+
+    // Calculate p-value
+    let chi_square_dist = ChiSquared::new((counts.len() - 1) as f64).unwrap();
+    let p_value = 1.0 - chi_square_dist.cdf(chi_square_stat);
+
+    (chi_square_stat, p_value)
+}

src/hazmat.rs

@@ -5,16 +5,19 @@ mod gcd;
 mod jacobi;
 mod lucas;
 mod miller_rabin;
-mod precomputed;
+pub(crate) mod precomputed;
 #[cfg(test)]
 pub(crate) mod primes;
 #[cfg(test)]
 pub(crate) mod pseudoprimes;
 mod sieve;
+mod uniform_sieve;
 
 pub use lucas::{lucas_test, AStarBase, BruteForceBase, LucasBase, LucasCheck, SelfridgeBase};
 pub use miller_rabin::MillerRabin;
+
 pub use sieve::{random_odd_integer, Sieve};
+pub use uniform_sieve::UniformSieve;
 
 /// Possible results of various primality tests.
 #[derive(Copy, Clone, Debug, PartialEq, Eq)]

src/hazmat/gcd.rs

@@ -42,7 +42,7 @@ pub(crate) fn gcd_vartime<T: Integer>(n: &T, m: NonZero<Word>) -> Word {
 // As GCD is commutative `gcd(n, m) = gcd(m, n)` those identities still apply if the operands are swapped.
 //
 // [1]: https://en.wikipedia.org/wiki/Binary_GCD_algorithm
-fn binary_gcd(mut n: Word, mut m: Word) -> Word {
+pub(crate) fn binary_gcd(mut n: Word, mut m: Word) -> Word {
     // Using identities 2 and 3:
     // gcd(2ⁱn, 2ʲm) = 2ᵏ gcd(n, m) with n, m odd and k = min(i, j)
     // 2ᵏ is the greatest power of two that divides both 2ⁱn and 2ʲm