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

//! Fast primality tests for big integers using Miller-Rabin tests
//!
//! Algorithm is optimized for deep checking one prime number, and attempts to
//! run all of the witness checks in parallel to reject composites as quickly
//! as possible.
//
// Copyright: Adam Cottrell (cottrela@gmail.com)

use num_bigint::{BigUint, RandBigInt};
use num_traits::{FromPrimitive, One, Zero};
use rayon::prelude::*;
use std::sync::Arc;

fn decompose(n: &BigUint) -> (BigUint, u128) {
    let mut d = n - BigUint::one();
    let mut s = 0;

    // While d is even
    while &d & BigUint::one() == BigUint::zero() {
        d >>= 1;
        s += 1;
    }
    (d, s)
}

fn is_witness(a: BigUint, d: &BigUint, s: u128, n: &BigUint) -> bool {
    let mut x = a.modpow(d, n);
    let two: BigUint = BigUint::from_u8(2).unwrap();

    if x == BigUint::one() || x == n - BigUint::one() {
        return false;
    }
    for _ in 0..s - 1 {
        x = x.modpow(&two, n);
        if x == n - BigUint::one() {
            return false;
        }
    }
    true
}

/// Non-deterministic primality test
///
/// Performs Miller-Rabin algorithm, but with a random witness in the range
/// 2..n-2 where n is the candidate number.
///
/// Note that this function is optimized to work on groups of large primes
/// but will by design not check primes less than 11 (which is the first root
/// of prime)
pub fn is_likely_big_prime(candidate: &BigUint, num_checks: u32) -> bool {
    let (d, s) = decompose(candidate);
    let candidate_arc = Arc::new(candidate.clone());
    let d_arc = Arc::new(d);
    let two: BigUint = BigUint::from_u8(2).unwrap();

    // Use parallel any to short-circuit on first composite proof
    !(0..num_checks).into_par_iter().any(|_| {
        let mut rng = rand::thread_rng();
        let candidate = Arc::clone(&candidate_arc);
        let d = Arc::clone(&d_arc);

        // Generate random witness in [2, n-2]
        let a = rng.gen_biguint_range(&two, &(candidate.as_ref() - &two));

        is_witness(a, d.as_ref(), s, candidate.as_ref())
    })
}

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

    #[test]
    fn test_decompose() {
        assert_eq!(decompose(&BigUint::from(13u8)), (BigUint::from(3u8), 2)); // 12 = 3 * 2^2
        assert_eq!(decompose(&BigUint::from(9u8)), (BigUint::from(1u8), 3)); // 8 = 1 * 2^3
        assert_eq!(decompose(&BigUint::from(17u8)), (BigUint::from(1u8), 4)); // 16 = 1 * 2^4
        assert_eq!(decompose(&BigUint::from(7u8)), (BigUint::from(3u8), 1)); // 6 = 3 * 2^1
    }

    #[test]
    fn test_probabilistic_primes() {
        // Test with 10 iterations for reliability
        assert!(is_likely_big_prime(&BigUint::from(7919u32), 10)); // Prime
        assert!(!is_likely_big_prime(&BigUint::from(7921u32), 10)); // 89^2 = 7921

        // Large prime
        assert!(is_likely_big_prime(
            &BigUint::from(0xFFFFFFFFFFFFFFC5u128),
            10
        ));
    }

    #[test]
    fn test_large_numbers() {
        // Mersenne prime (2^61 - 1)
        let mersenne61 = BigUint::from((1u128 << 61) - 1);
        assert!(is_likely_big_prime(&mersenne61, 5));

        // Large composite (m61 ^ 2)
        let large_composite = &mersenne61 * &mersenne61;
        assert!(!is_likely_big_prime(&large_composite, 5));
    }
}