prime-formula 0.3.1

High-performance prime number generation and constellation finding using novel wheel factorization
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117

//! Library for generating likely primes in the range 10^9 to 10^38
//!
//! The algorithm uses The Periodic Table of Primes to select candidates
//! for prime numbers within the given range, and then in parallel runs
//! a u128 optimized version of the Miller-Rabin tests to find likely primes
//!
//! For primes below 10^18 then a deterministic algorithm will be used
//! which greatly increases the likelihood of finding true primes. Above
//! this limit, the probabilistic variant of the algorithm will be used
//! which may return pseudoprimes unless the number of iterations (k) is
//! set sufficiently high.
//
// Copyright: Adam Cottrell (cottrela@gmail.com)

use rayon::prelude::*;

use crate::common::{PRIME_ROOTS, align_to_cycle};
pub mod primality;

/// Generate huge prime numbers using a primality test
///
/// Using primality tests is quicker for very large prime
/// numbers, but considerably longer than a sieve for smaller
/// prime numbers.

#[derive(Debug, PartialEq)]
pub enum PrimalityTest {
    MillerRabinDeterministic,
    MillerRabinProbabilistic,
}

/// Use the formula of primes to find large primes
pub struct PrimeFormula {
    k: u32,              // Number of tests to use (limited to 48 in deterministic mode)
    mode: PrimalityTest, // Primality test
}

impl PrimeFormula {
    pub fn new(mode: PrimalityTest, k: u32) -> Self {
        PrimeFormula { k, mode }
    }

    /// Get primes within the start/end range (inclusive)
    ///
    /// Note that primes below 11 will not be found using this
    /// formula, and so must be manually added by the caller if
    /// needed.
    ///
    /// This function returns the primes in order of their roots,
    /// therefore may need sorting before use. Set get_primes.
    pub fn get_unsorted_primes(&self, start: u128, end: u128) -> Vec<u128> {
        let (b_start, b_end) = align_to_cycle(start, end);
        let primes: Vec<u128> = PRIME_ROOTS
            .into_par_iter()
            .flat_map(|root| {
                (b_start..b_end)
                    .map(|i| root as u128 + i * 210)
                    .filter(|p| *p >= start && *p <= end) // Filter by range
                    .filter(|p| {
                        primality::is_prime(
                            p,
                            self.k,
                            self.mode == PrimalityTest::MillerRabinDeterministic,
                        )
                    })
                    .collect::<Vec<u128>>() // Collect inner results
            })
            .collect(); // Collect parallel results

        primes
    }

    /// Get primes within the start/end range (inclusive)
    ///
    /// Note that primes below 11 will not be found using this
    /// formula, and so must be manually added by the caller if
    /// needed.
    ///
    /// This function returns the primes sorted. To get the unsorted
    /// primes, which is faster, use get_unsorted_primes.
    pub fn get_primes(&self, start: u128, end: u128) -> Vec<u128> {
        let mut primes = self.get_unsorted_primes(start, end);

        primes.par_sort_unstable();
        primes
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_deterministic_primes() {
        let primes;
        let inst = PrimeFormula::new(PrimalityTest::MillerRabinDeterministic, 2);

        // Check the first million primes
        primes = inst.get_primes(11, 1_000_000);
        assert_eq!(primes.len(), 78498 - 4); // Number of primes below 1e6, excluding 2,3,5,7
        assert_eq!(primes[0], 11); // First prime from root 11
        assert_eq!(primes[78498 - 4 - 1], 999_983); // Last prime before 1e6
    }

    #[test]
    fn test_probabilistic_primes() {
        // Very low chance of failure for low order primes
        let primes;
        let inst = PrimeFormula::new(PrimalityTest::MillerRabinProbabilistic, 4);

        // Check the first million primes
        primes = inst.get_primes(11, 1_000);
        assert_eq!(primes.len(), 168 - 4); // Number of primes below 1e6, excluding 2,3,5,7
        assert_eq!(primes[0], 11); // First prime from root 11
        assert_eq!(primes[168 - 4 - 1], 997); // Last prime before 1e6
    }
}