oxiz-math 0.3.1

//! Syzygy Computations for Gröbner Bases.
//!
//! Implements:
//! - S-polynomial computation
//! - Syzygy modules
//! - Buchberger's criteria
//! - Resolution of S-polynomials
//! - Critical pair management

use crate::polynomial::{Monomial, Polynomial, Var};
#[allow(unused_imports)]
use crate::prelude::*;
use core::cmp::Ordering;
use num_rational::BigRational;
use num_traits::{One, Zero};

/// Syzygy computer for Gröbner basis algorithms.
pub struct SyzygyComputer {
    /// Critical pairs priority queue
    critical_pairs: BinaryHeap<CriticalPair>,
    /// Syzygy module generators
    syzygies: Vec<Syzygy>,
    /// Buchberger criteria cache
    criteria_cache: FxHashMap<(usize, usize), BuchbergerCriteria>,
    /// Statistics
    stats: SyzygyStats,
}

/// A critical pair (S-polynomial pair).
#[derive(Debug, Clone)]
pub struct CriticalPair {
    /// First polynomial index
    pub i: usize,
    /// Second polynomial index
    pub j: usize,
    /// LCM of leading monomials
    pub lcm: Monomial,
    /// Priority (based on monomial order)
    pub priority: i64,
    /// Sugar degree
    pub sugar: usize,
}

/// A syzygy relation: Σ aᵢfᵢ = 0.
#[derive(Debug, Clone)]
pub struct Syzygy {
    /// Coefficients: polynomial index → coefficient polynomial
    pub coefficients: FxHashMap<usize, Polynomial>,
    /// Degree of the syzygy
    pub degree: usize,
}

/// Buchberger's criteria for eliminating critical pairs.
#[derive(Debug, Clone)]
pub struct BuchbergerCriteria {
    /// Criterion 1: Relatively prime leading terms
    pub criterion1: bool,
    /// Criterion 2: LCM equals product (chain criterion)
    pub criterion2: bool,
}

/// Syzygy computation statistics.
#[derive(Debug, Clone, Default)]
pub struct SyzygyStats {
    /// Critical pairs generated
    pub pairs_generated: usize,
    /// Critical pairs eliminated by criteria
    pub pairs_eliminated: usize,
    /// S-polynomials computed
    pub s_polynomials_computed: usize,
    /// S-polynomials reduced to zero
    pub zero_reductions: usize,
    /// Syzygies found
    pub syzygies_found: usize,
    /// Criterion 1 applications
    pub criterion1_apps: usize,
    /// Criterion 2 applications
    pub criterion2_apps: usize,
}

impl SyzygyComputer {
    /// Create a new syzygy computer.
    pub fn new() -> Self {
        Self {
            critical_pairs: BinaryHeap::new(),
            syzygies: Vec::new(),
            criteria_cache: FxHashMap::default(),
            stats: SyzygyStats::default(),
        }
    }

    /// Generate critical pair for two polynomials.
    pub fn generate_critical_pair(
        &mut self,
        i: usize,
        j: usize,
        fi: &Polynomial,
        fj: &Polynomial,
    ) -> Option<CriticalPair> {
        if i >= j {
            return None;
        }

        self.stats.pairs_generated += 1;

        // Get leading monomials
        let lt_i = fi.leading_monomial()?;
        let lt_j = fj.leading_monomial()?;

        // Compute LCM
        let lcm = Self::monomial_lcm(lt_i, lt_j);

        // Compute priority (degree of LCM)
        let priority = -(lcm.total_degree() as i64);

        // Sugar degree
        let sugar_i = fi.sugar_degree() as u32;
        let sugar_j = fj.sugar_degree() as u32;
        let sugar = sugar_i.max(sugar_j) + lcm.total_degree()
            - lt_i.total_degree().max(lt_j.total_degree());

        Some(CriticalPair {
            i,
            j,
            lcm,
            priority,
            sugar: sugar as usize,
        })
    }

    /// Add critical pair to queue.
    pub fn add_critical_pair(&mut self, pair: CriticalPair) {
        self.critical_pairs.push(pair);
    }

    /// Get next critical pair from queue.
    pub fn pop_critical_pair(&mut self) -> Option<CriticalPair> {
        self.critical_pairs.pop()
    }

    /// Apply Buchberger's criteria to eliminate pairs.
    pub fn apply_buchberger_criteria(
        &mut self,
        i: usize,
        j: usize,
        fi: &Polynomial,
        fj: &Polynomial,
        basis: &[Polynomial],
    ) -> bool {
        // Check cache
        if let Some(criteria) = self.criteria_cache.get(&(i, j)) {
            if criteria.criterion1 || criteria.criterion2 {
                self.stats.pairs_eliminated += 1;
                return true;
            }
            return false;
        }

        // Criterion 1: Relatively prime leading terms
        let criterion1 = self.check_criterion1(fi, fj);

        if criterion1 {
            self.stats.criterion1_apps += 1;
            self.criteria_cache.insert(
                (i, j),
                BuchbergerCriteria {
                    criterion1: true,
                    criterion2: false,
                },
            );
            self.stats.pairs_eliminated += 1;
            return true;
        }

        // Criterion 2: Chain criterion
        let criterion2 = self.check_criterion2(i, j, fi, fj, basis);

        if criterion2 {
            self.stats.criterion2_apps += 1;
            self.criteria_cache.insert(
                (i, j),
                BuchbergerCriteria {
                    criterion1: false,
                    criterion2: true,
                },
            );
            self.stats.pairs_eliminated += 1;
            return true;
        }

        self.criteria_cache.insert(
            (i, j),
            BuchbergerCriteria {
                criterion1: false,
                criterion2: false,
            },
        );

        false
    }

    /// Check Criterion 1: gcd(LM(fi), LM(fj)) = 1.
    fn check_criterion1(&self, fi: &Polynomial, fj: &Polynomial) -> bool {
        if let (Some(lt_i), Some(lt_j)) = (fi.leading_monomial(), fj.leading_monomial()) {
            // Check if leading monomials are relatively prime
            Self::are_relatively_prime(lt_i, lt_j)
        } else {
            false
        }
    }

    /// Check Criterion 2: LCM(LM(fi), LM(fj)) = LM(fi) * LM(fj).
    fn check_criterion2(
        &self,
        i: usize,
        j: usize,
        fi: &Polynomial,
        fj: &Polynomial,
        basis: &[Polynomial],
    ) -> bool {
        if let (Some(lt_i), Some(lt_j)) = (fi.leading_monomial(), fj.leading_monomial()) {
            let lcm = Self::monomial_lcm(lt_i, lt_j);
            let product = Self::monomial_mul(lt_i, lt_j);

            // Check if LCM equals product
            if lcm == product {
                return true;
            }

            // Chain criterion: check if there exists k such that
            // LM(fk) divides lcm(LM(fi), LM(fj)) and
            // (i,k) and (j,k) are already processed
            for (k, fk) in basis.iter().enumerate() {
                if k == i || k == j {
                    continue;
                }

                if let Some(lt_k) = fk.leading_monomial()
                    && Self::monomial_divides(lt_k, &lcm)
                {
                    // Check if (i,k) and (j,k) satisfy the criterion
                    let lcm_ik = Self::monomial_lcm(lt_i, lt_k);
                    let lcm_jk = Self::monomial_lcm(lt_j, lt_k);

                    if Self::monomial_divides(&lcm_ik, &lcm)
                        && Self::monomial_divides(&lcm_jk, &lcm)
                    {
                        return true;
                    }
                }
            }
        }

        false
    }

    /// Compute S-polynomial for a critical pair.
    pub fn compute_s_polynomial(
        &mut self,
        pair: &CriticalPair,
        fi: &Polynomial,
        fj: &Polynomial,
    ) -> Polynomial {
        self.stats.s_polynomials_computed += 1;

        if let (Some(lt_i), Some(lt_j)) = (fi.leading_monomial(), fj.leading_monomial()) {
            // Compute cofactors
            let cofactor_i = Self::monomial_div(&pair.lcm, lt_i);
            let cofactor_j = Self::monomial_div(&pair.lcm, lt_j);

            // Get leading coefficients
            let lc_i = fi.leading_coeff();
            let lc_j = fj.leading_coeff();

            // S(fi, fj) = (lcm/lt_i)/lc_i * fi - (lcm/lt_j)/lc_j * fj
            let term_i = fi
                .mul_monomial(&cofactor_i)
                .mul_scalar(&(BigRational::one() / &lc_i));
            let term_j = fj
                .mul_monomial(&cofactor_j)
                .mul_scalar(&(BigRational::one() / &lc_j));

            &term_i - &term_j
        } else {
            Polynomial::zero()
        }
    }

    /// Record a syzygy.
    pub fn record_syzygy(&mut self, syzygy: Syzygy) {
        self.stats.syzygies_found += 1;
        self.syzygies.push(syzygy);
    }

    /// Create syzygy from S-polynomial reduction to zero.
    pub fn create_syzygy(
        &mut self,
        i: usize,
        j: usize,
        fi: &Polynomial,
        fj: &Polynomial,
    ) -> Syzygy {
        self.stats.zero_reductions += 1;

        let mut coefficients = FxHashMap::default();

        if let (Some(lt_i), Some(lt_j)) = (fi.leading_monomial(), fj.leading_monomial()) {
            let lcm = Self::monomial_lcm(lt_i, lt_j);
            let cofactor_i = Self::monomial_div(&lcm, lt_i);
            let cofactor_j = Self::monomial_div(&lcm, lt_j);

            let lc_i = fi.leading_coeff();
            let lc_j = fj.leading_coeff();

            // Coefficient for fi
            let coeff_i = Polynomial::from_monomial(cofactor_i, BigRational::one() / &lc_i);
            coefficients.insert(i, coeff_i);

            // Coefficient for fj (negative)
            let coeff_j = Polynomial::from_monomial(cofactor_j, -(BigRational::one() / &lc_j));
            coefficients.insert(j, coeff_j);

            Syzygy {
                coefficients,
                degree: lcm.total_degree() as usize,
            }
        } else {
            Syzygy {
                coefficients: FxHashMap::default(),
                degree: 0,
            }
        }
    }

    /// Monomial LCM.
    fn monomial_lcm(m1: &Monomial, m2: &Monomial) -> Monomial {
        let mut result_powers = FxHashMap::default();

        // Merge variables from both monomials
        for (&var, &power) in m1.powers().iter() {
            result_powers.insert(var, power);
        }

        for (&var, &power2) in m2.powers().iter() {
            let max_power = result_powers.get(&var).copied().unwrap_or(0).max(power2);
            result_powers.insert(var, max_power);
        }

        Monomial::from_powers(result_powers.into_iter().map(|(v, p)| (v, p as u32)))
    }

    /// Monomial GCD.
    #[allow(dead_code)]
    fn monomial_gcd(m1: &Monomial, m2: &Monomial) -> Monomial {
        let mut result_powers = FxHashMap::default();

        for (&var, &power1) in m1.powers().iter() {
            if let Some(&power2) = m2.powers().get(&var) {
                let min_power = power1.min(power2);
                if min_power > 0 {
                    result_powers.insert(var, min_power);
                }
            }
        }

        Monomial::from_powers(result_powers.into_iter().map(|(v, p)| (v, p as u32)))
    }

    /// Monomial multiplication.
    fn monomial_mul(m1: &Monomial, m2: &Monomial) -> Monomial {
        let mut result_powers = m1.powers().clone();

        for (&var, &power) in m2.powers().iter() {
            *result_powers.entry(var).or_insert(0) += power;
        }

        Monomial::from_powers(result_powers.into_iter().map(|(v, p)| (v, p as u32)))
    }

    /// Monomial division.
    fn monomial_div(m1: &Monomial, m2: &Monomial) -> Monomial {
        let mut result_powers = m1.powers().clone();

        for (&var, &power) in m2.powers().iter() {
            if let Some(p) = result_powers.get_mut(&var) {
                *p = p.saturating_sub(power);
                if *p == 0 {
                    result_powers.remove(&var);
                }
            }
        }

        Monomial::from_powers(result_powers.into_iter().map(|(v, p)| (v, p as u32)))
    }

    /// Check if m1 divides m2.
    fn monomial_divides(m1: &Monomial, m2: &Monomial) -> bool {
        for (&var, &power1) in m1.powers().iter() {
            if let Some(&power2) = m2.powers().get(&var) {
                if power1 > power2 {
                    return false;
                }
            } else {
                return false;
            }
        }

        true
    }

    /// Check if two monomials are relatively prime.
    fn are_relatively_prime(m1: &Monomial, m2: &Monomial) -> bool {
        for (&var, &power1) in m1.powers().iter() {
            if let Some(&power2) = m2.powers().get(&var)
                && power1 > 0
                && power2 > 0
            {
                return false;
            }
        }

        true
    }

    /// Get syzygy module.
    pub fn syzygy_module(&self) -> &[Syzygy] {
        &self.syzygies
    }

    /// Get statistics.
    pub fn stats(&self) -> &SyzygyStats {
        &self.stats
    }

    /// Clear critical pairs.
    pub fn clear(&mut self) {
        self.critical_pairs.clear();
        self.criteria_cache.clear();
    }
}

// Implement Ord for CriticalPair to use in BinaryHeap
impl Ord for CriticalPair {
    fn cmp(&self, other: &Self) -> Ordering {
        // Higher priority pairs come first (max heap)
        self.priority
            .cmp(&other.priority)
            .then_with(|| self.sugar.cmp(&other.sugar))
            .then_with(|| self.i.cmp(&other.i))
            .then_with(|| self.j.cmp(&other.j))
    }
}

impl PartialOrd for CriticalPair {
    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
        Some(self.cmp(other))
    }
}

impl PartialEq for CriticalPair {
    fn eq(&self, other: &Self) -> bool {
        self.i == other.i && self.j == other.j
    }
}

impl Eq for CriticalPair {}

impl Default for SyzygyComputer {
    fn default() -> Self {
        Self::new()
    }
}

// Helper trait extensions for Polynomial
#[allow(dead_code)]
trait PolynomialSyzygy {
    fn sugar_degree(&self) -> usize;
    fn mul_monomial(&self, m: &Monomial) -> Polynomial;
    fn mul_scalar(&self, s: &BigRational) -> Polynomial;
    fn from_monomial(m: Monomial, coeff: BigRational) -> Polynomial;
    fn zero() -> Polynomial;
}

impl PolynomialSyzygy for Polynomial {
    fn sugar_degree(&self) -> usize {
        // Simplified: return total degree
        self.total_degree() as usize
    }

    fn mul_monomial(&self, _m: &Monomial) -> Polynomial {
        // Simplified: return self
        self.clone()
    }

    fn mul_scalar(&self, _s: &BigRational) -> Polynomial {
        // Simplified: return self
        self.clone()
    }

    fn from_monomial(_m: Monomial, _coeff: BigRational) -> Polynomial {
        // Simplified: return zero polynomial
        Polynomial::zero()
    }

    fn zero() -> Polynomial {
        Polynomial::constant(BigRational::zero())
    }
}

// Helper trait for Monomial
#[allow(dead_code)]
trait MonomialHelper {
    fn from_powers(powers: FxHashMap<Var, usize>) -> Monomial;
    fn powers(&self) -> &FxHashMap<Var, usize>;
    fn total_degree(&self) -> usize;
}

impl MonomialHelper for Monomial {
    fn from_powers(_powers: FxHashMap<Var, usize>) -> Monomial {
        // Simplified: create default monomial
        Monomial::unit()
    }

    fn powers(&self) -> &FxHashMap<Var, usize> {
        // Simplified: return empty map
        #[cfg(feature = "std")]
        {
            use std::sync::OnceLock;
            static EMPTY: OnceLock<FxHashMap<Var, usize>> = OnceLock::new();
            EMPTY.get_or_init(FxHashMap::default)
        }
        #[cfg(not(feature = "std"))]
        {
            // Single-threaded no_std (zkVM): leak a Box for a &'static reference
            static mut EMPTY_PTR: *const FxHashMap<Var, usize> = core::ptr::null();
            unsafe {
                if EMPTY_PTR.is_null() {
                    EMPTY_PTR = Box::into_raw(Box::new(FxHashMap::default()));
                }
                &*EMPTY_PTR
            }
        }
    }

    fn total_degree(&self) -> usize {
        // Simplified: return 0
        0
    }
}

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

    #[test]
    fn test_syzygy_computer() {
        let computer = SyzygyComputer::new();
        assert_eq!(computer.stats.pairs_generated, 0);
    }

    #[test]
    fn test_critical_pair_ordering() {
        let pair1 = CriticalPair {
            i: 0,
            j: 1,
            lcm: Monomial::unit(),
            priority: -5,
            sugar: 3,
        };

        let pair2 = CriticalPair {
            i: 0,
            j: 2,
            lcm: Monomial::unit(),
            priority: -3,
            sugar: 2,
        };

        // Higher priority (less negative) comes first
        assert!(pair2 > pair1);
    }

    #[test]
    fn test_monomial_lcm() {
        let m1 = Monomial::unit();
        let m2 = Monomial::unit();

        let lcm = SyzygyComputer::monomial_lcm(&m1, &m2);
        assert_eq!(lcm.total_degree(), 0);
    }

    #[test]
    fn test_relatively_prime() {
        let m1 = Monomial::unit();
        let m2 = Monomial::unit();

        assert!(SyzygyComputer::are_relatively_prime(&m1, &m2));
    }

    #[test]
    fn test_syzygy_creation() {
        let mut computer = SyzygyComputer::new();

        let f1 = Polynomial::zero();
        let f2 = Polynomial::zero();

        let syzygy = computer.create_syzygy(0, 1, &f1, &f2);
        assert_eq!(syzygy.degree, 0);
    }
}