Increase the Prime wheel size

- From a 30 spoke wheel before (2*3*5)
- Up to a 210 spokes wheel (2*3*5*7)
- Add more tests
- Improve and add benchmarks
- Add performance numbers for the 210 spoke wheel to the README
- Bump version to `0.5.1`

Cargo.toml

@@ -1,7 +1,7 @@
 [package]
 name = "prime-factor"
 description = "A prime number factorizer written in Rust"
-version = "0.5.0"
+version = "0.5.1"
 authors = ["Stefan Lindblad <[email protected]>"]
 license = "LGPL-2.1-or-later"
 edition = "2021"

README.md

@@ -26,14 +26,21 @@ We want this generator to be fast and give reasonably good guesses for prime num
 
 ## Factorization performance
 
-On an old system (i7-6700):
+On an old system (i7-6700), with a 30-spoke wheel:
 - 32-bit, random number in about 32 ms and worst case 300 ms
 - 64-bit, random number in about 1.4 s (± 1.1 s) with worst case about 20 s
-- full test suite completes in about 7 minutes
+- full benchmark completes in about 7 minutes
 
-On a modern system (i7-12700):
+On a modern system (i7-12700), with a 30-spoke wheel:
 - 32-bit, random number in about 6.5 us and worst case in 68 us
 - 64-bit, random number in about 140 ms ([3 .. 340] ms) and worst case in 4.6 s
-- full test suite completes in less than 3 minutes
+- full benchmark completes in less than 3 minutes
+
+Modern system (i7-12700) with a 210-spoke Prime Wheel:
+- 2..8-bit prime numbers in 12..34 ns
+- 9..16-bit prime numbers in 33..252 ns
+- 17..32-bit prime numbers in 0.25..57 us, on average 5.6 us
+- 33..64-bit prime numbers in 0.056..3704 ms
+- 65+ bits prime numbers from 3.65 s
 
 The above numbers are taken from the included benchmark test, which you can run with the command: `cargo bench`. Note that it will take a few minutes to run the full suite, in which time you need to close all other applications and leave it unattended, to give the benchmark the most processing power possible.

benches/benchmark.rs

@@ -39,15 +39,25 @@ fn criterion_benchmark(c: &mut Criterion) {
 
     let mut fixed_grp = c.benchmark_group("worst-case");
     fixed_grp.sample_size(10);
+    fixed_grp.bench_function("prime-factor   lowest prime", |b| b.iter(||
+        pf_number(2)));
     fixed_grp.bench_function("prime-factor   highest 8-bit prime", |b| b.iter(||
         pf_number(251)));
+    fixed_grp.bench_function("prime-factor   lowest 9-bit prime", |b| b.iter(||
+        pf_number(257)));
     fixed_grp.bench_function("prime-factor   highest 16-bit prime", |b| b.iter(||
         pf_number(65521)));
+    fixed_grp.bench_function("prime-factor   lowest 17-bit prime", |b| b.iter(||
+        pf_number(65537)));
     fixed_grp.bench_function("prime-factor   highest 32-bit prime", |b| b.iter(||
         pf_number(4294967291)));
+    fixed_grp.bench_function("prime-factor   lowest 33-bit prime", |b| b.iter(||
+        pf_number(4294967311)));
     fixed_grp.measurement_time(Duration::new(60, 0));
     fixed_grp.bench_function("prime-factor   highest 64-bit prime", |b| b.iter(||
         pf_number(18446744073709551557)));
+    fixed_grp.bench_function("prime-factor   lowest 65-bit prime", |b| b.iter(||
+        pf_number(18446744073709551629)));
     fixed_grp.finish();
 
     let mut rand_grp = c.benchmark_group("random-numbers");

src/candidates.rs

@@ -1,10 +1,11 @@
+//! Implementations of Prime wheels for number factorization
+//! https://en.wikipedia.org/wiki/Wheel_factorization
 #![deny(unsafe_code)]
 #![allow(dead_code)]
 use genawaiter::yield_;
 use genawaiter::stack::producer_fn;
 
-/// Wheel factorization algorithm with base {2, 3, 5}
-/// https://en.wikipedia.org/wiki/Wheel_factorization
+/// Wheel factorization algorithm with base {2, 3, 5} (30 spokes)
 #[producer_fn(u128)]
 pub async fn prime_wheel_30() {
     yield_!(2);
@@ -18,3 +19,47 @@ pub async fn prime_wheel_30() {
         base += 30;
     }
 }
+
+pub const SPOKES: [u128; 48] = [
+    11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73,
+    79, 83, 89, 97, 101, 103, 107, 109, 113, 121, 127, 131, 137, 139, 143,
+    149, 151, 157, 163, 167, 169, 173, 179, 181, 187, 191, 193, 197, 199,
+    209, 211];
+
+/// Wheel factorization algorithm with base {2, 3, 5, 7} (210 spokes)
+#[producer_fn(u128)]
+pub async fn prime_wheel_210() {
+    yield_!(2);
+    yield_!(3);
+    yield_!(5);
+    yield_!(7);
+    // Valid spokes after removing multiples of 2, 3, 5, 7
+    let mut base = 0u128;
+    loop {
+        for n in SPOKES {
+            yield_!(base + n);
+        }
+        base += 210;
+    }
+}
+
+// Bit-map: 0x0200a28828228820a08a08820228a20208828828208a20a08a2802
+const PW210_BITMAP: [u8; 27] = [
+    0x02, 0x28, 0x8a, 0xa0, 0x20, 0x8a, 0x20, 0x28,
+    0x88, 0x82, 0x08, 0x02, 0xa2, 0x28, 0x02, 0x82,
+    0x08, 0x8a, 0xa0, 0x20, 0x88, 0x22, 0x28, 0x88,
+    0xa2, 0x00, 0x02];
+
+pub fn is_pw210_candidate(num: u128) -> bool {
+    if num < 11 {
+        match num {
+            2 | 3 | 5 | 7 => true,
+            _ => false,
+        }
+    } else {
+        let index = (num % 210) as usize; // Calculate bit position (0 to 209)
+        let byte_index = index / 8; // Calculate byte index within the array
+        let bit_mask = 1 << (index % 8); // Calculate bit-mask within the byte
+        PW210_BITMAP[byte_index] & bit_mask > 0
+    }
+}

src/lib.rs

@@ -5,7 +5,7 @@ use std::cmp::{min, Ordering};
 use std::convert::From;
 use std::fmt;
 use genawaiter::stack::let_gen_using;
-use candidates::prime_wheel_30;
+use candidates::{prime_wheel_210, is_pw210_candidate};
 use num::integer::Roots;
 
 #[derive(Clone, Copy, Debug, Eq, Ord, PartialEq, PartialOrd)]
@@ -96,7 +96,7 @@ impl From<u128> for PrimeFactors {
         if n < 2 { return pf; }
         // A factor of n must have a value less than or equal to sqrt(n)
         let mut maxf = n.sqrt() + 1;
-        let_gen_using!(mpgen, prime_wheel_30);
+        let_gen_using!(mpgen, prime_wheel_210);
         let mut x = n;
         for f in mpgen.into_iter() {
             if f >= maxf {
@@ -145,17 +145,12 @@ impl<'a> IntoIterator for &'a PrimeFactors {
 
 /// Test if the value is a prime number, or not
 pub fn u128_is_prime(n: u128) -> bool {
-    if n < 2 { return false; }
-    if n > 30 {
-        // check spoke in the prime wheel (base 30)
-        match n % 30 {
-            1|7|11|13|17|19|23|29 => (), // may be prime
-            _ => return false, // cannot be prime
-        }
+    if !is_pw210_candidate(n) {
+        return false;
     }
     // A factor of n must have a value less than or equal to sqrt(n)
     let maxf = n.sqrt() + 1;
-    let_gen_using!(mpgen, prime_wheel_30);
+    let_gen_using!(mpgen, prime_wheel_210);
     for f in mpgen.into_iter() {
         if f >= maxf {
             break;

tests/test.rs

@@ -5,14 +5,15 @@ use reikna;
 use genawaiter::stack::let_gen_using;
 use primefactor::{
     candidates::prime_wheel_30,
+    candidates::{prime_wheel_210, is_pw210_candidate},
     primefactor_gcd,
     PrimeFactors,
     u128_gcd,
     u128_is_prime,
     u128_lcm};
 
 #[test]
-fn test_early_prime_wheel_numbers() {
+fn test_early_prime_wheel_30_numbers() {
     let testvec = vec![
         2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49, 53, 59,
         61, 67, 71, 73, 77, 79, 83, 89, 91, 97, 101, 103, 107, 109, 113
@@ -26,21 +27,52 @@ fn test_early_prime_wheel_numbers() {
 }
 
 #[test]
-fn test_early_prime_candidate_numbers() {
-    let testvec = vec![
-        2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49, 53, 59,
-        61, 67, 71, 73, 77, 79, 83, 89, 91, 97, 101, 103, 107, 109, 113
-    ];
-    let_gen_using!(mpgen, prime_wheel_30);
+fn test_1000th_prime_with_pw30() {
+    let mut primes = 0;
+    let mut prime = 0;
+    for i in 0..7920 {
+        if u128_is_prime(i as u128) {
+            primes += 1;
+            prime = i;
+        }
+    }
+    assert_eq!(primes, 1000);
+    assert_eq!(prime, 7919);
+}
+
+#[test]
+fn test_early_prime_wheel_210_numbers() {
+    let testvec: Vec<u128> = vec![
+        2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61,
+        67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 121, 127, 131,
+        137, 139, 143, 149, 151, 157, 163, 167, 169, 173, 179, 181, 187, 191,
+        193, 197, 199, 209, 211, 221, 223, 227, 229, 233, 239, 241, 247, 251,
+        253, 257, 263, 269, 271, 277, 281, 283, 289, 293, 299, 307, 311, 313];
+    let_gen_using!(mpgen, prime_wheel_210);
     let mut mp = mpgen.into_iter();
     for i in 0..testvec.len() {
         let p = mp.next().unwrap();
         assert_eq!(testvec[i], p);
+        assert!(is_pw210_candidate(p));
+    }
+}
+
+#[test]
+fn test_1000th_prime_with_pw210() {
+    let mut primes = 0;
+    let mut prime = 0;
+    for i in 0..7920 {
+        if u128_is_prime(i as u128) {
+            primes += 1;
+            prime = i;
+        }
     }
+    assert_eq!(primes, 1000);
+    assert_eq!(prime, 7919);
 }
 
 #[test]
-fn test_prime_wheel_quality() {
+fn test_prime_wheel_30_quality() {
     let mut primes: u128 = 0;
     let mut others: u128 = 0;
     let_gen_using!(mpgen, prime_wheel_30);
@@ -59,6 +91,26 @@ fn test_prime_wheel_quality() {
     assert!(percent > 25.0);
 }
 
+#[test]
+fn test_prime_wheel_210_quality() {
+    let mut primes: u128 = 0;
+    let mut others: u128 = 0;
+    let_gen_using!(mpgen, prime_wheel_210);
+    let mut mp = mpgen.into_iter();
+    for _ in 0..1000000 {
+        let p = mp.next().unwrap();
+        if u128_is_prime(p) {
+            primes += 1;
+        } else {
+            others += 1;
+        }
+    }
+    let percent = primes as f64 / (primes + others) as f64 * 100.0;
+    println!("Prime wheel generated {}/{} ({:.3}%) primes",
+             primes, primes+others, percent);
+    assert!(percent > 30.0);
+}
+
 #[test]
 fn test_is_prime() {
     for num in 2..=1000 {