quantrs2-tytan 0.2.0

//! Variational Quantum Factoring (VQF) implementation.
//!
//! This module provides quantum algorithms for integer factorization
//! using variational approaches.

#![allow(dead_code)]

use crate::hybrid_algorithms::{AnsatzType, ClassicalOptimizer, Hamiltonian, PauliTerm, VQE};
use scirs2_core::random::prelude::*;
use scirs2_core::random::prelude::*;

/// Variational Quantum Factoring solver
pub struct VQF {
    /// Number to factor
    n: u64,
    /// Number of qubits for factors
    n_qubits_p: usize,
    n_qubits_q: usize,
    /// Optimization method
    optimizer: ClassicalOptimizer,
    /// Penalty strength for constraints
    penalty_strength: f64,
    /// Use preprocessing
    use_preprocessing: bool,
    /// Symmetry reduction
    use_symmetry_reduction: bool,
}

impl VQF {
    /// Create new VQF solver
    pub fn new(n: u64, optimizer: ClassicalOptimizer) -> Result<Self, String> {
        if n < 2 {
            return Err("Number must be >= 2".to_string());
        }

        // Estimate number of qubits needed
        let n_bits = 64 - n.leading_zeros() as usize;
        let n_qubits_p = n_bits.div_ceil(2);
        let n_qubits_q = n_bits - n_qubits_p;

        Ok(Self {
            n,
            n_qubits_p,
            n_qubits_q,
            optimizer,
            penalty_strength: 100.0,
            use_preprocessing: true,
            use_symmetry_reduction: true,
        })
    }

    /// Set penalty strength
    pub const fn with_penalty_strength(mut self, strength: f64) -> Self {
        self.penalty_strength = strength;
        self
    }

    /// Enable/disable preprocessing
    pub const fn with_preprocessing(mut self, use_preprocessing: bool) -> Self {
        self.use_preprocessing = use_preprocessing;
        self
    }

    /// Enable/disable symmetry reduction
    pub const fn with_symmetry_reduction(mut self, use_symmetry: bool) -> Self {
        self.use_symmetry_reduction = use_symmetry;
        self
    }

    /// Factor the number
    pub fn factor(&mut self) -> Result<FactorizationResult, String> {
        // Preprocessing
        if self.use_preprocessing {
            if let Some(factors) = self.preprocess() {
                return Ok(factors);
            }
        }

        // Build factorization Hamiltonian
        let hamiltonian = self.build_hamiltonian()?;

        // Create VQE solver
        let total_qubits = self.n_qubits_p + self.n_qubits_q;
        let mut vqe = VQE::new(
            total_qubits,
            AnsatzType::HardwareEfficient {
                depth: 3,
                entangling_gate: "CZ".to_string(),
            },
            self.optimizer.clone(),
        );

        // Run VQE
        let vqe_result = vqe.solve(&hamiltonian)?;

        // Extract factors from solution
        self.extract_factors(&vqe_result.optimal_parameters)
    }

    /// Preprocess for trivial cases
    fn preprocess(&self) -> Option<FactorizationResult> {
        // Check if even
        if self.n % 2 == 0 {
            return Some(FactorizationResult {
                p: 2,
                q: self.n / 2,
                confidence: 1.0,
                iterations: 0,
                success: true,
            });
        }

        // Check small primes
        for &prime in &[3, 5, 7, 11, 13, 17, 19, 23, 29, 31] {
            if self.n % prime == 0 {
                return Some(FactorizationResult {
                    p: prime,
                    q: self.n / prime,
                    confidence: 1.0,
                    iterations: 0,
                    success: true,
                });
            }
        }

        // Check if perfect square
        let sqrt_n = (self.n as f64).sqrt() as u64;
        if sqrt_n * sqrt_n == self.n {
            return Some(FactorizationResult {
                p: sqrt_n,
                q: sqrt_n,
                confidence: 1.0,
                iterations: 0,
                success: true,
            });
        }

        None
    }

    /// Build Hamiltonian for factorization
    fn build_hamiltonian(&self) -> Result<Hamiltonian, String> {
        let mut terms = Vec::new();

        // Binary representation constraints
        // N = p * q where p and q are represented in binary

        // Build multiplication constraints
        for i in 0..self.n_qubits_p {
            for j in 0..self.n_qubits_q {
                let bit_product = (1u64 << i) * (1u64 << j);

                // Check which bits of N this contributes to
                let mut k = 0;
                let mut carry = bit_product;

                while carry > 0 && k < 64 {
                    let n_bit = (self.n >> k) & 1;
                    let expected = carry & 1;

                    if expected == 1 {
                        // Add constraint term
                        let mut pauli_string = vec!['I'; self.n_qubits_p + self.n_qubits_q];
                        pauli_string[i] = 'Z';
                        pauli_string[self.n_qubits_p + j] = 'Z';

                        let coefficient = if n_bit == 1 {
                            -self.penalty_strength
                        } else {
                            self.penalty_strength
                        };

                        terms.push(PauliTerm {
                            coefficient,
                            pauli_string,
                        });
                    }

                    carry >>= 1;
                    k += 1;
                }
            }
        }

        // Add terms to enforce p, q > 1
        if self.use_symmetry_reduction {
            // Fix least significant bit of p to 1 (p is odd)
            let mut pauli_string = vec!['I'; self.n_qubits_p + self.n_qubits_q];
            pauli_string[0] = 'Z';
            terms.push(PauliTerm {
                coefficient: -self.penalty_strength,
                pauli_string,
            });

            // Enforce p <= q by comparing binary representations
            for i in (0..self.n_qubits_p.min(self.n_qubits_q)).rev() {
                let mut pauli_string = vec!['I'; self.n_qubits_p + self.n_qubits_q];
                pauli_string[i] = 'Z';
                pauli_string[self.n_qubits_p + i] = 'Z';

                terms.push(PauliTerm {
                    coefficient: self.penalty_strength * 0.1,
                    pauli_string,
                });
            }
        }

        Ok(Hamiltonian { terms })
    }

    /// Extract factors from VQE solution
    fn extract_factors(&self, params: &[f64]) -> Result<FactorizationResult, String> {
        // Simulate measurement (simplified)
        let mut rng = thread_rng();
        let mut best_p = 0u64;
        let mut best_q = 0u64;
        let mut best_error = self.n as f64;

        // Sample multiple times
        for _ in 0..100 {
            let mut p = 0u64;
            let mut q = 0u64;

            // Extract p
            for (i, &param) in params.iter().enumerate().take(self.n_qubits_p) {
                if rng.random_bool(0.4f64.mul_add(param.sin(), 0.5)) {
                    p |= 1u64 << i;
                }
            }

            // Extract q
            for j in 0..self.n_qubits_q {
                if rng.random_bool(0.4f64.mul_add(params[self.n_qubits_p + j].sin(), 0.5)) {
                    q |= 1u64 << j;
                }
            }

            // Check if valid factorization
            if p > 1 && q > 1 {
                let product = p * q;
                let error = (product as f64 - self.n as f64).abs();

                if error < best_error {
                    best_error = error;
                    best_p = p;
                    best_q = q;

                    if product == self.n {
                        return Ok(FactorizationResult {
                            p: best_p,
                            q: best_q,
                            confidence: 1.0,
                            iterations: 100,
                            success: true,
                        });
                    }
                }
            }
        }

        // Return best approximation
        Ok(FactorizationResult {
            p: best_p,
            q: best_q,
            confidence: 1.0 - (best_error / self.n as f64),
            iterations: 100,
            success: best_p * best_q == self.n,
        })
    }
}

#[derive(Debug, Clone)]
pub struct FactorizationResult {
    pub p: u64,
    pub q: u64,
    pub confidence: f64,
    pub iterations: usize,
    pub success: bool,
}

/// Enhanced VQF with additional techniques
pub struct EnhancedVQF {
    /// Base VQF solver
    base_vqf: VQF,
    /// Use carry handling
    use_carry_handling: bool,
    /// Use modular arithmetic constraints
    use_modular_constraints: bool,
    /// Classical preprocessing depth
    preprocessing_depth: usize,
    /// Multi-level optimization
    use_multilevel: bool,
}

impl EnhancedVQF {
    /// Create new enhanced VQF
    pub fn new(n: u64, optimizer: ClassicalOptimizer) -> Result<Self, String> {
        Ok(Self {
            base_vqf: VQF::new(n, optimizer)?,
            use_carry_handling: true,
            use_modular_constraints: true,
            preprocessing_depth: 3,
            use_multilevel: false,
        })
    }

    /// Enable carry handling
    pub const fn with_carry_handling(mut self, use_carry: bool) -> Self {
        self.use_carry_handling = use_carry;
        self
    }

    /// Enable modular constraints
    pub const fn with_modular_constraints(mut self, use_modular: bool) -> Self {
        self.use_modular_constraints = use_modular;
        self
    }

    /// Enable/disable multi-level optimization
    pub const fn with_multilevel(mut self, use_multilevel: bool) -> Self {
        self.use_multilevel = use_multilevel;
        self
    }

    /// Factor with enhanced techniques
    pub fn factor(&mut self) -> Result<EnhancedFactorizationResult, String> {
        // For multilevel mode, skip classical preprocessing to get full prime factorization.
        if self.use_multilevel {
            return self.multilevel_factorization();
        }

        // Classical preprocessing (single-level path only)
        if let Some(factors) = self.classical_preprocessing() {
            return Ok(EnhancedFactorizationResult {
                factors: vec![factors.p, factors.q],
                method: "Classical preprocessing".to_string(),
                quantum_advantage: 0.0,
                circuit_depth: 0,
            });
        }

        self.single_level_factorization()
    }

    /// Classical preprocessing with multiple techniques
    fn classical_preprocessing(&self) -> Option<FactorizationResult> {
        let n = self.base_vqf.n;

        // Pollard's rho (simplified)
        if self.preprocessing_depth >= 1 {
            if let Some(factor) = self.pollard_rho(n, 1000) {
                return Some(FactorizationResult {
                    p: factor,
                    q: n / factor,
                    confidence: 1.0,
                    iterations: 0,
                    success: true,
                });
            }
        }

        // Trial division with wheel
        if self.preprocessing_depth >= 2 {
            if let Some(factor) = self.wheel_factorization(n) {
                return Some(FactorizationResult {
                    p: factor,
                    q: n / factor,
                    confidence: 1.0,
                    iterations: 0,
                    success: true,
                });
            }
        }

        // Fermat's method for numbers close to perfect square
        if self.preprocessing_depth >= 3 {
            if let Some(factor) = self.fermat_factorization(n) {
                return Some(FactorizationResult {
                    p: factor,
                    q: n / factor,
                    confidence: 1.0,
                    iterations: 0,
                    success: true,
                });
            }
        }

        None
    }

    /// Pollard's rho algorithm
    fn pollard_rho(&self, n: u64, max_iter: usize) -> Option<u64> {
        let mut x = 2u64;
        let mut y = 2u64;
        #[allow(unused_assignments)]
        let mut d = 1u64;

        for _ in 0..max_iter {
            x = (x * x + 1) % n;
            y = (y * y + 1) % n;
            y = (y * y + 1) % n;

            d = gcd(x.abs_diff(y), n);

            if d != 1 && d != n {
                return Some(d);
            }
        }

        None
    }

    /// Wheel factorization
    const fn wheel_factorization(&self, n: u64) -> Option<u64> {
        // Use 2-3-5 wheel
        let wheel = [4, 2, 4, 2, 4, 6, 2, 6];
        let mut k = 7u64;
        let mut i = 0;

        while k * k <= n {
            if n % k == 0 {
                return Some(k);
            }
            k += wheel[i % 8];
            i += 1;
        }

        None
    }

    /// Fermat's factorization method
    fn fermat_factorization(&self, n: u64) -> Option<u64> {
        let mut a = ((n as f64).sqrt().ceil()) as u64;
        let mut b2 = a * a - n;

        for _ in 0..1000 {
            let b = (b2 as f64).sqrt() as u64;
            if b * b == b2 {
                return Some(a - b);
            }
            a += 1;
            b2 = a * a - n;
        }

        None
    }

    /// Single-level quantum factorization
    fn single_level_factorization(&mut self) -> Result<EnhancedFactorizationResult, String> {
        let result = self.base_vqf.factor()?;

        Ok(EnhancedFactorizationResult {
            factors: vec![result.p, result.q],
            method: "Quantum VQF".to_string(),
            quantum_advantage: self.estimate_quantum_advantage(),
            circuit_depth: self.estimate_circuit_depth(),
        })
    }

    /// Multi-level factorization for larger numbers
    fn multilevel_factorization(&mut self) -> Result<EnhancedFactorizationResult, String> {
        let n = self.base_vqf.n;
        let mut factors = Self::factorize_recursive(self, n)?;
        factors.sort_unstable();
        Ok(EnhancedFactorizationResult {
            factors,
            method: "Multilevel recursive".to_string(),
            quantum_advantage: self.estimate_quantum_advantage(),
            circuit_depth: self.estimate_circuit_depth(),
        })
    }

    /// Recursive prime factorization helper
    fn factorize_recursive(&mut self, n: u64) -> Result<Vec<u64>, String> {
        if n <= 1 {
            return Ok(vec![]);
        }
        if is_prime(n) {
            return Ok(vec![n]);
        }
        // Try Pollard's rho
        if let Some(f) = self.pollard_rho(n, 1000) {
            let mut left = self.factorize_recursive(f)?;
            let right = self.factorize_recursive(n / f)?;
            left.extend(right);
            return Ok(left);
        }
        // Try Fermat
        if let Some(f) = self.fermat_factorization(n) {
            if f > 1 && f < n {
                let mut left = self.factorize_recursive(f)?;
                let right = self.factorize_recursive(n / f)?;
                left.extend(right);
                return Ok(left);
            }
        }
        // Try wheel
        if let Some(f) = self.wheel_factorization(n) {
            if f > 1 && f < n {
                let mut left = self.factorize_recursive(f)?;
                let right = self.factorize_recursive(n / f)?;
                left.extend(right);
                return Ok(left);
            }
        }
        // Fall back to base VQF
        self.base_vqf.n = n;
        let result = self.base_vqf.factor()?;
        if result.success && result.p > 1 && result.q > 1 && result.p * result.q == n {
            let mut left = self.factorize_recursive_plain(result.p)?;
            let right = self.factorize_recursive_plain(result.q)?;
            left.extend(right);
            Ok(left)
        } else {
            // Cannot factor further — treat as prime
            Ok(vec![n])
        }
    }

    /// Recursive factorization without VQF fallback (avoids borrowing self.base_vqf)
    fn factorize_recursive_plain(&self, n: u64) -> Result<Vec<u64>, String> {
        if n <= 1 {
            return Ok(vec![]);
        }
        if is_prime(n) {
            return Ok(vec![n]);
        }
        if let Some(f) = self.pollard_rho(n, 1000) {
            let mut left = self.factorize_recursive_plain(f)?;
            let right = self.factorize_recursive_plain(n / f)?;
            left.extend(right);
            return Ok(left);
        }
        if let Some(f) = self.fermat_factorization(n) {
            if f > 1 && f < n {
                let mut left = self.factorize_recursive_plain(f)?;
                let right = self.factorize_recursive_plain(n / f)?;
                left.extend(right);
                return Ok(left);
            }
        }
        if let Some(f) = self.wheel_factorization(n) {
            if f > 1 && f < n {
                let mut left = self.factorize_recursive_plain(f)?;
                let right = self.factorize_recursive_plain(n / f)?;
                left.extend(right);
                return Ok(left);
            }
        }
        Ok(vec![n])
    }

    /// Estimate quantum advantage
    fn estimate_quantum_advantage(&self) -> f64 {
        let n = self.base_vqf.n;
        let classical_complexity = (n as f64).ln() * (n as f64).ln().ln();
        let quantum_complexity = (n as f64).ln().powi(3);

        classical_complexity / quantum_complexity
    }

    /// Estimate circuit depth
    const fn estimate_circuit_depth(&self) -> usize {
        let n_qubits = self.base_vqf.n_qubits_p + self.base_vqf.n_qubits_q;
        n_qubits * 10 // Rough estimate
    }
}

#[derive(Debug, Clone)]
pub struct EnhancedFactorizationResult {
    pub factors: Vec<u64>,
    pub method: String,
    pub quantum_advantage: f64,
    pub circuit_depth: usize,
}

/// Shor's algorithm implementation (simulated)
pub struct ShorsAlgorithm {
    /// Number to factor
    n: u64,
    /// Use quantum period finding
    use_quantum_period_finding: bool,
    /// Classical simulation accuracy
    simulation_shots: usize,
}

impl ShorsAlgorithm {
    /// Create new Shor's algorithm instance
    pub const fn new(n: u64) -> Self {
        Self {
            n,
            use_quantum_period_finding: true,
            simulation_shots: 1000,
        }
    }

    /// Run Shor's algorithm
    pub fn factor(&self) -> Result<ShorsResult, String> {
        // Check if even
        if self.n % 2 == 0 {
            return Ok(ShorsResult {
                factors: vec![2, self.n / 2],
                period: 2,
                success_probability: 1.0,
            });
        }

        // Check if perfect power
        if let Some((base, exp)) = self.is_perfect_power() {
            return Ok(ShorsResult {
                factors: vec![base; exp],
                period: 0,
                success_probability: 1.0,
            });
        }

        let mut rng = thread_rng();

        // Try random bases
        for attempt in 0..10 {
            // Choose random a coprime to n
            let a = rng.random_range(2..self.n);
            if gcd(a, self.n) != 1 {
                let factor = gcd(a, self.n);
                return Ok(ShorsResult {
                    factors: vec![factor, self.n / factor],
                    period: 0,
                    success_probability: 1.0,
                });
            }

            // Find period of a^x mod n
            let period = if self.use_quantum_period_finding {
                self.quantum_period_finding(a)?
            } else {
                self.classical_period_finding(a)
            };

            if period > 0 && period % 2 == 0 {
                let factor1 = gcd(mod_exp(a, period / 2, self.n) - 1, self.n);
                let factor2 = gcd(mod_exp(a, period / 2, self.n) + 1, self.n);

                if factor1 > 1 && factor1 < self.n {
                    return Ok(ShorsResult {
                        factors: vec![factor1, self.n / factor1],
                        period,
                        success_probability: 0.05f64.mul_add(-(attempt as f64), 0.9),
                    });
                }

                if factor2 > 1 && factor2 < self.n {
                    return Ok(ShorsResult {
                        factors: vec![factor2, self.n / factor2],
                        period,
                        success_probability: 0.05f64.mul_add(-(attempt as f64), 0.9),
                    });
                }
            }
        }

        Err("Failed to find factors".to_string())
    }

    /// Check if n is a perfect power
    fn is_perfect_power(&self) -> Option<(u64, usize)> {
        for exp in 2..64 {
            let base = (self.n as f64).powf(1.0 / exp as f64) as u64;

            if base.pow(exp) == self.n {
                return Some((base, exp as usize));
            }

            if base <= 1 {
                break;
            }
        }

        None
    }

    /// Quantum period finding (simulated)
    fn quantum_period_finding(&self, a: u64) -> Result<u64, String> {
        // Simulate quantum Fourier transform
        // This is a classical simulation of the quantum algorithm

        let max_period = self.n.min(1000);
        let mut amplitudes = vec![0.0; max_period as usize];

        // Prepare superposition and apply modular exponentiation
        let amplitudes_len = amplitudes.len();
        for x in 0..max_period {
            let state = mod_exp(a, x, self.n);
            amplitudes[state as usize % amplitudes_len] += 1.0;
        }

        // Find period from amplitude pattern
        for period in 1..max_period {
            let mut is_period = true;

            for i in 0..amplitudes.len() {
                if amplitudes[i] > 0.0
                    && amplitudes[(i + period as usize) % amplitudes.len()] == 0.0
                {
                    is_period = false;
                    break;
                }
            }

            if is_period && mod_exp(a, period, self.n) == 1 {
                return Ok(period);
            }
        }

        Ok(0)
    }

    /// Classical period finding
    fn classical_period_finding(&self, a: u64) -> u64 {
        let mut x = 1u64;

        for period in 1..self.n.min(10000) {
            x = (x * a) % self.n;

            if x == 1 {
                return period;
            }
        }

        0
    }
}

#[derive(Debug, Clone)]
pub struct ShorsResult {
    pub factors: Vec<u64>,
    pub period: u64,
    pub success_probability: f64,
}

/// Primality test via trial division up to √n
fn is_prime(n: u64) -> bool {
    if n < 2 {
        return false;
    }
    if n == 2 || n == 3 {
        return true;
    }
    if n % 2 == 0 || n % 3 == 0 {
        return false;
    }
    let mut k = 5u64;
    while k * k <= n {
        if n % k == 0 || n % (k + 2) == 0 {
            return false;
        }
        k += 6;
    }
    true
}

/// Helper functions
const fn gcd(mut a: u64, mut b: u64) -> u64 {
    while b != 0 {
        let temp = b;
        b = a % b;
        a = temp;
    }
    a
}

const fn mod_exp(base: u64, exp: u64, modulus: u64) -> u64 {
    let mut result = 1u64;
    let mut base = base % modulus;
    let mut exp = exp;

    while exp > 0 {
        if exp % 2 == 1 {
            result = (result * base) % modulus;
        }
        exp >>= 1;
        base = (base * base) % modulus;
    }

    result
}

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

    #[test]
    fn test_vqf_small() {
        let optimizer = ClassicalOptimizer::GradientDescent {
            learning_rate: 0.1,
            momentum: 0.0,
        };

        let vqf = VQF::new(15, optimizer).expect("Failed to create VQF for n=15");
        assert_eq!(vqf.n, 15);

        // Test preprocessing
        let mut result = vqf.preprocess();
        assert!(result.is_some());
        if let Some(factors) = result {
            assert_eq!(factors.p * factors.q, 15);
        }
    }

    #[test]
    fn test_enhanced_vqf() {
        let optimizer = ClassicalOptimizer::SPSA {
            a: 0.1,
            c: 0.1,
            alpha: 0.602,
            gamma: 0.101,
        };

        let enhanced =
            EnhancedVQF::new(21, optimizer).expect("Failed to create EnhancedVQF for n=21");
        assert!(enhanced.use_carry_handling);
    }

    #[test]
    fn test_shors_algorithm() {
        let shors = ShorsAlgorithm::new(15);
        let mut result = shors.factor();

        assert!(result.is_ok());
        let factors = result.expect("Shor's algorithm should factor 15").factors;
        assert_eq!(factors[0] * factors[1], 15);
    }

    #[test]
    fn test_helper_functions() {
        assert_eq!(gcd(48, 18), 6);
        assert_eq!(mod_exp(3, 4, 7), 4); // 3^4 mod 7 = 81 mod 7 = 4
    }

    #[test]
    fn test_multilevel_factor_15() {
        let optimizer = ClassicalOptimizer::GradientDescent {
            learning_rate: 0.1,
            momentum: 0.0,
        };
        let mut vqf = EnhancedVQF::new(15, optimizer)
            .expect("Failed to create EnhancedVQF for n=15")
            .with_multilevel(true);
        let result = vqf
            .factor()
            .expect("multilevel_factorization(15) should succeed");
        assert!(result.factors.contains(&3), "factors should contain 3");
        assert!(result.factors.contains(&5), "factors should contain 5");
    }

    #[test]
    fn test_multilevel_factor_prime() {
        let optimizer = ClassicalOptimizer::GradientDescent {
            learning_rate: 0.1,
            momentum: 0.0,
        };
        let mut vqf = EnhancedVQF::new(13, optimizer)
            .expect("Failed to create EnhancedVQF for n=13")
            .with_multilevel(true);
        let result = vqf
            .factor()
            .expect("multilevel_factorization(13) should succeed");
        assert_eq!(result.factors, vec![13], "prime 13 should return [13]");
    }

    #[test]
    fn test_multilevel_factor_composite() {
        let optimizer = ClassicalOptimizer::GradientDescent {
            learning_rate: 0.1,
            momentum: 0.0,
        };
        let mut vqf = EnhancedVQF::new(105, optimizer)
            .expect("Failed to create EnhancedVQF for n=105")
            .with_multilevel(true);
        let result = vqf
            .factor()
            .expect("multilevel_factorization(105) should succeed");
        assert!(result.factors.contains(&3), "factors should contain 3");
        assert!(result.factors.contains(&5), "factors should contain 5");
        assert!(result.factors.contains(&7), "factors should contain 7");
    }
}