prime-formula 0.3.1

//! Common library for generating prime numbers
//
// Copyright: Adam Cottrell (cottrela@gmail.com)

/// Base primes excluded from root calculations (2, 3, 5, 7)
///
/// These form the foundation of the prime period wheel factorization
pub const NON_ROOT_PRIMES: [u8; 4] = [2, 3, 5, 7];

/// The prime period (210) - product of base primes
///
/// All prime roots >7 will be coprime with this value
pub const PRIME_PERIOD: u8 = {
    let mut product = 1;
    let mut i = 0;
    while i < NON_ROOT_PRIMES.len() {
        product *= NON_ROOT_PRIMES[i];
        i += 1;
    }
    product
};

/// Prime roots coprime with the base primes (2,3,5,7)
///
/// These 48 roots form the basis for all prime number generation in this sieve.
/// Any prime number >7 can be expressed as: root + k×210
pub const PRIME_ROOTS: [u8; 48] = {
    let mut primes = [0u8; 48];
    let mut num = 2;
    let mut idx = 0;

    // Iterate through all candidates in range
    while num < PRIME_PERIOD + 2 {
        // Check if number is coprime with base primes
        if num % 2 != 0 && num % 3 != 0 && num % 5 != 0 && num % 7 != 0 {
            primes[idx] = num;
            idx += 1;
        }
        num += 1;
    }

    // Compile-time validation of array size
    if idx != 48 {
        panic!("Prime roots count mismatch");
    }

    primes
};

/// The number of prime roots (48)
pub const NUM_ROOTS: usize = PRIME_ROOTS.len();

/// Public function to align start and end on 210 boundaries
pub fn align_to_cycle(start: u128, end: u128) -> (u128, u128) {
    // Find w/c aligned range for a given start/end window
    //
    // For any prime number (P):
    //    P = r + b * 210
    //    b = (P - r) / 210
    //
    // Where:
    //   r is the prime root (one of 48)
    //   b is the multiplier (1..N)
    let r_min = PRIME_ROOTS[0] as u128; // 13
    let r_max = PRIME_ROOTS[NUM_ROOTS - 1] as u128; // 211

    // Get ceil(b_start) and floor(b_end) to find the tightest range
    let b_start = (start.saturating_sub(r_max)).div_ceil(210); // Use integer ceil
    let b_end = 1 + (end - r_min) / 210; // Use integer floor

    // If b_start == b_end, then there is nothing to do
    (b_start, b_end)
}

/// Helper function to find the inverse prime root values.
/// The inverse of a prime root (q) is the congruent
/// of q * q_hat - r = 0, solved for q_hat.
/// Where r is a different prime root.
///
/// n.b. there is one solution per prime root that satisfies
/// this equation.
pub fn find_inv_prime_root_vec(root_idx: usize) -> Vec<u8> {
    let r = PRIME_ROOTS[root_idx] as u32;
    PRIME_ROOTS
        .iter()
        .copied()
        .map(|q| {
            let period = PRIME_PERIOD as u32;
            let q = q as u32;
            PRIME_ROOTS
                .iter()
                .copied()
                .find(|&h_candidate| (q * h_candidate as u32) % period == r % period)
                .expect("No inverse prime root found")
        })
        .collect()
}

/// Function to generate cyclic table of composites (CTC)
/// These offsets are used as look-ups for each
/// prime root.
///
/// Constructing the CTC is the first step in creating
/// the Periodic Table of Primes (PTP).
pub fn get_cyclic_composite_vec(root_index: usize, inv_prime_root_table_row: &[u8]) -> Vec<u8> {
    let r = PRIME_ROOTS[root_index] as u32;
    let period = PRIME_PERIOD as u32;
    PRIME_ROOTS
        .iter()
        .enumerate()
        .map(|(j, &q)| {
            let q_hat = inv_prime_root_table_row[j] as u32;
            let q = q as u32;
            let val = if r > (q * q_hat) {
                1 + (period + q * q_hat - r) / period
            } else {
                1 + (q * q_hat - r) / period
            };
            val as u8
        })
        .collect()
}

/// Get primes 2,3,5,7 within range [start, end]
pub fn get_small_primes(start: u128, end: u128) -> Vec<u128> {
    [2, 3, 5, 7]
        .iter()
        .filter(|&&p| p >= start && p <= end)
        .cloned()
        .collect()
}

/// High-level mode for the selecting algorithm to use
#[derive(Debug, Copy, Clone, PartialEq)]
pub enum ConfigMode {
    /// Novel sieve method - best for medium ranges
    Sieve,
    /// Windowed search - optimized for small number ranges
    Window,
    /// Use Miller-Rabin test with deterministic witnesses
    MillerRabinDeterministic,
    /// Use Miller-Rabin test with random witnesses
    MillerRabinProbabilistic,
}

/// Helper function to select ConfigMode based on the
/// start/end range from big-O approximations. The goal
/// being to select the fastest algorithm for the range.
pub fn get_auto_mode(_start: u128, end: u128) -> ConfigMode {
    let big_o_sieve = end.isqrt(); // Dominated by sqrt(n)
    let big_o_miller = end.ilog10().pow(6) as u128;

    if big_o_sieve > big_o_miller {
        // For very large numbers, use best-effort
        // to reduce the number of options - this
        // may return pseudoprimes.
        if end > u64::MAX as u128 {
            ConfigMode::MillerRabinProbabilistic
        } else {
            ConfigMode::MillerRabinDeterministic
        }
    } else {
        // For everything else, sieve mode should give
        // the best overall performance. It is exhaustive
        // and will quickly eliminate a lot of composites.
        ConfigMode::Sieve
    }
}

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

    #[test]
    fn test_constants() {
        // Verify base primes
        assert_eq!(NON_ROOT_PRIMES, [2, 3, 5, 7]);

        // Verify prime period calculation
        assert_eq!(PRIME_PERIOD, 2 * 3 * 5 * 7);

        // Validate prime roots count and properties
        assert_eq!(PRIME_ROOTS.len(), 48);
        for &root in &PRIME_ROOTS {
            assert!(root >= 11);
            assert!(root <= 211);
            assert_ne!(root % 2, 0);
            assert_ne!(root % 3, 0);
            assert_ne!(root % 5, 0);
            assert_ne!(root % 7, 0);
        }
    }

    #[test]
    fn test_align_to_cycle() {
        // Test exact boundaries
        assert_eq!(align_to_cycle(0, 210), (0, 1));
        assert_eq!(align_to_cycle(210, 420), (0, 2)); // 210 < 211 so should include first cycle as well as second

        // Test values within cycle
        assert_eq!(align_to_cycle(100, 200), (0, 1));
        assert_eq!(align_to_cycle(300, 500), (1, 3));

        // Test edge cases
        assert_eq!(align_to_cycle(211, 211), (0, 1)); // Upper boundary
        assert_eq!(align_to_cycle(13, 13), (0, 1)); // Lower boundary
    }

    #[test]
    fn test_inverse_prime_roots() {
        // Test first and last roots
        let first_inverses = find_inv_prime_root_vec(0);
        assert_eq!(first_inverses.len(), 48);

        let last_inverses = find_inv_prime_root_vec(47);
        assert_eq!(last_inverses.len(), 48);

        // Verify inverse property for sample root
        let test_root = PRIME_ROOTS[10]; // Arbitrary index
        let inverses = find_inv_prime_root_vec(10);
        for (i, &h) in inverses.iter().enumerate() {
            let q = PRIME_ROOTS[i];
            assert_eq!(
                (q as u32 * h as u32) % PRIME_PERIOD as u32,
                test_root as u32 % PRIME_PERIOD as u32
            );
        }
    }

    #[test]
    fn test_cyclic_composites() {
        // Test known values
        let test_index = PRIME_ROOTS.iter().position(|&r| r == 11).unwrap();
        let inverses = find_inv_prime_root_vec(test_index);
        let composites = get_cyclic_composite_vec(test_index, &inverses);

        // Should generate predictable offsets
        assert_eq!(composites.len(), 48);
        assert!(composites.iter().all(|&v| v > 0));
    }

    #[test]
    fn test_small_primes() {
        // Full range
        assert_eq!(get_small_primes(2, 10), vec![2, 3, 5, 7]);

        // Partial range
        assert_eq!(get_small_primes(5, 7), vec![5, 7]);

        // Empty range
        assert!(get_small_primes(8, 10).is_empty());
    }

    #[test]
    fn test_auto_mode_selection() {
        // Small numbers should use sieve
        assert_eq!(get_auto_mode(0, 1_000), ConfigMode::Sieve);

        // Very large numbers should use probabilistic MR
        assert_eq!(
            get_auto_mode(0, u128::MAX),
            ConfigMode::MillerRabinProbabilistic
        );

        // Borderline case (sqrt(n) vs log^6(n))
        let borderline = 10u128.pow(15);
        let mode = get_auto_mode(0, borderline);
        assert!(matches!(
            mode,
            ConfigMode::Sieve | ConfigMode::MillerRabinDeterministic
        ));
    }

    #[test]
    fn test_prime_root_validity() {
        // Verify all roots are in correct range
        assert!(PRIME_ROOTS.iter().all(|&r| (11..=211).contains(&r)));

        // Verify no duplicates
        let mut sorted = PRIME_ROOTS.to_vec();
        sorted.sort_unstable();
        sorted.dedup();
        assert_eq!(sorted.len(), 48);
    }

    #[test]
    fn test_boundary_alignment() {
        // Test alignment near cycle boundaries
        assert_eq!(align_to_cycle(209, 211), (0, 1));
        assert_eq!(align_to_cycle(211, 421), (0, 2));

        // Test with extremely large numbers
        let big_num = 210u128.pow(10);
        assert_eq!(
            align_to_cycle(big_num, big_num + 100),
            (big_num / 210 - 1, big_num / 210 + 1)
        );
    }
}