tx2_iff/

prime.rs

//! 360 Prime Pattern and PPF-based quantization
//!
//! This module implements the 360 Prime Pattern for optimal quantization
//! of wavelet coefficients. The pattern is based on the mathematical discovery
//! that every prime can be located within distance ≤180 from either:
//! 1. A factor of m×360, OR
//! 2. A term in the recursive sequence N_i
//!
//! ## Mathematical Foundation
//!
//! The number 360 = 2³ × 3² × 5 has special properties:
//! - φ(360) = 96 residue classes coprime to 360
//! - 48 residue classes are always factor-covered
//! - 46 residue classes are always sequence-covered
//! - 2 residue classes (261, 269) are scale-dependent
//!
//! ## Convergence Properties
//!
//! The 360-prime pattern exhibits remarkable convergence:
//! - Precision gain: 2.56 decimal digits per iteration
//! - d_n ≈ -log₁₀(K) + n × log₁₀(360)
//! - After 22 iterations: 64 correct digits of π

use serde::{Deserialize, Serialize};

/// Number of residue classes modulo 360 that can contain primes
pub const RESIDUE_CLASSES: usize = 96;

/// The 96 residue classes modulo 360 that are coprime to 360
/// These are the only residue classes that can contain primes > 5
pub const PRIME_RESIDUE_CLASSES: [u16; RESIDUE_CLASSES] = [
    1, 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, 119,
    121, 127, 131, 133, 137, 139, 143, 149, 151, 157, 161, 163, 167, 169, 173, 179,
    181, 187, 191, 193, 197, 199, 203, 209, 211, 217, 221, 223, 227, 229, 233, 239,
    241, 247, 251, 253, 257, 259, 263, 269, 271, 277, 281, 283, 287, 289, 293, 299,
    301, 307, 311, 313, 317, 319, 323, 329, 331, 337, 341, 343, 347, 349, 353, 359,
];

/// Category 1: Residue classes that are always factor-covered (48 total)
pub const FACTOR_COVERED: [u16; 48] = [
    1, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49, 53, 59,
    61, 67, 71, 73, 77, 79, 83, 89, 271, 277, 281, 283, 287, 289, 293, 299,
    301, 307, 311, 313, 317, 319, 323, 329, 331, 337, 341, 343, 347, 349, 353, 359,
];

/// Category 2: Residue classes that are always sequence-covered (46 total)
pub const SEQUENCE_COVERED: [u16; 46] = [
    91, 97, 101, 103, 107, 109, 113, 119, 121, 127, 131, 133, 137, 139, 143, 149,
    151, 157, 161, 163, 167, 169, 173, 179, 181, 187, 191, 193, 197, 199, 203, 209,
    211, 217, 221, 223, 227, 229, 233, 239, 241, 247, 251, 253, 257, 259,
];

/// Category 3: Residue classes with scale-dependent coverage (2 total)
pub const SCALE_DEPENDENT: [u16; 2] = [261, 269];

/// Extended prime set including -1 as the Sign Prime
#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
pub enum ExtendedPrime {
    /// The Sign Prime: -1 (unique negative prime)
    SignPrime,
    /// Magnitude primes: 2, 3, 5, 7, 11, ...
    MagnitudePrime(u32),
}

impl ExtendedPrime {
    /// Check if a number is prime in the PPF framework
    pub fn is_prime(n: i32) -> bool {
        if n == -1 {
            return true; // Sign prime
        }
        if n < 2 {
            return false;
        }
        if n == 2 || n == 3 {
            return true;
        }
        if n % 2 == 0 || n % 3 == 0 {
            return false;
        }

        let mut i = 5;
        while i * i <= n {
            if n % i == 0 || n % (i + 2) == 0 {
                return false;
            }
            i += 6;
        }
        true
    }

    /// Get the value as i32
    pub fn value(&self) -> i32 {
        match self {
            ExtendedPrime::SignPrime => -1,
            ExtendedPrime::MagnitudePrime(p) => *p as i32,
        }
    }
}

/// PPF quantization table based on 360-prime pattern
#[derive(Debug, Clone)]
pub struct QuantizationTable {
    /// Base quantization value
    pub base_q: u16,
    /// Per-residue class weights (96 entries)
    pub weights: [f32; RESIDUE_CLASSES],
}

impl QuantizationTable {
    /// Create a new quantization table with given base Q
    pub fn new(base_q: u16) -> Self {
        let mut weights = [1.0f32; RESIDUE_CLASSES];

        // Apply weighting based on category
        // Factor-covered classes get slightly higher weight (sharper quantization)
        for &residue in &FACTOR_COVERED {
            if let Some(idx) = PRIME_RESIDUE_CLASSES.iter().position(|&r| r == residue) {
                weights[idx] = 0.9; // 10% finer quantization
            }
        }

        // Sequence-covered classes get standard weight
        for &residue in &SEQUENCE_COVERED {
            if let Some(idx) = PRIME_RESIDUE_CLASSES.iter().position(|&r| r == residue) {
                weights[idx] = 1.0;
            }
        }

        // Scale-dependent classes get adaptive weight
        for &residue in &SCALE_DEPENDENT {
            if let Some(idx) = PRIME_RESIDUE_CLASSES.iter().position(|&r| r == residue) {
                weights[idx] = 1.1; // 10% coarser quantization
            }
        }

        QuantizationTable { base_q, weights }
    }

    /// Get quantization step for a given position
    pub fn get_step(&self, x: usize, y: usize) -> u16 {
        // Map position to residue class using 360-pattern
        let residue = ((x + y * 360) % 360) as u16;

        // Find the closest prime residue class
        let idx = self.find_closest_residue(residue);

        // Apply weight
        let weighted_q = self.base_q as f32 * self.weights[idx];
        weighted_q.round() as u16
    }

    /// Find the index of the closest prime residue class
    fn find_closest_residue(&self, residue: u16) -> usize {
        let mut min_dist = u16::MAX;
        let mut min_idx = 0;

        for (idx, &prime_residue) in PRIME_RESIDUE_CLASSES.iter().enumerate() {
            let dist = if residue > prime_residue {
                residue - prime_residue
            } else {
                prime_residue - residue
            };

            // Also consider wraparound distance
            let wraparound_dist = 360 - dist;
            let effective_dist = dist.min(wraparound_dist);

            if effective_dist < min_dist {
                min_dist = effective_dist;
                min_idx = idx;
            }
        }

        min_idx
    }
}

/// PPF-based hash function using 360-prime pattern
pub struct PpfHash {
    /// Seed value
    seed: u32,
}

impl PpfHash {
    /// Create a new PPF hash with given seed
    pub fn new(seed: u32) -> Self {
        PpfHash { seed }
    }

    /// Hash a coordinate pair to a value in [0, 1)
    pub fn hash(&self, x: u32, y: u32) -> f32 {
        // Combine coordinates with seed using prime factorization properties
        let combined = self.combine_coords(x, y);

        // Apply Fibonacci hashing (golden ratio constant)
        let hash = combined.wrapping_mul(2654435761u32);

        // Map to 360-prime residue
        let residue = (hash % 360) as usize;

        // Find closest prime residue class
        let idx = self.find_prime_residue(residue);

        // Normalize using the prime residue value
        PRIME_RESIDUE_CLASSES[idx] as f32 / 360.0
    }

    /// Combine coordinates using PPF properties
    fn combine_coords(&self, x: u32, y: u32) -> u32 {
        // Use the fact that (-1) × (-1) = (+1) in PPF
        // Combine x and y with seed through XOR (preserves prime structure)
        let x_prime = self.prime_factorize_approx(x);
        let y_prime = self.prime_factorize_approx(y);
        x_prime ^ y_prime ^ self.seed
    }

    /// Approximate prime factorization hash
    /// Maps input to a value that preserves residue class properties
    fn prime_factorize_approx(&self, n: u32) -> u32 {
        if n == 0 {
            return 0;
        }

        // Extract factors of 2, 3, 5 (divisors of 360)
        let mut residue = n;
        let mut factor_hash = 1u32;

        while residue % 2 == 0 {
            factor_hash = factor_hash.wrapping_mul(2);
            residue /= 2;
        }
        while residue % 3 == 0 {
            factor_hash = factor_hash.wrapping_mul(3);
            residue /= 3;
        }
        while residue % 5 == 0 {
            factor_hash = factor_hash.wrapping_mul(5);
            residue /= 5;
        }

        // Combine with residue
        factor_hash.wrapping_add(residue)
    }

    /// Find the index of the closest prime residue class
    fn find_prime_residue(&self, residue: usize) -> usize {
        let residue = residue as u16;
        let mut min_dist = u16::MAX;
        let mut min_idx = 0;

        for (idx, &prime_residue) in PRIME_RESIDUE_CLASSES.iter().enumerate() {
            let dist = if residue > prime_residue {
                residue - prime_residue
            } else {
                prime_residue - residue
            };

            let wraparound_dist = 360 - dist;
            let effective_dist = dist.min(wraparound_dist);

            if effective_dist < min_dist {
                min_dist = effective_dist;
                min_idx = idx;
            }
        }

        min_idx
    }
}

/// Recursive sequence for 360-prime pattern
pub struct RecursiveSequence {
    /// Current scale (m)
    scale: u32,
    /// Current index in sequence
    index: u32,
    /// Current value
    value: u32,
}

impl RecursiveSequence {
    /// Create a new recursive sequence for given scale
    pub fn new(scale: u32) -> Self {
        // N_1 = (m-1) × 360 + 181
        let value = (scale - 1) * 360 + 181;
        RecursiveSequence {
            scale,
            index: 1,
            value,
        }
    }

    /// Get the current value
    pub fn value(&self) -> u32 {
        self.value
    }

    /// Advance to next term: N_i = N_(i-1) + i
    pub fn next(&mut self) -> u32 {
        self.index += 1;
        self.value += self.index;
        self.value
    }

    /// Check if sequence has exceeded the scale range
    pub fn in_range(&self) -> bool {
        self.value <= self.scale * 360 + 180
    }
}

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

    #[test]
    fn test_residue_classes() {
        // φ(360) = 96
        assert_eq!(PRIME_RESIDUE_CLASSES.len(), RESIDUE_CLASSES);

        // Verify all are coprime to 360
        for &residue in &PRIME_RESIDUE_CLASSES {
            assert_eq!(gcd(residue as u32, 360), 1);
        }
    }

    #[test]
    fn test_category_distribution() {
        // 48 + 46 + 2 = 96
        assert_eq!(
            FACTOR_COVERED.len() + SEQUENCE_COVERED.len() + SCALE_DEPENDENT.len(),
            RESIDUE_CLASSES
        );
    }

    #[test]
    fn test_extended_prime() {
        assert!(ExtendedPrime::is_prime(-1)); // Sign prime
        assert!(ExtendedPrime::is_prime(2));
        assert!(ExtendedPrime::is_prime(3));
        assert!(ExtendedPrime::is_prime(5));
        assert!(ExtendedPrime::is_prime(7));

        assert!(!ExtendedPrime::is_prime(0));
        assert!(!ExtendedPrime::is_prime(1));
        assert!(!ExtendedPrime::is_prime(4));
        assert!(!ExtendedPrime::is_prime(6));

        // -2, -3, -5 are not prime (can be factored as (-1) × 2, etc.)
        assert!(!ExtendedPrime::is_prime(-2));
        assert!(!ExtendedPrime::is_prime(-3));
    }

    #[test]
    fn test_quantization_table() {
        let table = QuantizationTable::new(10);

        // Base Q should be 10
        assert_eq!(table.base_q, 10);

        // Should have 96 weights
        assert_eq!(table.weights.len(), RESIDUE_CLASSES);

        // Get step for various positions
        let step1 = table.get_step(0, 0);
        let step2 = table.get_step(100, 100);

        // Steps should be close to base_q with weighting
        assert!(step1 >= 8 && step1 <= 12);
        assert!(step2 >= 8 && step2 <= 12);
    }

    #[test]
    fn test_ppf_hash_determinism() {
        let hash = PpfHash::new(42);

        // Same input should give same output
        let h1 = hash.hash(10, 20);
        let h2 = hash.hash(10, 20);
        assert_eq!(h1, h2);

        // Different inputs should give different outputs
        let h3 = hash.hash(10, 21);
        assert_ne!(h1, h3);

        // Hash should be in [0, 1)
        assert!(h1 >= 0.0 && h1 < 1.0);
    }

    #[test]
    fn test_recursive_sequence() {
        // Test scale m=1
        let mut seq = RecursiveSequence::new(1);
        assert_eq!(seq.value(), 181); // N_1 = (1-1)×360 + 181

        assert_eq!(seq.next(), 183); // N_2 = 181 + 2
        assert_eq!(seq.next(), 186); // N_3 = 183 + 3
        assert_eq!(seq.next(), 190); // N_4 = 186 + 4

        // Test scale m=2
        let mut seq2 = RecursiveSequence::new(2);
        assert_eq!(seq2.value(), 541); // N_1 = (2-1)×360 + 181 = 541
    }

    #[test]
    fn test_sequence_coverage() {
        // Verify sequence stays in range for scale m=1
        let mut seq = RecursiveSequence::new(1);
        let mut count = 0;

        while seq.in_range() && count < 1000 {
            seq.next();
            count += 1;
        }

        // Should have generated a reasonable number of terms
        assert!(count > 10);
        assert!(count < 1000); // Should terminate
    }

    // Helper function for GCD
    fn gcd(mut a: u32, mut b: u32) -> u32 {
        while b != 0 {
            let temp = b;
            b = a % b;
            a = temp;
        }
        a
    }
}