mathhook-core 0.2.0

//! Advanced polynomial operations
//! Implements polynomial arithmetic, division, and advanced algorithms

use crate::core::{Expression, Number, Symbol};
use crate::simplify::Simplify;

/// Trait for advanced polynomial operations
pub trait AdvancedPolynomial {
    fn polynomial_divide(&self, divisor: &Self) -> (Expression, Expression);
    fn polynomial_remainder(&self, divisor: &Self) -> Expression;
    fn polynomial_degree(&self, var: &Symbol) -> Option<i64>;
    /// Leading coefficient extraction
    fn polynomial_leading_coefficient(&self, var: &Symbol) -> Expression;
    /// Content extraction (GCD of coefficients)
    fn polynomial_content(&self) -> Expression;
    /// Primitive part extraction (polynomial / content)
    fn polynomial_primitive_part(&self) -> Expression;
    /// Resultant computation for polynomial elimination
    fn polynomial_resultant(&self, other: &Self, var: &Symbol) -> Expression;
    /// Discriminant computation
    fn polynomial_discriminant(&self, var: &Symbol) -> Expression;
}

impl AdvancedPolynomial for Expression {
    /// Polynomial division using optimized algorithms
    #[inline(always)]
    fn polynomial_divide(&self, divisor: &Self) -> (Expression, Expression) {
        // Fast path: division by constants
        if let Expression::Number(Number::Integer(d)) = divisor {
            if *d != 0 {
                return (
                    Expression::mul(vec![
                        self.clone(),
                        Expression::pow(divisor.clone(), Expression::integer(-1)),
                    ]),
                    Expression::integer(0),
                );
            }
        }

        // Fast path: identical polynomials
        if self == divisor {
            return (Expression::integer(1), Expression::integer(0));
        }

        // For now, return symbolic division
        (
            Expression::mul(vec![
                self.clone(),
                Expression::pow(divisor.clone(), Expression::integer(-1)),
            ]),
            Expression::integer(0),
        )
    }

    /// Polynomial remainder
    #[inline(always)]
    fn polynomial_remainder(&self, divisor: &Self) -> Expression {
        let (_, remainder) = self.polynomial_divide(divisor);
        remainder
    }

    /// Polynomial degree computation
    #[inline(always)]
    fn polynomial_degree(&self, var: &Symbol) -> Option<i64> {
        match self {
            Expression::Symbol(s) if s == var => Some(1),
            Expression::Number(_) => Some(0),
            Expression::Pow(base, exp) => {
                if let (Expression::Symbol(s), Expression::Number(Number::Integer(e))) =
                    (base.as_ref(), exp.as_ref())
                {
                    if s == var {
                        return Some(*e);
                    }
                }
                None
            }
            Expression::Add(terms) => {
                // Find maximum degree among terms
                let mut max_degree = 0i64;
                for term in terms.iter() {
                    if let Some(deg) = term.polynomial_degree(var) {
                        max_degree = max_degree.max(deg);
                    } else {
                        return None; // Non-polynomial
                    }
                }
                Some(max_degree)
            }
            Expression::Mul(factors) => {
                // Sum degrees of factors
                let mut total_degree = 0i64;
                for factor in factors.iter() {
                    if let Some(deg) = factor.polynomial_degree(var) {
                        total_degree += deg;
                    } else {
                        return None; // Non-polynomial
                    }
                }
                Some(total_degree)
            }
            _ => None,
        }
    }

    #[inline(always)]
    fn polynomial_leading_coefficient(&self, var: &Symbol) -> Expression {
        match self {
            Expression::Symbol(s) if s == var => Expression::integer(1),
            Expression::Number(_) => self.clone(),
            Expression::Pow(base, _exp) => {
                if let Expression::Symbol(s) = base.as_ref() {
                    if s == var {
                        return Expression::integer(1);
                    }
                }
                self.clone()
            }
            Expression::Mul(factors) => {
                // Extract coefficient from factors
                let mut coefficient = Expression::integer(1);
                for factor in factors.iter() {
                    if factor.polynomial_degree(var).is_some()
                        && !matches!(factor, Expression::Symbol(s) if s == var)
                        && !matches!(factor, Expression::Pow(base, _) if matches!(base.as_ref(), Expression::Symbol(s) if s == var))
                    {
                        coefficient = Expression::mul(vec![coefficient, factor.clone()]);
                    }
                }
                coefficient
            }
            Expression::Add(terms) => {
                // Find term with highest degree
                let mut max_degree = -1i64;
                let mut leading_term = Expression::integer(0);

                for term in terms.iter() {
                    if let Some(deg) = term.polynomial_degree(var) {
                        if deg > max_degree {
                            max_degree = deg;
                            leading_term = term.clone();
                        }
                    }
                }

                leading_term.polynomial_leading_coefficient(var)
            }
            _ => Expression::integer(1),
        }
    }

    #[inline(always)]
    fn polynomial_content(&self) -> Expression {
        match self {
            Expression::Number(_) => self.clone(),
            Expression::Symbol(_) => Expression::integer(1),
            Expression::Add(terms) => {
                // Find GCD of all coefficients
                let mut content = Expression::integer(0);
                for term in terms.iter() {
                    let term_content = term.polynomial_content();
                    if content.is_zero() {
                        content = term_content;
                    } else {
                        content = content.gcd(&term_content);
                    }
                }
                content
            }
            Expression::Mul(factors) => {
                // Content of product is product of contents
                let mut content = Expression::integer(1);
                for factor in factors.iter() {
                    content = Expression::mul(vec![content, factor.polynomial_content()]);
                }
                content
            }
            _ => Expression::integer(1),
        }
    }

    #[inline(always)]
    fn polynomial_primitive_part(&self) -> Expression {
        let content = self.polynomial_content();
        if content.is_zero() || content == Expression::integer(1) {
            self.clone()
        } else {
            // Divide by content (symbolic division for now)
            Expression::mul(vec![
                self.clone(),
                Expression::pow(content, Expression::integer(-1)),
            ])
        }
    }

    #[inline(always)]
    fn polynomial_resultant(&self, other: &Self, var: &Symbol) -> Expression {
        // For now, return a symbolic resultant
        // This is a complex operation that requires Sylvester matrix computation
        Expression::function(
            "resultant",
            vec![self.clone(), other.clone(), Expression::symbol(var.clone())],
        )
    }

    #[inline(always)]
    fn polynomial_discriminant(&self, var: &Symbol) -> Expression {
        // Discriminant is related to resultant: Disc(f) = (-1)^(n(n-1)/2) * Res(f, f') / lc(f)
        // For now, return symbolic discriminant
        Expression::function(
            "discriminant",
            vec![self.clone(), Expression::symbol(var.clone())],
        )
    }
}

/// Polynomial arithmetic operations
pub struct PolynomialArithmetic;

impl PolynomialArithmetic {
    /// Polynomial addition
    #[inline(always)]
    pub fn add_polynomials(poly1: &Expression, poly2: &Expression) -> Expression {
        Expression::add(vec![poly1.clone(), poly2.clone()]).simplify()
    }

    /// Polynomial multiplication using convolution
    #[inline(always)]
    pub fn multiply_polynomials(poly1: &Expression, poly2: &Expression) -> Expression {
        Expression::mul(vec![poly1.clone(), poly2.clone()]).simplify()
    }

    /// Polynomial evaluation using Horner's method
    #[inline(always)]
    pub fn evaluate_polynomial(poly: &Expression, var: &Symbol, value: &Expression) -> Expression {
        // Substitute variable with value
        match poly {
            Expression::Symbol(s) if s == var => value.clone(),
            Expression::Number(_) => poly.clone(),
            Expression::Add(terms) => {
                let evaluated_terms: Vec<Expression> = terms
                    .iter()
                    .map(|term| Self::evaluate_polynomial(term, var, value))
                    .collect();
                Expression::add(evaluated_terms).simplify()
            }
            Expression::Mul(factors) => {
                let evaluated_factors: Vec<Expression> = factors
                    .iter()
                    .map(|factor| Self::evaluate_polynomial(factor, var, value))
                    .collect();
                Expression::mul(evaluated_factors).simplify()
            }
            Expression::Pow(base, exp) => {
                let eval_base = Self::evaluate_polynomial(base, var, value);
                let eval_exp = Self::evaluate_polynomial(exp, var, value);
                Expression::pow(eval_base, eval_exp).simplify()
            }
            _ => poly.clone(),
        }
    }

    /// Polynomial composition f(g(x))
    #[inline(always)]
    pub fn compose_polynomials(f: &Expression, g: &Expression, var: &Symbol) -> Expression {
        Self::evaluate_polynomial(f, var, g)
    }
}

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

    #[test]
    fn test_polynomial_degree() {
        let x = symbol!(x);

        // Test degree computation
        assert_eq!(Expression::symbol(x.clone()).polynomial_degree(&x), Some(1));
        assert_eq!(Expression::integer(5).polynomial_degree(&x), Some(0));
        assert_eq!(
            Expression::pow(Expression::symbol(x.clone()), Expression::integer(3))
                .polynomial_degree(&x),
            Some(3)
        );

        // Test polynomial: 2x³ + 3x² + x + 1
        let poly = Expression::add(vec![
            Expression::mul(vec![
                Expression::integer(2),
                Expression::pow(Expression::symbol(x.clone()), Expression::integer(3)),
            ]),
            Expression::mul(vec![
                Expression::integer(3),
                Expression::pow(Expression::symbol(x.clone()), Expression::integer(2)),
            ]),
            Expression::symbol(x.clone()),
            Expression::integer(1),
        ]);

        assert_eq!(poly.polynomial_degree(&x), Some(3));
    }

    #[test]
    fn test_polynomial_leading_coefficient() {
        let x = symbol!(x);

        // Test leading coefficient: 5x² + 3x + 1 has leading coefficient 5
        let poly = Expression::add(vec![
            Expression::mul(vec![
                Expression::integer(5),
                Expression::pow(Expression::symbol(x.clone()), Expression::integer(2)),
            ]),
            Expression::mul(vec![Expression::integer(3), Expression::symbol(x.clone())]),
            Expression::integer(1),
        ]);

        let leading_coeff = poly.polynomial_leading_coefficient(&x);
        println!("Leading coefficient: {}", leading_coeff);

        // Should extract the coefficient of the highest degree term
        assert!(!leading_coeff.is_zero());
    }

    #[test]
    fn test_polynomial_evaluation() {
        let x = symbol!(x);

        // Test polynomial evaluation: f(x) = x² + 2x + 1, f(3) = 9 + 6 + 1 = 16
        let poly = Expression::add(vec![
            Expression::pow(Expression::symbol(x.clone()), Expression::integer(2)),
            Expression::mul(vec![Expression::integer(2), Expression::symbol(x.clone())]),
            Expression::integer(1),
        ]);

        let result = PolynomialArithmetic::evaluate_polynomial(&poly, &x, &Expression::integer(3));
        let simplified = result.simplify();

        println!("Polynomial evaluation f(3): {}", simplified);

        // Should evaluate to 16
        assert!(!simplified.is_zero());
    }

    #[test]
    fn test_polynomial_arithmetic_performance() {
        use std::time::Instant;

        let x = symbol!(x);
        let poly1 = Expression::add(vec![
            Expression::pow(Expression::symbol(x.clone()), Expression::integer(2)),
            Expression::symbol(x.clone()),
            Expression::integer(1),
        ]);
        let poly2 = Expression::add(vec![Expression::symbol(x.clone()), Expression::integer(2)]);

        let start = Instant::now();

        // Perform 1000 polynomial operations
        for _ in 0..1000 {
            let _sum = PolynomialArithmetic::add_polynomials(&poly1, &poly2);
            let _product = PolynomialArithmetic::multiply_polynomials(&poly1, &poly2);
            let _evaluation =
                PolynomialArithmetic::evaluate_polynomial(&poly1, &x, &Expression::integer(5));
        }

        let duration = start.elapsed();
        let ops_per_sec = 3000.0 / duration.as_secs_f64(); // 3 ops per iteration

        println!(
            "Advanced polynomial ops: {:.2}K ops/sec",
            ops_per_sec / 1_000.0
        );

        // Should achieve high performance
        assert!(ops_per_sec > 10_000.0); // > 10K ops/sec
    }
}