thales 0.4.2

//! Algebraic equation solver with symbolic manipulation.
//!
//! This module provides a comprehensive framework for solving algebraic equations
//! symbolically. It supports linear, quadratic, polynomial, and transcendental
//! equations, with automatic method selection via the [`SmartSolver`].
//!
//! # Overview
//!
//! The solver works by:
//! 1. Analyzing the equation structure to determine appropriate solving method
//! 2. Applying symbolic transformations to isolate the target variable
//! 3. Simplifying and evaluating the result
//! 4. Recording all steps in a [`ResolutionPath`] for display
//!
//! # Solver Types
//!
//! - [`LinearSolver`]: Solves equations of the form `ax + b = c`
//! - [`QuadraticSolver`]: Solves equations with x² terms (not yet implemented)
//! - [`PolynomialSolver`]: General polynomial equations (not yet implemented)
//! - [`TranscendentalSolver`]: Equations with sin, cos, tan, exp, ln, log functions
//! - [`SmartSolver`]: Automatically selects the appropriate solver
//!
//! # Solution Types
//!
//! Solutions can be:
//! - [`Solution::Unique`]: Single solution (e.g., x = 5)
//! - [`Solution::Multiple`]: Discrete solutions (e.g., x = 2 or x = -2)
//! - [`Solution::Parametric`]: Solution depends on other variables
//! - [`Solution::None`]: No solution exists (inconsistent equation)
//! - [`Solution::Infinite`]: All values satisfy the equation (identity)
//!
//! # Examples
//!
//! ## Basic Linear Equation
//!
//! ```
//! use thales::solver::{LinearSolver, Solver};
//! use thales::ast::{Equation, Expression, Variable, BinaryOp};
//!
//! // Solve: 2x + 3 = 11
//! let x = Expression::Variable(Variable::new("x"));
//! let left = Expression::Binary(
//!     BinaryOp::Add,
//!     Box::new(Expression::Binary(
//!         BinaryOp::Mul,
//!         Box::new(Expression::Integer(2)),
//!         Box::new(x),
//!     )),
//!     Box::new(Expression::Integer(3)),
//! );
//! let right = Expression::Integer(11);
//! let equation = Equation::new("linear_eq", left, right);
//!
//! let solver = LinearSolver::new();
//! let (solution, path) = solver.solve(&equation, &Variable::new("x")).unwrap();
//!
//! // Solution is x = 4
//! # use thales::solver::Solution;
//! # match solution {
//! #     Solution::Unique(expr) => {
//! #         assert_eq!(expr.evaluate(&std::collections::HashMap::new()), Some(4.0));
//! #     }
//! #     _ => panic!("Expected unique solution"),
//! # }
//! ```
//!
//! ## Using SmartSolver
//!
//! ```
//! use thales::solver::{SmartSolver, Solver};
//! use thales::ast::{Equation, Expression, Variable, BinaryOp};
//!
//! // SmartSolver automatically picks the right method
//! let solver = SmartSolver::new();
//!
//! // Solve: 3x = 12
//! let x = Expression::Variable(Variable::new("x"));
//! let left = Expression::Binary(
//!     BinaryOp::Mul,
//!     Box::new(Expression::Integer(3)),
//!     Box::new(x),
//! );
//! let equation = Equation::new("simple", left, Expression::Integer(12));
//!
//! let (solution, _path) = solver.solve(&equation, &Variable::new("x")).unwrap();
//! // Solution is x = 4
//! ```
//!
//! ## High-Level API with Known Values
//!
//! ```
//! use thales::solver::solve_for;
//! use thales::ast::{Equation, Expression, Variable, BinaryOp};
//! use std::collections::HashMap;
//!
//! // Solve: ax + b = c for x, given a=2, b=3, c=11
//! let a = Expression::Variable(Variable::new("a"));
//! let x = Expression::Variable(Variable::new("x"));
//! let b = Expression::Variable(Variable::new("b"));
//! let c = Expression::Variable(Variable::new("c"));
//!
//! let ax = Expression::Binary(BinaryOp::Mul, Box::new(a), Box::new(x));
//! let left = Expression::Binary(BinaryOp::Add, Box::new(ax), Box::new(b));
//! let equation = Equation::new("parametric", left, c);
//!
//! let mut known = HashMap::new();
//! known.insert("a".to_string(), 2.0);
//! known.insert("b".to_string(), 3.0);
//! known.insert("c".to_string(), 11.0);
//!
//! let path = solve_for(&equation, "x", &known).unwrap();
//! // Result is x = 4.0
//! # assert_eq!(path.result.evaluate(&HashMap::new()), Some(4.0));
//! ```

use crate::ast::{BinaryOp, Equation, Expression, UnaryOp, Variable};
use crate::resolution_path::{Operation, ResolutionPath, ResolutionPathBuilder, ResolutionStep};
use std::collections::HashMap;

/// Error types for equation solving.
///
/// These errors represent different failure modes when attempting to solve
/// an equation symbolically.
///
/// # Examples
///
/// ```
/// use thales::solver::{SolverError, LinearSolver, Solver};
/// use thales::ast::{Equation, Expression, Variable, BinaryOp};
///
/// // NoSolution: 0 = 5 (inconsistent equation)
/// let eq = Equation::new("bad", Expression::Integer(0), Expression::Integer(5));
/// // This would fail with an error during solving
///
/// // CannotSolve: x² = 4 (not linear, LinearSolver can't handle it)
/// let x = Expression::Variable(Variable::new("x"));
/// let x_squared = Expression::Power(
///     Box::new(x),
///     Box::new(Expression::Integer(2)),
/// );
/// let eq = Equation::new("quadratic", x_squared, Expression::Integer(4));
/// let solver = LinearSolver::new();
/// let result = solver.solve(&eq, &Variable::new("x"));
/// assert!(result.is_err());
/// ```
#[derive(Debug, Clone, PartialEq)]
#[non_exhaustive]
pub enum SolverError {
    /// Equation has no solution (inconsistent).
    ///
    /// Example: `0 = 5` or `x + 1 = x + 2`
    NoSolution,

    /// Equation has infinite solutions (identity).
    ///
    /// Example: `x = x` or `2(x + 1) = 2x + 2`
    InfiniteSolutions,

    /// Cannot solve for the given variable with this solver.
    ///
    /// This typically means the equation is too complex for the solver,
    /// or the variable doesn't appear in a solvable form. The message
    /// provides specific details about why solving failed.
    ///
    /// Example: Variable not in equation, or pattern not recognized
    CannotSolve(String),

    /// Equation type is not supported by this solver.
    ///
    /// Example: Trying to solve a quadratic equation with LinearSolver
    UnsupportedEquationType,

    /// Division by zero encountered during solving.
    ///
    /// Example: Attempting to divide by a coefficient that evaluates to zero
    DivisionByZero,

    /// Other error with description.
    ///
    /// Used for errors that don't fit other categories, such as
    /// domain errors (e.g., asin(2)) or not-yet-implemented features.
    Other(String),
}

impl std::fmt::Display for SolverError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            SolverError::NoSolution => write!(f, "Equation has no solution"),
            SolverError::InfiniteSolutions => write!(f, "Equation has infinite solutions"),
            SolverError::CannotSolve(msg) => write!(f, "Cannot solve: {}", msg),
            SolverError::UnsupportedEquationType => write!(f, "Equation type is not supported"),
            SolverError::DivisionByZero => write!(f, "Division by zero encountered"),
            SolverError::Other(msg) => write!(f, "{}", msg),
        }
    }
}

impl std::error::Error for SolverError {}

/// Result type for solver operations.
pub type SolverResult<T> = Result<T, SolverError>;

/// Solution to an equation.
///
/// Represents the different types of solutions an equation can have.
/// Each variant captures a different solution structure.
///
/// # Examples
///
/// ```
/// use thales::solver::{Solution, LinearSolver, Solver};
/// use thales::ast::{Equation, Expression, Variable, BinaryOp};
///
/// // Unique solution: 2x = 8 → x = 4
/// let x = Expression::Variable(Variable::new("x"));
/// let left = Expression::Binary(
///     BinaryOp::Mul,
///     Box::new(Expression::Integer(2)),
///     Box::new(x),
/// );
/// let eq = Equation::new("simple", left, Expression::Integer(8));
///
/// let solver = LinearSolver::new();
/// let (solution, _) = solver.solve(&eq, &Variable::new("x")).unwrap();
///
/// match solution {
///     Solution::Unique(expr) => {
///         // expr evaluates to 4
///         assert_eq!(expr.evaluate(&std::collections::HashMap::new()), Some(4.0));
///     }
///     _ => panic!("Expected unique solution"),
/// }
/// ```
#[derive(Debug, Clone, PartialEq)]
pub enum Solution {
    /// Single unique solution.
    ///
    /// The equation has exactly one solution, represented as an expression.
    ///
    /// # Examples
    ///
    /// - Linear: `2x + 3 = 11` → `x = 4`
    /// - Transcendental: `sin(x) = 0.5` → `x = asin(0.5)`
    Unique(Expression),

    /// Multiple discrete solutions.
    ///
    /// The equation has a finite number of distinct solutions.
    ///
    /// # Examples
    ///
    /// - Quadratic: `x² - 4 = 0` → `x = 2` or `x = -2`
    /// - Trigonometric: `sin(x) = 0` on [0, 2π] → `x = 0, π, 2π`
    Multiple(Vec<Expression>),

    /// Parametric solution with constraints.
    ///
    /// The solution depends on other variables, with optional constraints.
    /// Useful for underdetermined systems or equations with parameters.
    ///
    /// # Examples
    ///
    /// - `x + y = 5` solving for x → `x = 5 - y` (y is a parameter)
    /// - `sqrt(x) = 2` → `x = 4` with constraint `x ≥ 0`
    Parametric {
        /// The solution expression, potentially containing other variables
        expression: Expression,
        /// Constraints that must be satisfied
        constraints: Vec<Constraint>,
    },

    /// No solution exists.
    ///
    /// The equation is inconsistent and has no values that satisfy it.
    ///
    /// # Examples
    ///
    /// - `0 = 5` (contradiction)
    /// - `x + 1 = x + 2` (no solution)
    None,

    /// Infinite solutions (identity).
    ///
    /// The equation is satisfied by all values (tautology).
    ///
    /// # Examples
    ///
    /// - `x = x` (trivial identity)
    /// - `2(x + 1) = 2x + 2` (identity after simplification)
    Infinite,
}

/// Constraint on a solution.
///
/// Represents a condition that must be satisfied for a solution to be valid.
/// Typically used with parametric solutions to specify domain restrictions.
///
/// # Examples
///
/// ```
/// use thales::ast::{Variable, Expression};
/// use thales::solver::Constraint;
///
/// // Constraint: x != 0 (for denominators)
/// // Note: The condition expression format depends on application needs
/// let constraint = Constraint {
///     variable: Variable::new("x"),
///     condition: Expression::Variable(Variable::new("x")),  // Placeholder for non-zero condition
/// };
/// ```
///
/// # Note
///
/// The exact representation of constraints is application-specific. Common uses include:
/// - Domain restrictions (e.g., x > 0 for sqrt, log)
/// - Non-zero denominators
/// - Parameter ranges
#[derive(Debug, Clone, PartialEq)]
pub struct Constraint {
    /// The variable being constrained
    pub variable: Variable,
    /// The condition that must hold (e.g., x >= 0)
    pub condition: Expression,
}

/// Trait for equation solvers.
///
/// Implementors of this trait provide methods to solve equations symbolically.
/// Each solver specializes in a particular type of equation (linear, quadratic, etc.).
///
/// # Design
///
/// The trait has two methods:
/// - [`can_solve`](Solver::can_solve): Quick check if equation is suitable for this solver
/// - [`solve`](Solver::solve): Perform the actual solving and return solution with steps
///
/// This design allows for solver selection (see [`SmartSolver`]) and error handling.
///
/// # Examples
///
/// ```
/// use thales::solver::{Solver, LinearSolver};
/// use thales::ast::{Equation, Expression, Variable, BinaryOp};
///
/// let solver = LinearSolver::new();
///
/// // Build equation: 5x = 20
/// let x = Expression::Variable(Variable::new("x"));
/// let left = Expression::Binary(
///     BinaryOp::Mul,
///     Box::new(Expression::Integer(5)),
///     Box::new(x),
/// );
/// let eq = Equation::new("test", left, Expression::Integer(20));
///
/// // Check if solver can handle it
/// assert!(solver.can_solve(&eq));
///
/// // Solve it
/// let (solution, path) = solver.solve(&eq, &Variable::new("x")).unwrap();
/// // Solution is x = 4
/// ```
pub trait Solver {
    /// Solve an equation for the specified variable.
    ///
    /// Returns the solution(s) and a [`ResolutionPath`] showing the steps taken.
    ///
    /// # Arguments
    ///
    /// * `equation` - The equation to solve
    /// * `variable` - The variable to solve for
    ///
    /// # Returns
    ///
    /// A tuple containing:
    /// - [`Solution`]: The solution (unique, multiple, none, etc.)
    /// - [`ResolutionPath`]: Step-by-step record of solving process
    ///
    /// # Errors
    ///
    /// Returns [`SolverError`] if:
    /// - Variable not found in equation
    /// - Equation type not supported by this solver
    /// - Equation has no solution or infinite solutions
    /// - Other solving failures (see [`SolverError`] variants)
    fn solve(
        &self,
        equation: &Equation,
        variable: &Variable,
    ) -> SolverResult<(Solution, ResolutionPath)>;

    /// Check if this solver can handle the given equation.
    ///
    /// This is a fast pre-check that examines the equation structure without
    /// actually solving it. It's used by [`SmartSolver`] to select the
    /// appropriate solver.
    ///
    /// # Note
    ///
    /// Returning `true` doesn't guarantee successful solving - the equation
    /// might still have no solution or be too complex. This method only
    /// checks if the equation type matches this solver's capabilities.
    ///
    /// # Examples
    ///
    /// ```
    /// use thales::solver::{Solver, LinearSolver};
    /// use thales::ast::{Equation, Expression, Variable};
    ///
    /// let solver = LinearSolver::new();
    ///
    /// // Linear equation: x + 5 = 10
    /// let eq1 = Equation::new(
    ///     "linear",
    ///     Expression::Variable(Variable::new("x")),
    ///     Expression::Integer(5),
    /// );
    /// assert!(solver.can_solve(&eq1));
    ///
    /// // Quadratic equation: x² = 4
    /// let x = Expression::Variable(Variable::new("x"));
    /// let x_squared = Expression::Power(
    ///     Box::new(x),
    ///     Box::new(Expression::Integer(2)),
    /// );
    /// let eq2 = Equation::new("quadratic", x_squared, Expression::Integer(4));
    /// assert!(!solver.can_solve(&eq2)); // LinearSolver rejects quadratics
    /// ```
    fn can_solve(&self, equation: &Equation) -> bool;
}

// ============================================================================
// Helper Functions
// ============================================================================

/// Check if expression contains the given variable.
///
/// This is a convenience wrapper around [`Expression::contains_variable`].
///
/// # Examples
///
/// ```ignore
/// let x = Expression::Variable(Variable::new("x"));
/// let y = Expression::Variable(Variable::new("y"));
/// let expr = Expression::Binary(BinaryOp::Add, Box::new(x), Box::new(Expression::Integer(5)));
///
/// assert!(contains_variable(&expr, "x"));
/// assert!(!contains_variable(&expr, "y"));
/// ```
fn contains_variable(expr: &Expression, var: &str) -> bool {
    expr.contains_variable(var)
}

/// Extract the coefficient of a variable from an expression.
///
/// Recognizes patterns like:
/// - `x` → coefficient is 1
/// - `3 * x` → coefficient is 3
/// - `x * 3` → coefficient is 3
/// - `a * x` → coefficient is a (where a doesn't contain x)
///
/// Returns `None` if the variable is not found or appears in a non-linear way.
///
/// # Examples
///
/// ```ignore
/// // 3 * x → Some(3)
/// let expr = Expression::Binary(
///     BinaryOp::Mul,
///     Box::new(Expression::Integer(3)),
///     Box::new(Expression::Variable(Variable::new("x"))),
/// );
/// let coeff = extract_coefficient(&expr, "x");
/// assert_eq!(coeff, Some(Expression::Integer(3)));
///
/// // x → Some(1)
/// let expr = Expression::Variable(Variable::new("x"));
/// let coeff = extract_coefficient(&expr, "x");
/// assert_eq!(coeff, Some(Expression::Integer(1)));
///
/// // x + 5 → None (not a pure coefficient pattern)
/// let expr = Expression::Binary(
///     BinaryOp::Add,
///     Box::new(Expression::Variable(Variable::new("x"))),
///     Box::new(Expression::Integer(5)),
/// );
/// let coeff = extract_coefficient(&expr, "x");
/// assert_eq!(coeff, None);
/// ```
fn extract_coefficient(expr: &Expression, var: &str) -> Option<Expression> {
    match expr {
        // x -> coefficient is 1
        Expression::Variable(v) if v.name == var => Some(Expression::Integer(1)),

        // a * x or x * a
        Expression::Binary(BinaryOp::Mul, left, right) => {
            if let Expression::Variable(v) = left.as_ref() {
                if v.name == var && !contains_variable(right, var) {
                    return Some(right.as_ref().clone());
                }
            }
            if let Expression::Variable(v) = right.as_ref() {
                if v.name == var && !contains_variable(left, var) {
                    return Some(left.as_ref().clone());
                }
            }
            None
        }

        _ => None,
    }
}

/// Evaluate constant expressions to their numeric values.
/// If the expression contains only constants, evaluate it completely.
fn evaluate_constants(expr: &Expression) -> Expression {
    // First simplify
    let simplified = expr.simplify();

    // Try to evaluate if it's all constants
    if !has_any_variable(&simplified) {
        if let Some(value) = simplified.evaluate(&HashMap::new()) {
            // Check if it's an integer value
            if value.fract().abs() < 1e-10 {
                return Expression::Integer(value.round() as i64);
            } else {
                return Expression::Float(value);
            }
        }
    }

    simplified
}

/// Isolate a variable in an equation.
///
/// Rearranges the equation to solve for the target variable, returning the
/// expression that equals the variable. This is the core solving logic for
/// linear equations.
///
/// # Algorithm
///
/// Recognizes and solves several linear patterns:
/// 1. Variable already isolated: `x = expr` or `expr = x`
/// 2. Coefficient pattern: `a*x = c` → `x = c/a`
/// 3. Addition pattern: `x + b = c` → `x = c - b`
/// 4. Combined pattern: `a*x + b = c` → `x = (c - b)/a`
///
/// All patterns are checked in both left-to-right and right-to-left orientations.
///
/// # Arguments
///
/// * `equation` - The equation to solve
/// * `var` - The variable name to isolate
/// * `path` - Resolution path builder to record solving steps
///
/// # Returns
///
/// An [`Expression`] representing the isolated variable's value.
///
/// # Errors
///
/// Returns [`SolverError`] if:
/// - Variable not found in equation
/// - Equation pattern not recognized (too complex for Phase 1)
///
/// # Examples
///
/// ```ignore
/// // 2x + 3 = 11
/// // Step 1: Recognize pattern a*x + b = c
/// // Step 2: Compute (c - b) / a = (11 - 3) / 2 = 4
/// // Returns: Expression::Integer(4)
///
/// let equation = /* ... */;
/// let mut path = ResolutionPathBuilder::new(/* ... */);
/// let result = isolate_variable(&equation, "x", &mut path)?;
/// ```
fn isolate_variable(
    equation: &Equation,
    var: &str,
    _path: &mut ResolutionPathBuilder,
) -> Result<Expression, SolverError> {
    let left = &equation.left;
    let right = &equation.right;

    // Check if variable exists in equation
    if !contains_variable(left, var) && !contains_variable(right, var) {
        return Err(SolverError::CannotSolve(format!(
            "Variable '{}' not found in equation",
            var
        )));
    }

    // Special case: variable already isolated (x = expr or expr = x)
    if let Expression::Variable(v) = left {
        if v.name == var && !contains_variable(right, var) {
            return Ok(right.clone());
        }
    }
    if let Expression::Variable(v) = right {
        if v.name == var && !contains_variable(left, var) {
            return Ok(left.clone());
        }
    }

    // Try to solve simple patterns

    // Pattern: a * x = c  =>  x = c / a
    if let Some(coeff) = extract_coefficient(left, var) {
        if !contains_variable(right, var) {
            let result = Expression::Binary(
                BinaryOp::Div,
                Box::new(right.clone()),
                Box::new(coeff.clone()),
            )
            .simplify();
            let evaluated = evaluate_constants(&result);
            return Ok(evaluated);
        }
    }

    // Pattern: c = a * x  =>  x = c / a
    if let Some(coeff) = extract_coefficient(right, var) {
        if !contains_variable(left, var) {
            let result = Expression::Binary(
                BinaryOp::Div,
                Box::new(left.clone()),
                Box::new(coeff.clone()),
            )
            .simplify();
            let evaluated = evaluate_constants(&result);
            return Ok(evaluated);
        }
    }

    // Pattern: x + b = c  =>  x = c - b
    if let Expression::Binary(BinaryOp::Add, l, r) = left {
        if let Expression::Variable(v) = l.as_ref() {
            if v.name == var && !contains_variable(r, var) && !contains_variable(right, var) {
                let result = Expression::Binary(
                    BinaryOp::Sub,
                    Box::new(right.clone()),
                    Box::new(r.as_ref().clone()),
                )
                .simplify();
                let evaluated = evaluate_constants(&result);
                return Ok(evaluated);
            }
        }
        if let Expression::Variable(v) = r.as_ref() {
            if v.name == var && !contains_variable(l, var) && !contains_variable(right, var) {
                let result = Expression::Binary(
                    BinaryOp::Sub,
                    Box::new(right.clone()),
                    Box::new(l.as_ref().clone()),
                )
                .simplify();
                let evaluated = evaluate_constants(&result);
                return Ok(evaluated);
            }
        }
    }

    // Pattern: c = x + b  =>  x = c - b
    if let Expression::Binary(BinaryOp::Add, l, r) = right {
        if let Expression::Variable(v) = l.as_ref() {
            if v.name == var && !contains_variable(r, var) && !contains_variable(left, var) {
                let result = Expression::Binary(
                    BinaryOp::Sub,
                    Box::new(left.clone()),
                    Box::new(r.as_ref().clone()),
                )
                .simplify();
                let evaluated = evaluate_constants(&result);
                return Ok(evaluated);
            }
        }
        if let Expression::Variable(v) = r.as_ref() {
            if v.name == var && !contains_variable(l, var) && !contains_variable(left, var) {
                let result = Expression::Binary(
                    BinaryOp::Sub,
                    Box::new(left.clone()),
                    Box::new(l.as_ref().clone()),
                )
                .simplify();
                let evaluated = evaluate_constants(&result);
                return Ok(evaluated);
            }
        }
    }

    // Pattern: a * x + b = c  =>  x = (c - b) / a
    if let Expression::Binary(BinaryOp::Add, l, r) = left {
        if let Some(coeff) = extract_coefficient(l, var) {
            if !contains_variable(r, var) && !contains_variable(right, var) {
                let numerator = Expression::Binary(
                    BinaryOp::Sub,
                    Box::new(right.clone()),
                    Box::new(r.as_ref().clone()),
                );
                let result =
                    Expression::Binary(BinaryOp::Div, Box::new(numerator), Box::new(coeff))
                        .simplify();
                let evaluated = evaluate_constants(&result);
                return Ok(evaluated);
            }
        }
        if let Some(coeff) = extract_coefficient(r, var) {
            if !contains_variable(l, var) && !contains_variable(right, var) {
                let numerator = Expression::Binary(
                    BinaryOp::Sub,
                    Box::new(right.clone()),
                    Box::new(l.as_ref().clone()),
                );
                let result =
                    Expression::Binary(BinaryOp::Div, Box::new(numerator), Box::new(coeff))
                        .simplify();
                let evaluated = evaluate_constants(&result);
                return Ok(evaluated);
            }
        }
    }

    // More complex cases not yet supported
    Err(SolverError::CannotSolve(
        "Equation pattern not yet supported for Phase 1".to_string(),
    ))
}

/// Linear equation solver for equations of the form `ax + b = c`.
///
/// Solves first-degree polynomial equations in one variable by pattern matching
/// and algebraic manipulation. Handles various linear patterns including:
/// - Simple variable: `x = 5`
/// - Multiplication: `3x = 12`
/// - Addition: `x + 7 = 10`
/// - Combined: `2x + 3 = 11`
///
/// # Mathematical Foundation
///
/// A linear equation in one variable has the general form:
/// ```text
/// ax + b = c
/// ```
///
/// The solution is obtained by:
/// 1. Subtracting `b` from both sides: `ax = c - b`
/// 2. Dividing both sides by `a`: `x = (c - b) / a`
///
/// The solver recognizes these patterns automatically and applies the
/// appropriate transformations.
///
/// # Limitations
///
/// - Only handles linear equations (degree 1)
/// - Cannot solve equations with the variable in denominators (e.g., `1/x = 2`)
/// - Cannot solve equations with the variable in exponents (e.g., `2^x = 8`)
/// - Cannot handle products of variables (e.g., `x*y = 5`)
///
/// For more complex equations, use [`TranscendentalSolver`] or [`SmartSolver`].
///
/// # Examples
///
/// ## Simple Linear Equation
///
/// ```
/// use thales::solver::{LinearSolver, Solver};
/// use thales::ast::{Equation, Expression, Variable, BinaryOp};
///
/// // Solve: 2x + 3 = 11
/// let x = Expression::Variable(Variable::new("x"));
/// let two_x = Expression::Binary(
///     BinaryOp::Mul,
///     Box::new(Expression::Integer(2)),
///     Box::new(x),
/// );
/// let left = Expression::Binary(
///     BinaryOp::Add,
///     Box::new(two_x),
///     Box::new(Expression::Integer(3)),
/// );
/// let equation = Equation::new("linear", left, Expression::Integer(11));
///
/// let solver = LinearSolver::new();
/// let (solution, path) = solver.solve(&equation, &Variable::new("x")).unwrap();
///
/// // Verify solution: x = 4
/// # use thales::solver::Solution;
/// # use std::collections::HashMap;
/// # match solution {
/// #     Solution::Unique(expr) => {
/// #         assert_eq!(expr.evaluate(&HashMap::new()), Some(4.0));
/// #     }
/// #     _ => panic!("Expected unique solution"),
/// # }
/// ```
///
/// ## Equation with Parametric Coefficients
///
/// ```
/// use thales::solver::{LinearSolver, Solver};
/// use thales::ast::{Equation, Expression, Variable, BinaryOp};
/// use std::collections::HashMap;
///
/// // Solve: ax = b for x (symbolic)
/// let a = Expression::Variable(Variable::new("a"));
/// let x = Expression::Variable(Variable::new("x"));
/// let b = Expression::Variable(Variable::new("b"));
///
/// let left = Expression::Binary(BinaryOp::Mul, Box::new(a), Box::new(x));
/// let equation = Equation::new("parametric", left, b);
///
/// let solver = LinearSolver::new();
/// let (solution, _path) = solver.solve(&equation, &Variable::new("x")).unwrap();
///
/// // Solution is symbolic: x = b/a
/// # use thales::solver::Solution;
/// # match solution {
/// #     Solution::Unique(expr) => {
/// #         // Can substitute values later
/// #         let mut values = HashMap::new();
/// #         values.insert("a".to_string(), 3.0);
/// #         values.insert("b".to_string(), 12.0);
/// #         // Result would be 4.0
/// #     }
/// #     _ => panic!("Expected unique solution"),
/// # }
/// ```
///
/// ## Checking Solver Applicability
///
/// ```
/// use thales::solver::{LinearSolver, Solver};
/// use thales::ast::{Equation, Expression, Variable};
///
/// let solver = LinearSolver::new();
///
/// // Can solve linear equation
/// let linear = Equation::new(
///     "linear",
///     Expression::Variable(Variable::new("x")),
///     Expression::Integer(5),
/// );
/// assert!(solver.can_solve(&linear));
///
/// // Cannot solve quadratic equation
/// let x = Expression::Variable(Variable::new("x"));
/// let x_squared = Expression::Power(
///     Box::new(x),
///     Box::new(Expression::Integer(2)),
/// );
/// let quadratic = Equation::new("quadratic", x_squared, Expression::Integer(4));
/// assert!(!solver.can_solve(&quadratic));
/// ```
///
/// # See Also
///
/// - [`SmartSolver`]: Automatically selects LinearSolver for linear equations
/// - [`TranscendentalSolver`]: For equations with sin, cos, exp, ln, etc.
/// - [`solve_for`]: High-level API that uses SmartSolver and substitutes known values
#[derive(Debug, Default)]
pub struct LinearSolver;

impl LinearSolver {
    /// Create a new linear equation solver.
    ///
    /// # Examples
    ///
    /// ```
    /// use thales::solver::LinearSolver;
    ///
    /// let solver = LinearSolver::new();
    /// ```
    pub fn new() -> Self {
        Self
    }
}

impl Solver for LinearSolver {
    fn solve(
        &self,
        equation: &Equation,
        variable: &Variable,
    ) -> SolverResult<(Solution, ResolutionPath)> {
        let var_name = &variable.name;

        // Initialize resolution path
        let initial_expr = Expression::Binary(
            BinaryOp::Sub,
            Box::new(equation.left.clone()),
            Box::new(equation.right.clone()),
        );
        let mut path = ResolutionPathBuilder::new(initial_expr.clone());

        // Check if variable appears in equation
        let left_has_var = contains_variable(&equation.left, var_name);
        let right_has_var = contains_variable(&equation.right, var_name);

        if !left_has_var && !right_has_var {
            return Err(SolverError::CannotSolve(format!(
                "Variable '{}' not found in equation",
                var_name
            )));
        }

        // Check if equation is linear in the target variable
        if !is_linear_in_variable(&equation.left, var_name)
            || !is_linear_in_variable(&equation.right, var_name)
        {
            return Err(SolverError::UnsupportedEquationType);
        }

        // Isolate the variable
        let result_expr = isolate_variable(equation, var_name, &mut path)?;

        // Add isolation step
        path = path.step(
            Operation::Isolate(variable.clone()),
            format!("Isolate {} on one side", variable),
            result_expr.clone(),
        );

        // Build final resolution path
        let resolution_path = path.finish(result_expr.clone());

        Ok((Solution::Unique(result_expr), resolution_path))
    }

    fn can_solve(&self, equation: &Equation) -> bool {
        // Check if equation has obvious non-linear features (powers > 1 with variables)
        // We're more permissive here since we don't know the target variable yet,
        // but we can still reject clearly quadratic/polynomial equations.
        !has_obvious_nonlinearity(&equation.left) && !has_obvious_nonlinearity(&equation.right)
    }
}

/// Check if an expression has obvious non-linear features like x^2.
fn has_obvious_nonlinearity(expr: &Expression) -> bool {
    match expr {
        Expression::Power(base, exp) => {
            // x^2 or any variable raised to power > 1
            if has_any_variable(base) {
                // Check if exponent is > 1
                if let Some(exp_val) = exp.evaluate(&HashMap::new()) {
                    if exp_val > 1.0 {
                        return true;
                    }
                }
            }
            has_obvious_nonlinearity(base) || has_obvious_nonlinearity(exp)
        }
        Expression::Unary(_, inner) => has_obvious_nonlinearity(inner),
        Expression::Binary(_, left, right) => {
            has_obvious_nonlinearity(left) || has_obvious_nonlinearity(right)
        }
        Expression::Function(_, args) => args.iter().any(|arg| has_obvious_nonlinearity(arg)),
        _ => false,
    }
}

/// Check if expression contains any variables.
fn has_any_variable(expr: &Expression) -> bool {
    match expr {
        Expression::Variable(_) => true,
        Expression::Unary(_, inner) => has_any_variable(inner),
        Expression::Binary(_, left, right) => has_any_variable(left) || has_any_variable(right),
        Expression::Function(_, args) => args.iter().any(has_any_variable),
        Expression::Power(base, exp) => has_any_variable(base) || has_any_variable(exp),
        _ => false,
    }
}

/// Extract coefficients (a, b, c) from a quadratic expression ax² + bx + c.
fn extract_quadratic_coefficients(expr: &Expression, var: &str) -> (f64, f64, f64) {
    let mut a = 0.0;
    let mut b = 0.0;
    let mut c = 0.0;

    extract_poly_coefficients_recursive(expr, var, 1.0, &mut a, &mut b, &mut c);
    (a, b, c)
}

/// Recursively extract polynomial coefficients.
fn extract_poly_coefficients_recursive(
    expr: &Expression,
    var: &str,
    multiplier: f64,
    a: &mut f64,
    b: &mut f64,
    c: &mut f64,
) {
    match expr {
        Expression::Integer(n) => *c += (*n as f64) * multiplier,
        Expression::Float(f) => *c += f * multiplier,
        Expression::Rational(r) => *c += (*r.numer() as f64 / *r.denom() as f64) * multiplier,
        Expression::Variable(v) if v.name == var => *b += multiplier,
        Expression::Variable(_) | Expression::Constant(_) => {
            // Other variable or constant treated as part of constant term
            if let Some(val) = expr.evaluate(&std::collections::HashMap::new()) {
                *c += val * multiplier;
            }
        }
        Expression::Unary(UnaryOp::Neg, inner) => {
            extract_poly_coefficients_recursive(inner, var, -multiplier, a, b, c);
        }
        Expression::Binary(BinaryOp::Add, left, right) => {
            extract_poly_coefficients_recursive(left, var, multiplier, a, b, c);
            extract_poly_coefficients_recursive(right, var, multiplier, a, b, c);
        }
        Expression::Binary(BinaryOp::Sub, left, right) => {
            extract_poly_coefficients_recursive(left, var, multiplier, a, b, c);
            extract_poly_coefficients_recursive(right, var, -multiplier, a, b, c);
        }
        Expression::Binary(BinaryOp::Mul, left, right) => {
            // Check for coefficient * variable or coefficient * x^2
            let left_val = left.evaluate(&std::collections::HashMap::new());
            let right_val = right.evaluate(&std::collections::HashMap::new());

            match (left_val, right_val) {
                (Some(lv), None) => {
                    extract_poly_coefficients_recursive(right, var, multiplier * lv, a, b, c);
                }
                (None, Some(rv)) => {
                    extract_poly_coefficients_recursive(left, var, multiplier * rv, a, b, c);
                }
                (Some(lv), Some(rv)) => {
                    *c += lv * rv * multiplier;
                }
                (None, None) => {
                    // Could be x * x = x^2 or other variable products
                    if matches!(&**left, Expression::Variable(v) if v.name == var)
                        && matches!(&**right, Expression::Variable(v) if v.name == var)
                    {
                        *a += multiplier;
                    } else if matches!(&**left, Expression::Variable(v) if v.name == var) {
                        if let Some(rv) = right.evaluate(&std::collections::HashMap::new()) {
                            *b += multiplier * rv;
                        }
                    } else if matches!(&**right, Expression::Variable(v) if v.name == var) {
                        if let Some(lv) = left.evaluate(&std::collections::HashMap::new()) {
                            *b += multiplier * lv;
                        }
                    }
                }
            }
        }
        Expression::Power(base, exp) => {
            // Check for x^2 or x^n
            if matches!(&**base, Expression::Variable(v) if v.name == var) {
                if let Some(exp_val) = exp.evaluate(&std::collections::HashMap::new()) {
                    if (exp_val - 2.0).abs() < 1e-10 {
                        *a += multiplier;
                    } else if (exp_val - 1.0).abs() < 1e-10 {
                        *b += multiplier;
                    } else if exp_val.abs() < 1e-10 {
                        *c += multiplier;
                    }
                }
            }
        }
        _ => {
            // For other expressions, try to evaluate as constant
            if let Some(val) = expr.evaluate(&std::collections::HashMap::new()) {
                *c += val * multiplier;
            }
        }
    }
}

/// Simplify a numeric value to the best Expression representation.
fn simplify_numeric_expression(val: f64) -> Expression {
    // Check if it's close to an integer
    let rounded = val.round();
    if (val - rounded).abs() < 1e-10 && rounded.abs() < i64::MAX as f64 {
        Expression::Integer(rounded as i64)
    } else {
        Expression::Float(val)
    }
}

/// Extract coefficients for a general polynomial.
/// Returns vector of coefficients [a0, a1, a2, ..., an] for a0 + a1*x + a2*x^2 + ...
fn extract_polynomial_coefficients(expr: &Expression, var: &str, max_degree: usize) -> Vec<f64> {
    let mut coeffs = vec![0.0; max_degree + 1];
    extract_general_poly_coefficients(expr, var, 1.0, &mut coeffs);
    coeffs
}

/// Recursively extract general polynomial coefficients.
fn extract_general_poly_coefficients(
    expr: &Expression,
    var: &str,
    multiplier: f64,
    coeffs: &mut [f64],
) {
    match expr {
        Expression::Integer(n) => coeffs[0] += (*n as f64) * multiplier,
        Expression::Float(f) => coeffs[0] += f * multiplier,
        Expression::Rational(r) => {
            coeffs[0] += (*r.numer() as f64 / *r.denom() as f64) * multiplier
        }
        Expression::Variable(v) if v.name == var => {
            if coeffs.len() > 1 {
                coeffs[1] += multiplier;
            }
        }
        Expression::Variable(_) | Expression::Constant(_) => {
            if let Some(val) = expr.evaluate(&std::collections::HashMap::new()) {
                coeffs[0] += val * multiplier;
            }
        }
        Expression::Unary(UnaryOp::Neg, inner) => {
            extract_general_poly_coefficients(inner, var, -multiplier, coeffs);
        }
        Expression::Binary(BinaryOp::Add, left, right) => {
            extract_general_poly_coefficients(left, var, multiplier, coeffs);
            extract_general_poly_coefficients(right, var, multiplier, coeffs);
        }
        Expression::Binary(BinaryOp::Sub, left, right) => {
            extract_general_poly_coefficients(left, var, multiplier, coeffs);
            extract_general_poly_coefficients(right, var, -multiplier, coeffs);
        }
        Expression::Binary(BinaryOp::Mul, left, right) => {
            let left_val = left.evaluate(&std::collections::HashMap::new());
            let right_val = right.evaluate(&std::collections::HashMap::new());

            match (left_val, right_val) {
                (Some(lv), None) => {
                    extract_general_poly_coefficients(right, var, multiplier * lv, coeffs);
                }
                (None, Some(rv)) => {
                    extract_general_poly_coefficients(left, var, multiplier * rv, coeffs);
                }
                (Some(lv), Some(rv)) => {
                    coeffs[0] += lv * rv * multiplier;
                }
                (None, None) => {
                    // Variable * variable case
                    if matches!(&**left, Expression::Variable(v) if v.name == var)
                        && matches!(&**right, Expression::Variable(v) if v.name == var)
                    {
                        if coeffs.len() > 2 {
                            coeffs[2] += multiplier;
                        }
                    }
                }
            }
        }
        Expression::Power(base, exp) => {
            if matches!(&**base, Expression::Variable(v) if v.name == var) {
                if let Some(exp_val) = exp.evaluate(&std::collections::HashMap::new()) {
                    let degree = exp_val.round() as usize;
                    if degree < coeffs.len() {
                        coeffs[degree] += multiplier;
                    }
                }
            }
        }
        _ => {
            if let Some(val) = expr.evaluate(&std::collections::HashMap::new()) {
                coeffs[0] += val * multiplier;
            }
        }
    }
}

/// Get the degree of a polynomial expression with respect to a variable.
fn get_polynomial_degree(expr: &Expression, var: &str) -> usize {
    match expr {
        Expression::Integer(_)
        | Expression::Rational(_)
        | Expression::Float(_)
        | Expression::Complex(_)
        | Expression::Constant(_) => 0,

        Expression::Variable(v) if v.name == var => 1,
        Expression::Variable(_) => 0,

        Expression::Unary(UnaryOp::Neg, inner) => get_polynomial_degree(inner, var),

        Expression::Binary(BinaryOp::Add | BinaryOp::Sub, left, right) => {
            get_polynomial_degree(left, var).max(get_polynomial_degree(right, var))
        }

        Expression::Binary(BinaryOp::Mul, left, right) => {
            get_polynomial_degree(left, var) + get_polynomial_degree(right, var)
        }

        Expression::Binary(BinaryOp::Div, left, right) => {
            // Division by variable increases complexity, treat as special case
            if contains_variable(right, var) {
                0 // Not a polynomial
            } else {
                get_polynomial_degree(left, var)
            }
        }

        Expression::Power(base, exp) => {
            if let Expression::Variable(v) = base.as_ref() {
                if v.name == var {
                    if let Some(exp_val) = exp.evaluate(&HashMap::new()) {
                        if exp_val >= 0.0 && (exp_val - exp_val.round()).abs() < 1e-10 {
                            return exp_val.round() as usize;
                        }
                    }
                }
            }
            // For complex powers, multiply base degree by power
            let base_deg = get_polynomial_degree(base, var);
            if base_deg == 0 {
                0
            } else if let Some(exp_val) = exp.evaluate(&HashMap::new()) {
                if exp_val >= 0.0 && (exp_val - exp_val.round()).abs() < 1e-10 {
                    base_deg * (exp_val.round() as usize)
                } else {
                    0
                }
            } else {
                0
            }
        }

        _ => 0,
    }
}

/// Check if an expression is polynomial (contains no transcendental functions).
fn is_polynomial_expression(expr: &Expression) -> bool {
    match expr {
        Expression::Integer(_)
        | Expression::Rational(_)
        | Expression::Float(_)
        | Expression::Complex(_)
        | Expression::Constant(_)
        | Expression::Variable(_) => true,

        Expression::Unary(_, inner) => is_polynomial_expression(inner),

        Expression::Binary(_, left, right) => {
            is_polynomial_expression(left) && is_polynomial_expression(right)
        }

        Expression::Power(base, exp) => {
            // Power is polynomial if base is polynomial and exponent is a non-negative integer
            if !is_polynomial_expression(base) {
                return false;
            }
            if let Some(exp_val) = exp.evaluate(&HashMap::new()) {
                exp_val >= 0.0 && (exp_val - exp_val.round()).abs() < 1e-10
            } else {
                // If exponent contains variables, check if it's polynomial
                is_polynomial_expression(exp)
            }
        }

        Expression::Function(_, _) => false, // Functions are transcendental
    }
}

/// Solve cubic equation ax³ + bx² + cx + d = 0 using Cardano's formula.
/// coeffs = [d, c, b, a] (constant term first)
fn solve_cubic(
    coeffs: &[f64],
    _var: &str,
    mut path: ResolutionPathBuilder,
) -> SolverResult<(Solution, ResolutionPath)> {
    if coeffs.len() < 4 {
        return Err(SolverError::CannotSolve(
            "Not a cubic polynomial".to_string(),
        ));
    }

    let d = coeffs[0];
    let c = coeffs[1];
    let b = coeffs[2];
    let a = coeffs[3];

    if a.abs() < 1e-15 {
        // Not actually cubic, delegate to quadratic
        return Err(SolverError::CannotSolve(
            "Leading coefficient is zero".to_string(),
        ));
    }

    // Normalize to monic form: x³ + px² + qx + r = 0
    let p = b / a;
    let q = c / a;
    let r = d / a;

    path = path.step(
        Operation::Simplify,
        format!("Normalized cubic: x³ + {}x² + {}x + {} = 0", p, q, r),
        Expression::Integer(0),
    );

    // Depress the cubic: substitute x = t - p/3
    // t³ + pt + q = 0 where:
    // p = q - p²/3
    // q = r - pq/3 + 2p³/27
    let dep_p = q - p * p / 3.0;
    let dep_q = r - p * q / 3.0 + 2.0 * p * p * p / 27.0;

    path = path.step(
        Operation::Simplify,
        format!("Depressed cubic: t³ + {}t + {} = 0", dep_p, dep_q),
        Expression::Integer(0),
    );

    // Discriminant: Δ = -4p³ - 27q²
    let discriminant = -4.0 * dep_p * dep_p * dep_p - 27.0 * dep_q * dep_q;

    path = path.step(
        Operation::Simplify,
        format!("Discriminant: Δ = {}", discriminant),
        Expression::Integer(0),
    );

    let shift = -p / 3.0;
    let roots: Vec<Expression>;

    if discriminant.abs() < 1e-10 {
        // All roots are real, at least two are equal
        if dep_p.abs() < 1e-10 && dep_q.abs() < 1e-10 {
            // Triple root at t = 0
            let root = simplify_numeric_expression(shift);
            roots = vec![root.clone(), root.clone(), root];
        } else {
            // One single root and one double root
            let t1 = 3.0 * dep_q / dep_p;
            let t2 = -3.0 * dep_q / (2.0 * dep_p);
            roots = vec![
                simplify_numeric_expression(t1 + shift),
                simplify_numeric_expression(t2 + shift),
                simplify_numeric_expression(t2 + shift),
            ];
        }
    } else if discriminant > 0.0 {
        // Three distinct real roots (casus irreducibilis)
        // Use trigonometric method
        let m = 2.0 * (-dep_p / 3.0).sqrt();
        let theta = (3.0 * dep_q / (dep_p * m)).acos() / 3.0;

        let t1 = m * theta.cos();
        let t2 = m * (theta - 2.0 * std::f64::consts::PI / 3.0).cos();
        let t3 = m * (theta - 4.0 * std::f64::consts::PI / 3.0).cos();

        roots = vec![
            simplify_numeric_expression(t1 + shift),
            simplify_numeric_expression(t2 + shift),
            simplify_numeric_expression(t3 + shift),
        ];
    } else {
        // One real root and two complex conjugate roots
        // Use Cardano's formula
        let sqrt_term = (dep_q * dep_q / 4.0 + dep_p * dep_p * dep_p / 27.0).sqrt();
        let u = (-dep_q / 2.0 + sqrt_term).cbrt();
        let v = (-dep_q / 2.0 - sqrt_term).cbrt();

        let t_real = u + v;
        let real_part = -0.5 * (u + v) + shift;
        let imag_part = (3.0_f64).sqrt() / 2.0 * (u - v);

        roots = vec![
            simplify_numeric_expression(t_real + shift),
            Expression::Complex(num_complex::Complex64::new(real_part, imag_part)),
            Expression::Complex(num_complex::Complex64::new(real_part, -imag_part)),
        ];
    }

    path = path.step(
        Operation::Simplify,
        "Applied Cardano's formula".to_string(),
        roots[0].clone(),
    );

    let resolution_path = path.finish(roots[0].clone());
    Ok((Solution::Multiple(roots), resolution_path))
}

/// Solve quartic equation ax⁴ + bx³ + cx² + dx + e = 0 using Ferrari's method.
/// coeffs = [e, d, c, b, a] (constant term first)
fn solve_quartic(
    coeffs: &[f64],
    _var: &str,
    mut path: ResolutionPathBuilder,
) -> SolverResult<(Solution, ResolutionPath)> {
    if coeffs.len() < 5 {
        return Err(SolverError::CannotSolve(
            "Not a quartic polynomial".to_string(),
        ));
    }

    let e = coeffs[0];
    let d = coeffs[1];
    let c = coeffs[2];
    let b = coeffs[3];
    let a = coeffs[4];

    if a.abs() < 1e-15 {
        return Err(SolverError::CannotSolve(
            "Leading coefficient is zero".to_string(),
        ));
    }

    // Normalize to monic form: x⁴ + px³ + qx² + rx + s = 0
    let p = b / a;
    let q = c / a;
    let r = d / a;
    let s = e / a;

    path = path.step(
        Operation::Simplify,
        format!(
            "Normalized quartic: x⁴ + {}x³ + {}x² + {}x + {} = 0",
            p, q, r, s
        ),
        Expression::Integer(0),
    );

    // Depress the quartic: substitute x = y - p/4
    // y⁴ + αy² + βy + γ = 0
    let alpha = q - 3.0 * p * p / 8.0;
    let beta = r - p * q / 2.0 + p * p * p / 8.0;
    let gamma = s - p * r / 4.0 + p * p * q / 16.0 - 3.0 * p * p * p * p / 256.0;

    path = path.step(
        Operation::Simplify,
        format!(
            "Depressed quartic: y⁴ + {}y² + {}y + {} = 0",
            alpha, beta, gamma
        ),
        Expression::Integer(0),
    );

    let shift = -p / 4.0;

    // Handle special case: β = 0 (biquadratic)
    if beta.abs() < 1e-15 {
        // y⁴ + αy² + γ = 0, substitute u = y²
        let disc = alpha * alpha - 4.0 * gamma;
        if disc < -1e-15 {
            // Complex roots
            let u1_real = -alpha / 2.0;
            let u1_imag = (-disc).sqrt() / 2.0;

            // y² = u gives y = ±√u (complex square roots)
            let mut roots = Vec::new();
            for sign1 in [-1.0, 1.0] {
                let u_real = u1_real;
                let u_imag = sign1 * u1_imag;
                // √(a + bi) = ±(√((r+a)/2) + i*sign(b)*√((r-a)/2))
                let r = (u_real * u_real + u_imag * u_imag).sqrt();
                let sqrt_real = ((r + u_real) / 2.0).sqrt();
                let sqrt_imag = u_imag.signum() * ((r - u_real) / 2.0).sqrt();
                roots.push(Expression::Complex(num_complex::Complex64::new(
                    sqrt_real + shift,
                    sqrt_imag,
                )));
                roots.push(Expression::Complex(num_complex::Complex64::new(
                    -sqrt_real + shift,
                    -sqrt_imag,
                )));
            }
            let resolution_path = path.finish(roots[0].clone());
            return Ok((Solution::Multiple(roots), resolution_path));
        } else {
            let u1 = (-alpha + disc.sqrt()) / 2.0;
            let u2 = (-alpha - disc.sqrt()) / 2.0;

            let mut roots = Vec::new();
            for u in [u1, u2] {
                if u >= 0.0 {
                    roots.push(simplify_numeric_expression(u.sqrt() + shift));
                    roots.push(simplify_numeric_expression(-u.sqrt() + shift));
                } else {
                    let imag = (-u).sqrt();
                    roots.push(Expression::Complex(num_complex::Complex64::new(
                        shift, imag,
                    )));
                    roots.push(Expression::Complex(num_complex::Complex64::new(
                        shift, -imag,
                    )));
                }
            }
            let resolution_path = path.finish(roots[0].clone());
            return Ok((Solution::Multiple(roots), resolution_path));
        }
    }

    // Solve resolvent cubic: m³ + (α/2)m² + ((α² - 4γ)/16)m - β²/64 = 0
    let resolvent_coeffs = vec![
        -beta * beta / 64.0,
        (alpha * alpha - 4.0 * gamma) / 16.0,
        alpha / 2.0,
        1.0,
    ];

    // Get one real root of the resolvent cubic
    let dep_p = resolvent_coeffs[1] - resolvent_coeffs[2] * resolvent_coeffs[2] / 3.0;
    let dep_q = resolvent_coeffs[0] - resolvent_coeffs[2] * resolvent_coeffs[1] / 3.0
        + 2.0 * resolvent_coeffs[2] * resolvent_coeffs[2] * resolvent_coeffs[2] / 27.0;

    let disc_cubic = -4.0 * dep_p * dep_p * dep_p - 27.0 * dep_q * dep_q;

    let m: f64;
    if disc_cubic > 1e-10 {
        // Use trigonometric method for real root
        let sqrt_term = 2.0 * (-dep_p / 3.0).sqrt();
        let theta = (3.0 * dep_q / (dep_p * sqrt_term)).acos() / 3.0;
        m = sqrt_term * theta.cos() - resolvent_coeffs[2] / 3.0;
    } else {
        // Use Cardano's formula
        let sqrt_term = (dep_q * dep_q / 4.0 + dep_p * dep_p * dep_p / 27.0)
            .abs()
            .sqrt();
        let sign = if dep_q < 0.0 { 1.0 } else { -1.0 };
        let u = (sign * sqrt_term - dep_q / 2.0).abs().cbrt()
            * (sign * sqrt_term - dep_q / 2.0).signum();
        let v = if u.abs() > 1e-10 {
            -dep_p / (3.0 * u)
        } else {
            0.0
        };
        m = u + v - resolvent_coeffs[2] / 3.0;
    }

    path = path.step(
        Operation::Simplify,
        format!("Resolvent cubic root: m = {}", m),
        Expression::Integer(0),
    );

    // Factor quartic: (y² + m)² = (α + 2m)y² - βy + (m² + αm + γ - γ)
    // Using Ferrari's factorization into two quadratics
    let sqrt_2m_alpha = (2.0 * m + alpha).max(0.0).sqrt();

    // y² + sqrt(2m+α)y + (m + β/(2*sqrt(2m+α))) = 0
    // y² - sqrt(2m+α)y + (m - β/(2*sqrt(2m+α))) = 0
    let term = if sqrt_2m_alpha.abs() > 1e-10 {
        beta / (2.0 * sqrt_2m_alpha)
    } else {
        0.0
    };

    let mut roots = Vec::new();

    // First quadratic: y² + sqrt(2m+α)y + (m + term) = 0
    let a1 = 1.0;
    let b1 = sqrt_2m_alpha;
    let c1 = m + term;
    let disc1 = b1 * b1 - 4.0 * a1 * c1;

    if disc1 >= 0.0 {
        roots.push(simplify_numeric_expression(
            (-b1 + disc1.sqrt()) / 2.0 + shift,
        ));
        roots.push(simplify_numeric_expression(
            (-b1 - disc1.sqrt()) / 2.0 + shift,
        ));
    } else {
        let real = -b1 / 2.0 + shift;
        let imag = (-disc1).sqrt() / 2.0;
        roots.push(Expression::Complex(num_complex::Complex64::new(real, imag)));
        roots.push(Expression::Complex(num_complex::Complex64::new(
            real, -imag,
        )));
    }

    // Second quadratic: y² - sqrt(2m+α)y + (m - term) = 0
    let b2 = -sqrt_2m_alpha;
    let c2 = m - term;
    let disc2 = b2 * b2 - 4.0 * a1 * c2;

    if disc2 >= 0.0 {
        roots.push(simplify_numeric_expression(
            (-b2 + disc2.sqrt()) / 2.0 + shift,
        ));
        roots.push(simplify_numeric_expression(
            (-b2 - disc2.sqrt()) / 2.0 + shift,
        ));
    } else {
        let real = -b2 / 2.0 + shift;
        let imag = (-disc2).sqrt() / 2.0;
        roots.push(Expression::Complex(num_complex::Complex64::new(real, imag)));
        roots.push(Expression::Complex(num_complex::Complex64::new(
            real, -imag,
        )));
    }

    path = path.step(
        Operation::Simplify,
        "Applied Ferrari's method".to_string(),
        roots[0].clone(),
    );

    let resolution_path = path.finish(roots[0].clone());
    Ok((Solution::Multiple(roots), resolution_path))
}

/// Solve polynomial of degree 5+ using numerical methods (Durand-Kerner).
fn solve_polynomial_numerically(
    coeffs: &[f64],
    _var: &str,
    mut path: ResolutionPathBuilder,
) -> SolverResult<(Solution, ResolutionPath)> {
    let degree = coeffs.len() - 1;
    if degree < 1 {
        return Err(SolverError::CannotSolve("Invalid polynomial".to_string()));
    }

    // Find leading coefficient
    let leading = coeffs[degree];
    if leading.abs() < 1e-15 {
        return Err(SolverError::CannotSolve(
            "Leading coefficient is zero".to_string(),
        ));
    }

    path = path.step(
        Operation::Simplify,
        format!(
            "Solving degree {} polynomial numerically (Durand-Kerner method)",
            degree
        ),
        Expression::Integer(0),
    );

    // Initial guess: roots evenly spaced on a circle
    let radius = 1.0
        + coeffs
            .iter()
            .take(degree)
            .map(|c| (c / leading).abs())
            .fold(0.0, f64::max);
    let mut roots: Vec<num_complex::Complex64> = (0..degree)
        .map(|k| {
            let angle = 2.0 * std::f64::consts::PI * (k as f64) / (degree as f64) + 0.4;
            num_complex::Complex64::new(radius * angle.cos(), radius * angle.sin())
        })
        .collect();

    // Durand-Kerner iteration
    let max_iter = 100;
    let tolerance = 1e-12;

    for _ in 0..max_iter {
        let mut max_change: f64 = 0.0;

        for i in 0..degree {
            // Evaluate polynomial at roots[i]
            let mut p_val = num_complex::Complex64::new(0.0, 0.0);
            let mut power = num_complex::Complex64::new(1.0, 0.0);
            for &coeff in coeffs.iter() {
                p_val += num_complex::Complex64::new(coeff, 0.0) * power;
                power *= roots[i];
            }

            // Compute denominator product
            let mut denom = num_complex::Complex64::new(1.0, 0.0);
            for j in 0..degree {
                if i != j {
                    denom *= roots[i] - roots[j];
                }
            }

            if denom.norm() > 1e-15 {
                let delta = p_val / denom;
                roots[i] -= delta;
                max_change = max_change.max(delta.norm());
            }
        }

        if max_change < tolerance {
            break;
        }
    }

    // Convert to Expression
    let root_exprs: Vec<Expression> = roots
        .iter()
        .map(|r| {
            if r.im.abs() < 1e-10 {
                simplify_numeric_expression(r.re)
            } else {
                Expression::Complex(*r)
            }
        })
        .collect();

    path = path.step(
        Operation::Simplify,
        format!("Found {} roots numerically", degree),
        root_exprs[0].clone(),
    );

    let resolution_path = path.finish(root_exprs[0].clone());
    Ok((Solution::Multiple(root_exprs), resolution_path))
}

/// Check if an expression is linear with respect to a specific variable.
/// An expression is linear in variable x if:
/// - x appears to at most power 1
/// - x does not appear in denominators
/// - x does not appear multiplied by itself
/// - x does not appear in functions
fn is_linear_in_variable(expr: &Expression, var: &str) -> bool {
    match expr {
        Expression::Integer(_)
        | Expression::Rational(_)
        | Expression::Float(_)
        | Expression::Complex(_)
        | Expression::Constant(_) => true,

        Expression::Variable(_v) => {
            // The target variable itself is linear
            true
        }

        Expression::Unary(_, inner) => is_linear_in_variable(inner, var),

        Expression::Binary(op, left, right) => {
            let left_has_var = contains_variable(left, var);
            let right_has_var = contains_variable(right, var);

            match op {
                BinaryOp::Add | BinaryOp::Sub => {
                    // x + y and x - y are linear if both sides are linear
                    is_linear_in_variable(left, var) && is_linear_in_variable(right, var)
                }
                BinaryOp::Mul => {
                    // For multiplication to be linear in x, at most one side can contain x
                    if left_has_var && right_has_var {
                        // x * x or x * f(x) is not linear
                        false
                    } else {
                        // a * x is linear
                        is_linear_in_variable(left, var) && is_linear_in_variable(right, var)
                    }
                }
                BinaryOp::Div => {
                    // x / a is linear, but a / x is not
                    if right_has_var {
                        false // Variable in denominator makes it non-linear
                    } else {
                        is_linear_in_variable(left, var)
                    }
                }
                _ => false,
            }
        }

        Expression::Power(base, exp) => {
            // x^2 is not linear, but a^x could be (though we don't handle that in Phase 1)
            // For Phase 1, we only allow constant powers where base doesn't have the variable
            !contains_variable(base, var) && is_linear_in_variable(exp, var)
        }

        Expression::Function(_, _) => {
            // For Phase 1, functions are not supported
            false
        }
    }
}

/// Quadratic equation solver for equations of the form ax² + bx + c = 0.
///
/// Solves second-degree polynomial equations in one variable using the quadratic
/// formula and returns either zero, one, or two real solutions, or two complex solutions.
///
/// # Mathematical Foundation
///
/// A quadratic equation has the general form:
/// ```text
/// ax² + bx + c = 0    where a ≠ 0
/// ```
///
/// The solution is obtained using the quadratic formula:
/// ```text
/// x = (-b ± √(b² - 4ac)) / (2a)
/// ```
///
/// The discriminant Δ = b² - 4ac determines the nature of the roots:
/// - Δ > 0: Two distinct real roots
/// - Δ = 0: One repeated real root (multiplicity 2)
/// - Δ < 0: Two complex conjugate roots
///
/// # TODO: Planned Implementation
///
/// ## Phase 1: Real Roots Only
/// - Extract coefficients a, b, c from equation
/// - Compute discriminant Δ = b² - 4ac
/// - Handle discriminant cases:
///   - Δ > 0: Return [`Solution::Multiple`] with two roots
///   - Δ = 0: Return [`Solution::Unique`] with repeated root
///   - Δ < 0: Return [`SolverError::NoSolution`] (defer complex support)
/// - Apply quadratic formula: x = (-b ± √Δ) / (2a)
/// - Record resolution steps showing discriminant and formula application
///
/// ## Phase 2: Degenerate Cases
/// - Detect when a = 0 (linear equation, not quadratic)
///   - Delegate to [`LinearSolver`]
/// - Detect when a = b = 0 (constant equation)
///   - Return [`Solution::None`] if c ≠ 0
///   - Return [`Solution::Infinite`] if c = 0
///
/// ## Phase 3: Complex Root Support
/// - When Δ < 0, compute complex conjugate pairs
/// - Return [`Solution::Multiple`] with complex expressions
/// - Use [`Expression::Complex`] for representation
/// - Example: x² + 1 = 0 → x = ±i
///
/// ## Phase 4: Alternative Forms
/// - Vertex form: a(x - h)² + k = 0
/// - Factored form: a(x - r₁)(x - r₂) = 0
/// - Recognize perfect square trinomials
/// - Completing the square method (for educational purposes)
///
/// # Limitations (Current)
///
/// - **NOT YET IMPLEMENTED**: Always returns error
/// - Cannot handle equations where variable appears non-quadratically
/// - Cannot solve bivariate quadratics (e.g., x² + xy + y² = 0)
/// - Cannot handle parametric coefficients requiring symbolic computation
///
/// # See Also
///
/// - [`LinearSolver`]: For degenerate case when a = 0
/// - [`PolynomialSolver`]: General polynomial solver (uses QuadraticSolver for degree 2)
/// - [`SmartSolver`]: Automatically selects QuadraticSolver for quadratic equations
#[derive(Debug, Default)]
pub struct QuadraticSolver;

impl QuadraticSolver {
    /// Create a new quadratic equation solver.
    ///
    /// # Examples
    ///
    /// ```
    /// use thales::solver::QuadraticSolver;
    ///
    /// let solver = QuadraticSolver::new();
    /// ```
    pub fn new() -> Self {
        Self
    }
}

impl Solver for QuadraticSolver {
    fn solve(
        &self,
        equation: &Equation,
        variable: &Variable,
    ) -> SolverResult<(Solution, ResolutionPath)> {
        let var_name = &variable.name;

        // Initialize resolution path
        let initial_expr = Expression::Binary(
            BinaryOp::Sub,
            Box::new(equation.left.clone()),
            Box::new(equation.right.clone()),
        );
        let mut path = ResolutionPathBuilder::new(initial_expr.clone());

        // Check if variable appears in equation
        if !contains_variable(&equation.left, var_name)
            && !contains_variable(&equation.right, var_name)
        {
            return Err(SolverError::CannotSolve(format!(
                "Variable '{}' not found in equation",
                var_name
            )));
        }

        // Extract coefficients a, b, c from ax² + bx + c = 0
        // Move everything to left side: left - right = 0
        let combined = Expression::Binary(
            BinaryOp::Sub,
            Box::new(equation.left.clone()),
            Box::new(equation.right.clone()),
        )
        .simplify();

        // Extract polynomial coefficients
        let (a, b, c) = extract_quadratic_coefficients(&combined, var_name);

        // Add step showing coefficients
        path = path.step(
            Operation::Simplify,
            format!("Identified coefficients: a={}, b={}, c={}", a, b, c),
            combined.clone(),
        );

        // Check for degenerate case (a = 0)
        if a.abs() < 1e-15 {
            if b.abs() < 1e-15 {
                if c.abs() < 1e-15 {
                    let resolution_path = path.finish(Expression::Integer(0));
                    return Ok((Solution::Infinite, resolution_path));
                } else {
                    return Err(SolverError::NoSolution);
                }
            }
            // Linear equation: bx + c = 0 -> x = -c/b
            let solution = Expression::Float(-c / b);
            let resolution_path = path.finish(solution.clone());
            return Ok((Solution::Unique(solution), resolution_path));
        }

        // Compute discriminant Δ = b² - 4ac
        let discriminant = b * b - 4.0 * a * c;

        path = path.step(
            Operation::Simplify,
            format!("Computed discriminant: Δ = b² - 4ac = {}", discriminant),
            combined.clone(),
        );

        let epsilon = 1e-15;
        if discriminant > epsilon {
            // Two distinct real roots
            let sqrt_disc = discriminant.sqrt();
            let x1 = (-b + sqrt_disc) / (2.0 * a);
            let x2 = (-b - sqrt_disc) / (2.0 * a);

            let root1 = simplify_numeric_expression(x1);
            let root2 = simplify_numeric_expression(x2);

            path = path.step(
                Operation::Simplify,
                format!("Quadratic formula: x = (-b ± √Δ)/(2a) = {} or {}", x1, x2),
                root1.clone(),
            );

            let resolution_path = path.finish(root1.clone());
            Ok((Solution::Multiple(vec![root1, root2]), resolution_path))
        } else if discriminant.abs() <= epsilon {
            // One repeated real root
            let x = -b / (2.0 * a);
            let root = simplify_numeric_expression(x);

            path = path.step(
                Operation::Simplify,
                format!("Quadratic formula (Δ = 0): x = -b/(2a) = {}", x),
                root.clone(),
            );

            let resolution_path = path.finish(root.clone());
            Ok((Solution::Unique(root), resolution_path))
        } else {
            // Complex roots: x = -b/(2a) ± i√(-Δ)/(2a)
            let real_part = -b / (2.0 * a);
            let imag_part = (-discriminant).sqrt() / (2.0 * a);

            let root1 = Expression::Complex(num_complex::Complex64::new(real_part, imag_part));
            let root2 = Expression::Complex(num_complex::Complex64::new(real_part, -imag_part));

            path = path.step(
                Operation::Simplify,
                format!("Complex roots: x = {} ± {}i", real_part, imag_part),
                root1.clone(),
            );

            let resolution_path = path.finish(root1.clone());
            Ok((Solution::Multiple(vec![root1, root2]), resolution_path))
        }
    }

    fn can_solve(&self, equation: &Equation) -> bool {
        // Check if equation has quadratic terms
        has_obvious_nonlinearity(&equation.left) || has_obvious_nonlinearity(&equation.right)
    }
}

/// Polynomial equation solver for general degree n polynomials.
///
/// Solves polynomial equations in one variable using closed-form algebraic formulas
/// for degrees 1-4, and numerical methods for higher degrees.
///
/// # Mathematical Foundation
///
/// A polynomial equation has the general form:
/// ```text
/// aₙxⁿ + aₙ₋₁xⁿ⁻¹ + ... + a₂x² + a₁x + a₀ = 0
/// ```
/// where n is the degree and aₙ ≠ 0.
///
/// # Solution Methods by Degree
///
/// - **Degree 1 (Linear)**: Direct division: x = -a₀/a₁
/// - **Degree 2 (Quadratic)**: Quadratic formula: x = (-b ± √(b²-4ac))/(2a)
/// - **Degree 3 (Cubic)**: Cardano's formula or trigonometric method
/// - **Degree 4 (Quartic)**: Ferrari's method or resolvent cubic
/// - **Degree 5+ (Quintic and higher)**: Numerical root-finding methods
///
/// # TODO: Planned Implementation
///
/// ## Phase 1: Degree Detection and Delegation
/// - Analyze equation to extract polynomial form
/// - Determine degree of polynomial
/// - Delegate to specialized solver:
///   - Degree 1 → [`LinearSolver`]
///   - Degree 2 → [`QuadraticSolver`]
///   - Degree 3 → Cubic formula implementation
///   - Degree 4 → Quartic formula implementation
///   - Degree 5+ → Numerical solver (see below)
/// - Return appropriate [`Solution`] variant based on number of roots
///
/// ## Phase 2: Cubic Formula (Degree 3)
/// - Normalize to depressed cubic: t³ + pt + q = 0
/// - Compute discriminant: Δ = -4p³ - 27q²
/// - Handle three cases:
///   - Δ > 0: Three distinct real roots (trigonometric method)
///   - Δ = 0: Repeated roots (algebraic method)
///   - Δ < 0: One real root, two complex conjugate roots
/// - Transform back to original variable
/// - Record resolution steps
///
/// ## Phase 3: Quartic Formula (Degree 4)
/// - Normalize to depressed quartic: y⁴ + py² + qy + r = 0
/// - Solve resolvent cubic equation
/// - Use resolvent root to factor quartic
/// - Solve two quadratic equations
/// - Combine roots from both quadratics
/// - Transform back to original variable
///
/// ## Phase 4: Numerical Methods (Degree 5+)
/// - By Abel-Ruffini theorem, no general algebraic solution exists
/// - Implement numerical root-finding:
///   - Newton-Raphson method for initial root approximation
///   - Polynomial deflation to find remaining roots
///   - Durand-Kerner method for simultaneous approximation
///   - Aberth method for robust convergence
/// - Link to `numerical` module for implementation
/// - Return [`Solution::Multiple`] with approximate roots
/// - Note: Numerical solutions are approximate, not symbolic
///
/// ## Phase 5: Special Polynomial Forms
/// - Recognize and optimize special cases:
///   - Binomial: xⁿ - a = 0 (nth roots of a)
///   - Quadratic form: x²ⁿ + bxⁿ + c = 0 (substitute u = xⁿ)
///   - Palindromic: coefficients symmetric
///   - Reciprocal: x = 1/x symmetry
///   - Cyclotomic: roots of unity
/// - Use specialized algorithms for efficiency
///
/// # Limitations (Current)
///
/// - **NOT YET IMPLEMENTED**: Always returns error
/// - Cannot handle multivariate polynomials (e.g., x² + xy + y²)
/// - Cannot handle rational functions (polynomial ÷ polynomial)
/// - Cannot handle polynomials with transcendental coefficients
/// - Cannot handle symbolic/parametric coefficients
///
/// # See Also
///
/// - [`LinearSolver`]: Specialized solver for degree 1
/// - [`QuadraticSolver`]: Specialized solver for degree 2
/// - [`SmartSolver`]: Automatically selects PolynomialSolver for polynomial equations
///
/// # References
///
/// - [Cubic function](https://en.wikipedia.org/wiki/Cubic_function)
/// - [Quartic function](https://en.wikipedia.org/wiki/Quartic_function)
/// - [Abel-Ruffini theorem](https://en.wikipedia.org/wiki/Abel%E2%80%93Ruffini_theorem)
/// - [Newton's method](https://en.wikipedia.org/wiki/Newton%27s_method)
/// - [Durand-Kerner method](https://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method)
#[derive(Debug, Default)]
pub struct PolynomialSolver;

impl PolynomialSolver {
    /// Create a new polynomial equation solver.
    ///
    /// # Examples
    ///
    /// ```
    /// use thales::solver::PolynomialSolver;
    ///
    /// let solver = PolynomialSolver::new();
    /// ```
    pub fn new() -> Self {
        Self
    }
}

impl Solver for PolynomialSolver {
    fn solve(
        &self,
        equation: &Equation,
        variable: &Variable,
    ) -> SolverResult<(Solution, ResolutionPath)> {
        let var_name = &variable.name;

        // Initialize resolution path
        let initial_expr = Expression::Binary(
            BinaryOp::Sub,
            Box::new(equation.left.clone()),
            Box::new(equation.right.clone()),
        );
        let mut path = ResolutionPathBuilder::new(initial_expr.clone());

        // Check if variable appears in equation
        if !contains_variable(&equation.left, var_name)
            && !contains_variable(&equation.right, var_name)
        {
            return Err(SolverError::CannotSolve(format!(
                "Variable '{}' not found in equation",
                var_name
            )));
        }

        // Move everything to left side: left - right = 0
        let combined = Expression::Binary(
            BinaryOp::Sub,
            Box::new(equation.left.clone()),
            Box::new(equation.right.clone()),
        )
        .simplify();

        // Determine degree of polynomial
        let degree = get_polynomial_degree(&combined, var_name);

        path = path.step(
            Operation::Simplify,
            format!("Identified polynomial of degree {}", degree),
            combined.clone(),
        );

        match degree {
            0 => {
                // Constant equation: check if 0 = 0
                if let Some(val) = combined.evaluate(&HashMap::new()) {
                    if val.abs() < 1e-15 {
                        let resolution_path = path.finish(Expression::Integer(0));
                        return Ok((Solution::Infinite, resolution_path));
                    } else {
                        return Err(SolverError::NoSolution);
                    }
                }
                Err(SolverError::CannotSolve(
                    "Cannot evaluate constant expression".to_string(),
                ))
            }
            1 => {
                // Delegate to LinearSolver
                LinearSolver::new().solve(equation, variable)
            }
            2 => {
                // Delegate to QuadraticSolver
                QuadraticSolver::new().solve(equation, variable)
            }
            3 => {
                // Cubic: Cardano's formula
                let coeffs = extract_polynomial_coefficients(&combined, var_name, 3);
                solve_cubic(&coeffs, var_name, path)
            }
            4 => {
                // Quartic: Ferrari's method
                let coeffs = extract_polynomial_coefficients(&combined, var_name, 4);
                solve_quartic(&coeffs, var_name, path)
            }
            _ => {
                // Higher degree: numerical methods
                let coeffs = extract_polynomial_coefficients(&combined, var_name, degree);
                solve_polynomial_numerically(&coeffs, var_name, path)
            }
        }
    }

    fn can_solve(&self, equation: &Equation) -> bool {
        // Check if equation is polynomial (no transcendental functions)
        is_polynomial_expression(&equation.left) && is_polynomial_expression(&equation.right)
    }
}

/// Transcendental equation solver for equations with trig, exp, and log functions.
///
/// Solves equations involving transcendental functions - functions that cannot be
/// expressed in terms of algebraic operations alone. This includes trigonometric,
/// exponential, and logarithmic functions.
///
/// # Supported Equation Types
///
/// The solver recognizes and solves three categories of transcendental equations:
///
/// ## Trigonometric Equations
///
/// Equations with sin, cos, tan and their inverses:
/// - `sin(x) = a` → `x = asin(a)` (requires |a| ≤ 1)
/// - `cos(x) = a` → `x = acos(a)` (requires |a| ≤ 1)
/// - `tan(x) = a` → `x = atan(a)` (no domain restriction)
/// - `c * sin(x) = b` → `x = asin(b/c)`
/// - `sin(c*x) = a` → `x = asin(a)/c`
///
/// ## Logarithmic Equations
///
/// Equations with natural log, log base 10, and arbitrary base logarithms:
/// - `ln(x) = a` → `x = exp(a)`
/// - `log10(x) = a` → `x = 10^a`
/// - `log(x, b) = a` → `x = b^a`
/// - `c * ln(x) = a` → `x = exp(a/c)`
///
/// ## Exponential Equations
///
/// Equations with exponential functions:
/// - `exp(x) = a` → `x = ln(a)`
/// - `a^x = b` → `x = ln(b)/ln(a)` (change of base formula)
/// - `exp(c*x) = a` → `x = ln(a)/c`
/// - `c * exp(x) = a` → `x = ln(a/c)`
///
/// # Domain Validation
///
/// The solver automatically validates domain restrictions for inverse functions:
/// - **asin(x)** and **acos(x)** require `-1 ≤ x ≤ 1`
/// - **ln(x)** and **log(x)** require `x > 0` (validated during evaluation)
/// - **atan(x)** has no domain restrictions
///
/// When a domain restriction is violated, the solver returns a [`SolverError::Other`]
/// with a descriptive error message.
///
/// # Limitations
///
/// The solver currently handles only equations where:
/// 1. The variable appears in a single transcendental function call
/// 2. The function can be inverted by applying its inverse function
/// 3. No products or compositions of transcendental functions with the variable
///
/// For example, **cannot solve**:
/// - `sin(x) * cos(x) = 0.5` (product of functions)
/// - `sin(cos(x)) = 0.5` (composition of functions)
/// - `sin(x) + cos(x) = 1` (sum of different functions)
///
/// # Examples
///
/// ## Solving sin(x) = 0.5
///
/// ```
/// use thales::solver::{TranscendentalSolver, Solver};
/// use thales::ast::{Equation, Expression, Variable, Function};
///
/// // Build equation: sin(x) = 0.5
/// let x = Expression::Variable(Variable::new("x"));
/// let sin_x = Expression::Function(Function::Sin, vec![x]);
/// let half = Expression::Float(0.5);
/// let equation = Equation::new("trig", sin_x, half);
///
/// let solver = TranscendentalSolver::new();
/// let (solution, path) = solver.solve(&equation, &Variable::new("x")).unwrap();
///
/// // Solution is x = asin(0.5) ≈ 0.5236 radians (30 degrees)
/// # use thales::solver::Solution;
/// # use std::collections::HashMap;
/// # match solution {
/// #     Solution::Unique(expr) => {
/// #         let result = expr.evaluate(&HashMap::new()).unwrap();
/// #         assert!((result - 0.5235987755982989).abs() < 1e-10);
/// #     }
/// #     _ => panic!("Expected unique solution"),
/// # }
/// ```
///
/// ## Solving ln(x) = 2
///
/// ```
/// use thales::solver::{TranscendentalSolver, Solver};
/// use thales::ast::{Equation, Expression, Variable, Function};
///
/// // Build equation: ln(x) = 2
/// let x = Expression::Variable(Variable::new("x"));
/// let ln_x = Expression::Function(Function::Ln, vec![x]);
/// let two = Expression::Integer(2);
/// let equation = Equation::new("log", ln_x, two);
///
/// let solver = TranscendentalSolver::new();
/// let (solution, path) = solver.solve(&equation, &Variable::new("x")).unwrap();
///
/// // Solution is x = exp(2) ≈ 7.389
/// # use thales::solver::Solution;
/// # use std::collections::HashMap;
/// # match solution {
/// #     Solution::Unique(expr) => {
/// #         let result = expr.evaluate(&HashMap::new()).unwrap();
/// #         assert!((result - std::f64::consts::E.powi(2)).abs() < 1e-10);
/// #     }
/// #     _ => panic!("Expected unique solution"),
/// # }
/// ```
///
/// ## Solving 2^x = 8
///
/// ```
/// use thales::solver::{TranscendentalSolver, Solver};
/// use thales::ast::{Equation, Expression, Variable};
///
/// // Build equation: 2^x = 8
/// let x = Expression::Variable(Variable::new("x"));
/// let two_pow_x = Expression::Power(
///     Box::new(Expression::Integer(2)),
///     Box::new(x),
/// );
/// let eight = Expression::Integer(8);
/// let equation = Equation::new("exp", two_pow_x, eight);
///
/// let solver = TranscendentalSolver::new();
/// let (solution, path) = solver.solve(&equation, &Variable::new("x")).unwrap();
///
/// // Solution is x = ln(8)/ln(2) = 3
/// # use thales::solver::Solution;
/// # use std::collections::HashMap;
/// # match solution {
/// #     Solution::Unique(expr) => {
/// #         let result = expr.evaluate(&HashMap::new()).unwrap();
/// #         assert!((result - 3.0).abs() < 1e-10);
/// #     }
/// #     _ => panic!("Expected unique solution"),
/// # }
/// ```
///
/// ## Domain Validation
///
/// Domain restrictions are automatically validated. For example, attempting to solve
/// `sin(x) = 2` (which would require `asin(2)`) will fail because |sin(x)| ≤ 1 always.
/// The solver returns `Err(SolverError::CannotSolve(...))` when pattern matching fails
/// due to invalid domains.
///
/// # See Also
///
/// - [`LinearSolver`]: For linear equations (ax + b = c)
/// - [`SmartSolver`]: Automatically selects TranscendentalSolver for transcendental equations
/// - [`solve_for`]: High-level API that handles solver selection and value substitution
#[derive(Debug, Default)]
pub struct TranscendentalSolver;

impl TranscendentalSolver {
    /// Creates a new transcendental equation solver.
    pub fn new() -> Self {
        Self
    }

    /// Try to solve a trigonometric equation for the target variable.
    ///
    /// Attempts to match and solve equations involving sin, cos, or tan by applying
    /// the appropriate inverse function (asin, acos, atan). The method validates
    /// domain restrictions for asin and acos.
    ///
    /// # Supported Patterns
    ///
    /// - `sin(x) = a` or `a = sin(x)` → `x = asin(a)` (requires |a| ≤ 1)
    /// - `cos(x) = a` or `a = cos(x)` → `x = acos(a)` (requires |a| ≤ 1)
    /// - `tan(x) = a` or `a = tan(x)` → `x = atan(a)` (no restriction)
    /// - `c * sin(x) = b` → `x = asin(b/c)`
    /// - `sin(c*x) = a` → `x = asin(a)/c`
    ///
    /// # Domain Validation
    ///
    /// For asin and acos, the input value must satisfy |value| ≤ 1.
    /// If this constraint is violated, the method returns `None` to allow
    /// the caller to propagate the error.
    ///
    /// # Parameters
    ///
    /// - `equation`: The equation to solve
    /// - `variable`: The variable to solve for
    /// - `path`: Resolution path to record the solving steps
    ///
    /// # Returns
    ///
    /// - `Some(expression)` if a valid trigonometric pattern is matched
    /// - `None` if no pattern matches or domain validation fails
    ///
    /// # Examples
    ///
    /// This is a private helper method called by the public [`solve`](Solver::solve) method.
    /// See the [`TranscendentalSolver`] struct documentation for public API examples.
    fn solve_trig_equation(
        &self,
        equation: &Equation,
        variable: &Variable,
        path: &mut ResolutionPath,
    ) -> Option<Expression> {
        let var_name = &variable.name;

        // Pattern: sin(x) = a  →  x = asin(a)
        if let Some((result, func, value)) = self.match_trig_pattern_with_validation(
            &equation.left,
            &equation.right,
            var_name,
            crate::ast::Function::Sin,
            crate::ast::Function::Asin,
        ) {
            // Validate domain before creating result
            if let Err(_e) = Self::validate_trig_domain(value, &func) {
                return None; // Return None to allow error propagation at higher level
            }
            path.add_step(ResolutionStep::new(
                Operation::ApplyFunction("asin".to_string()),
                format!("Apply arcsine to solve sin({}) = value", variable),
                result.clone(),
            ));
            return Some(result);
        }

        // Pattern: a = sin(x)  →  x = asin(a)
        if let Some((result, func, value)) = self.match_trig_pattern_with_validation(
            &equation.right,
            &equation.left,
            var_name,
            crate::ast::Function::Sin,
            crate::ast::Function::Asin,
        ) {
            if let Err(_e) = Self::validate_trig_domain(value, &func) {
                return None;
            }
            path.add_step(ResolutionStep::new(
                Operation::ApplyFunction("asin".to_string()),
                format!("Apply arcsine to solve sin({}) = value", variable),
                result.clone(),
            ));
            return Some(result);
        }

        // Pattern: cos(x) = a  →  x = acos(a)
        if let Some((result, func, value)) = self.match_trig_pattern_with_validation(
            &equation.left,
            &equation.right,
            var_name,
            crate::ast::Function::Cos,
            crate::ast::Function::Acos,
        ) {
            if let Err(_e) = Self::validate_trig_domain(value, &func) {
                return None;
            }
            path.add_step(ResolutionStep::new(
                Operation::ApplyFunction("acos".to_string()),
                format!("Apply arccosine to solve cos({}) = value", variable),
                result.clone(),
            ));
            return Some(result);
        }

        // Pattern: a = cos(x)  →  x = acos(a)
        if let Some((result, func, value)) = self.match_trig_pattern_with_validation(
            &equation.right,
            &equation.left,
            var_name,
            crate::ast::Function::Cos,
            crate::ast::Function::Acos,
        ) {
            if let Err(_e) = Self::validate_trig_domain(value, &func) {
                return None;
            }
            path.add_step(ResolutionStep::new(
                Operation::ApplyFunction("acos".to_string()),
                format!("Apply arccosine to solve cos({}) = value", variable),
                result.clone(),
            ));
            return Some(result);
        }

        // Pattern: tan(x) = a  →  x = atan(a)
        if let Some(result) = self.match_trig_pattern(
            &equation.left,
            &equation.right,
            var_name,
            crate::ast::Function::Tan,
            crate::ast::Function::Atan,
        ) {
            path.add_step(ResolutionStep::new(
                Operation::ApplyFunction("atan".to_string()),
                format!("Apply arctangent to solve tan({}) = value", variable),
                result.clone(),
            ));
            return Some(result);
        }

        // Pattern: a = tan(x)  →  x = atan(a)
        if let Some(result) = self.match_trig_pattern(
            &equation.right,
            &equation.left,
            var_name,
            crate::ast::Function::Tan,
            crate::ast::Function::Atan,
        ) {
            path.add_step(ResolutionStep::new(
                Operation::ApplyFunction("atan".to_string()),
                format!("Apply arctangent to solve tan({}) = value", variable),
                result.clone(),
            ));
            return Some(result);
        }

        None
    }

    /// Match pattern with validation: returns (result, inverse_func, input_value)
    fn match_trig_pattern_with_validation(
        &self,
        left: &Expression,
        right: &Expression,
        var: &str,
        func: crate::ast::Function,
        inverse_func: crate::ast::Function,
    ) -> Option<(Expression, crate::ast::Function, f64)> {
        // Check if right side contains the variable
        if contains_variable(right, var) {
            return None;
        }

        // Try to evaluate the right side as a constant
        let value = match right.evaluate(&HashMap::new()) {
            Some(v) => v,
            None => return None, // Can't validate if not a constant
        };

        // Pattern 1: func(x) = a  →  x = inverse_func(a)
        if let Expression::Function(f, args) = left {
            if *f == func && args.len() == 1 {
                // Check if arg is exactly the variable
                if let Expression::Variable(v) = &args[0] {
                    if v.name == var {
                        let result =
                            Expression::Function(inverse_func.clone(), vec![right.clone()]);
                        return Some((result.simplify(), inverse_func, value));
                    }
                }

                // Check if arg is a linear expression like a*x
                if let Some(coeff) = extract_coefficient(&args[0], var) {
                    // func(a*x) = b  →  a*x = inverse_func(b)  →  x = inverse_func(b) / a
                    let inverse_applied =
                        Expression::Function(inverse_func.clone(), vec![right.clone()]);
                    let result = Expression::Binary(
                        BinaryOp::Div,
                        Box::new(inverse_applied),
                        Box::new(coeff),
                    );
                    return Some((result.simplify(), inverse_func, value));
                }
            }
        }

        // Pattern 2: a * func(x) = b  →  func(x) = b/a  →  x = inverse_func(b/a)
        if let Expression::Binary(BinaryOp::Mul, mul_left, mul_right) = left {
            // Check left side of multiplication
            if let Expression::Function(f, args) = mul_left.as_ref() {
                if *f == func && args.len() == 1 && !contains_variable(mul_right, var) {
                    if let Expression::Variable(v) = &args[0] {
                        if v.name == var {
                            // a * func(x) = b  →  func(x) = b/a  →  x = inverse_func(b/a)
                            let divided = Expression::Binary(
                                BinaryOp::Div,
                                Box::new(right.clone()),
                                Box::new(mul_right.as_ref().clone()),
                            );
                            // Need to evaluate the divided value
                            let divided_value = divided.evaluate(&HashMap::new()).unwrap_or(value);
                            let result = Expression::Function(inverse_func.clone(), vec![divided]);
                            return Some((result.simplify(), inverse_func, divided_value));
                        }
                    }
                }
            }

            // Check right side of multiplication
            if let Expression::Function(f, args) = mul_right.as_ref() {
                if *f == func && args.len() == 1 && !contains_variable(mul_left, var) {
                    if let Expression::Variable(v) = &args[0] {
                        if v.name == var {
                            // func(x) * a = b  →  func(x) = b/a  →  x = inverse_func(b/a)
                            let divided = Expression::Binary(
                                BinaryOp::Div,
                                Box::new(right.clone()),
                                Box::new(mul_left.as_ref().clone()),
                            );
                            let divided_value = divided.evaluate(&HashMap::new()).unwrap_or(value);
                            let result = Expression::Function(inverse_func.clone(), vec![divided]);
                            return Some((result.simplify(), inverse_func, divided_value));
                        }
                    }
                }
            }
        }

        None
    }

    /// Match pattern: func(var) = value or coeff * func(var) = value
    fn match_trig_pattern(
        &self,
        left: &Expression,
        right: &Expression,
        var: &str,
        func: crate::ast::Function,
        inverse_func: crate::ast::Function,
    ) -> Option<Expression> {
        // Check if right side contains the variable
        if contains_variable(right, var) {
            return None;
        }

        // Pattern 1: func(x) = a  →  x = inverse_func(a)
        if let Expression::Function(f, args) = left {
            if *f == func && args.len() == 1 {
                // Check if arg is exactly the variable
                if let Expression::Variable(v) = &args[0] {
                    if v.name == var {
                        let result = Expression::Function(inverse_func, vec![right.clone()]);
                        return Some(result.simplify());
                    }
                }

                // Check if arg is a linear expression like a*x
                if let Some(coeff) = extract_coefficient(&args[0], var) {
                    // func(a*x) = b  →  a*x = inverse_func(b)  →  x = inverse_func(b) / a
                    let inverse_applied = Expression::Function(inverse_func, vec![right.clone()]);
                    let result = Expression::Binary(
                        BinaryOp::Div,
                        Box::new(inverse_applied),
                        Box::new(coeff),
                    );
                    return Some(result.simplify());
                }
            }
        }

        // Pattern 2: a * func(x) = b  →  func(x) = b/a  →  x = inverse_func(b/a)
        if let Expression::Binary(BinaryOp::Mul, mul_left, mul_right) = left {
            // Check left side of multiplication
            if let Expression::Function(f, args) = mul_left.as_ref() {
                if *f == func && args.len() == 1 && !contains_variable(mul_right, var) {
                    if let Expression::Variable(v) = &args[0] {
                        if v.name == var {
                            // a * func(x) = b  →  func(x) = b/a  →  x = inverse_func(b/a)
                            let divided = Expression::Binary(
                                BinaryOp::Div,
                                Box::new(right.clone()),
                                Box::new(mul_right.as_ref().clone()),
                            );
                            let result = Expression::Function(inverse_func, vec![divided]);
                            return Some(result.simplify());
                        }
                    }
                }
            }

            // Check right side of multiplication
            if let Expression::Function(f, args) = mul_right.as_ref() {
                if *f == func && args.len() == 1 && !contains_variable(mul_left, var) {
                    if let Expression::Variable(v) = &args[0] {
                        if v.name == var {
                            // func(x) * a = b  →  func(x) = b/a  →  x = inverse_func(b/a)
                            let divided = Expression::Binary(
                                BinaryOp::Div,
                                Box::new(right.clone()),
                                Box::new(mul_left.as_ref().clone()),
                            );
                            let result = Expression::Function(inverse_func, vec![divided]);
                            return Some(result.simplify());
                        }
                    }
                }
            }
        }

        None
    }

    /// Try to solve a logarithmic equation for the target variable.
    ///
    /// Attempts to match and solve equations involving ln, log10, or log with
    /// arbitrary base by converting to exponential form.
    ///
    /// # Supported Patterns
    ///
    /// - `ln(x) = a` or `a = ln(x)` → `x = exp(a)`
    /// - `log10(x) = a` or `a = log10(x)` → `x = 10^a`
    /// - `log(x, b) = a` or `a = log(x, b)` → `x = b^a`
    /// - `c * ln(x) = a` → `x = exp(a/c)`
    ///
    /// # Mathematical Principle
    ///
    /// The solver uses the inverse relationship between logarithms and exponents:
    /// - If `log_b(x) = a`, then `x = b^a`
    /// - Natural log: `ln(x) = a` → `x = e^a` (using exp function)
    /// - Common log: `log10(x) = a` → `x = 10^a`
    ///
    /// # Parameters
    ///
    /// - `equation`: The equation to solve
    /// - `variable`: The variable to solve for
    /// - `path`: Resolution path to record the solving steps
    ///
    /// # Returns
    ///
    /// - `Some(expression)` if a valid logarithmic pattern is matched
    /// - `None` if no pattern matches
    ///
    /// # Examples
    ///
    /// This is a private helper method called by the public [`solve`](Solver::solve) method.
    /// See the [`TranscendentalSolver`] struct documentation for public API examples.
    fn solve_log_equation(
        &self,
        equation: &Equation,
        variable: &Variable,
        path: &mut ResolutionPath,
    ) -> Option<Expression> {
        let var_name = &variable.name;

        // Pattern: ln(x) = a  →  x = exp(a)
        if let Some(result) = self.match_log_pattern(&equation.left, &equation.right, var_name) {
            path.add_step(ResolutionStep::new(
                Operation::ApplyFunction("exp".to_string()),
                format!("Apply exponential to solve ln({}) = value", variable),
                result.clone(),
            ));
            return Some(result);
        }

        // Pattern: a = ln(x)  →  x = exp(a)
        if let Some(result) = self.match_log_pattern(&equation.right, &equation.left, var_name) {
            path.add_step(ResolutionStep::new(
                Operation::ApplyFunction("exp".to_string()),
                format!("Apply exponential to solve ln({}) = value", variable),
                result.clone(),
            ));
            return Some(result);
        }

        // Pattern: log10(x) = a  →  x = 10^a
        if let Some(result) = self.match_log10_pattern(&equation.left, &equation.right, var_name) {
            path.add_step(ResolutionStep::new(
                Operation::PowerBothSides(Expression::Integer(10)),
                format!("Apply 10^x to solve log10({}) = value", variable),
                result.clone(),
            ));
            return Some(result);
        }

        // Pattern: a = log10(x)  →  x = 10^a
        if let Some(result) = self.match_log10_pattern(&equation.right, &equation.left, var_name) {
            path.add_step(ResolutionStep::new(
                Operation::PowerBothSides(Expression::Integer(10)),
                format!("Apply 10^x to solve log10({}) = value", variable),
                result.clone(),
            ));
            return Some(result);
        }

        // Pattern: log(x, b) = a  →  x = b^a
        if let Some(result) = self.match_log_base_pattern(&equation.left, &equation.right, var_name)
        {
            path.add_step(ResolutionStep::new(
                Operation::ApplyLogProperty("exponential form".to_string()),
                format!(
                    "Convert logarithm to exponential form to solve for {}",
                    variable
                ),
                result.clone(),
            ));
            return Some(result);
        }

        // Pattern: a = log(x, b)  →  x = b^a
        if let Some(result) = self.match_log_base_pattern(&equation.right, &equation.left, var_name)
        {
            path.add_step(ResolutionStep::new(
                Operation::ApplyLogProperty("exponential form".to_string()),
                format!(
                    "Convert logarithm to exponential form to solve for {}",
                    variable
                ),
                result.clone(),
            ));
            return Some(result);
        }

        None
    }

    /// Match pattern: ln(var) = value or coeff * ln(var) = value
    fn match_log_pattern(
        &self,
        left: &Expression,
        right: &Expression,
        var: &str,
    ) -> Option<Expression> {
        // Check if right side contains the variable
        if contains_variable(right, var) {
            return None;
        }

        // Pattern 1: ln(x) = a  →  x = exp(a)
        if let Expression::Function(crate::ast::Function::Ln, args) = left {
            if args.len() == 1 {
                if let Expression::Variable(v) = &args[0] {
                    if v.name == var {
                        let result =
                            Expression::Function(crate::ast::Function::Exp, vec![right.clone()]);
                        return Some(result.simplify());
                    }
                }
            }
        }

        // Pattern 2: a * ln(x) = b  →  ln(x) = b/a  →  x = exp(b/a)
        if let Expression::Binary(BinaryOp::Mul, mul_left, mul_right) = left {
            if let Expression::Function(crate::ast::Function::Ln, args) = mul_left.as_ref() {
                if args.len() == 1 && !contains_variable(mul_right, var) {
                    if let Expression::Variable(v) = &args[0] {
                        if v.name == var {
                            let divided = Expression::Binary(
                                BinaryOp::Div,
                                Box::new(right.clone()),
                                Box::new(mul_right.as_ref().clone()),
                            );
                            let result =
                                Expression::Function(crate::ast::Function::Exp, vec![divided]);
                            return Some(result.simplify());
                        }
                    }
                }
            }

            if let Expression::Function(crate::ast::Function::Ln, args) = mul_right.as_ref() {
                if args.len() == 1 && !contains_variable(mul_left, var) {
                    if let Expression::Variable(v) = &args[0] {
                        if v.name == var {
                            let divided = Expression::Binary(
                                BinaryOp::Div,
                                Box::new(right.clone()),
                                Box::new(mul_left.as_ref().clone()),
                            );
                            let result =
                                Expression::Function(crate::ast::Function::Exp, vec![divided]);
                            return Some(result.simplify());
                        }
                    }
                }
            }
        }

        None
    }

    /// Match pattern: log10(var) = value
    fn match_log10_pattern(
        &self,
        left: &Expression,
        right: &Expression,
        var: &str,
    ) -> Option<Expression> {
        if contains_variable(right, var) {
            return None;
        }

        // Pattern: log10(x) = a  →  x = 10^a
        if let Expression::Function(crate::ast::Function::Log10, args) = left {
            if args.len() == 1 {
                if let Expression::Variable(v) = &args[0] {
                    if v.name == var {
                        let result = Expression::Power(
                            Box::new(Expression::Integer(10)),
                            Box::new(right.clone()),
                        );
                        return Some(result.simplify());
                    }
                }
            }
        }

        None
    }

    /// Match pattern: log(var, base) = value
    fn match_log_base_pattern(
        &self,
        left: &Expression,
        right: &Expression,
        var: &str,
    ) -> Option<Expression> {
        if contains_variable(right, var) {
            return None;
        }

        // Pattern: log(x, b) = a  →  x = b^a
        if let Expression::Function(crate::ast::Function::Log, args) = left {
            if args.len() == 2 {
                if let Expression::Variable(v) = &args[0] {
                    if v.name == var && !contains_variable(&args[1], var) {
                        let result =
                            Expression::Power(Box::new(args[1].clone()), Box::new(right.clone()));
                        return Some(result.simplify());
                    }
                }
            }
        }

        None
    }

    /// Try to solve an exponential equation for the target variable.
    ///
    /// Attempts to match and solve equations with the variable in an exponent
    /// by applying logarithms to isolate the variable.
    ///
    /// # Supported Patterns
    ///
    /// - `exp(x) = a` or `a = exp(x)` → `x = ln(a)`
    /// - `a^x = b` or `b = a^x` → `x = ln(b)/ln(a)` (change of base formula)
    /// - `exp(c*x) = a` → `x = ln(a)/c`
    /// - `c * exp(x) = a` → `x = ln(a/c)`
    /// - `a^(c*x) = b` → `x = ln(b)/(c*ln(a))`
    ///
    /// # Mathematical Principle
    ///
    /// The solver uses the inverse relationship between exponents and logarithms:
    /// - If `e^x = a`, then `x = ln(a)`
    /// - If `b^x = a`, then `x = log_b(a) = ln(a)/ln(b)` (change of base formula)
    ///
    /// # Parameters
    ///
    /// - `equation`: The equation to solve
    /// - `variable`: The variable to solve for
    /// - `path`: Resolution path to record the solving steps
    ///
    /// # Returns
    ///
    /// - `Some(expression)` if a valid exponential pattern is matched
    /// - `None` if no pattern matches
    ///
    /// # Examples
    ///
    /// This is a private helper method called by the public [`solve`](Solver::solve) method.
    /// See the [`TranscendentalSolver`] struct documentation for public API examples.
    fn solve_exp_equation(
        &self,
        equation: &Equation,
        variable: &Variable,
        path: &mut ResolutionPath,
    ) -> Option<Expression> {
        let var_name = &variable.name;

        // Pattern: exp(x) = a  →  x = ln(a)
        if let Some(result) = self.match_exp_pattern(&equation.left, &equation.right, var_name) {
            path.add_step(ResolutionStep::new(
                Operation::ApplyFunction("ln".to_string()),
                format!("Apply natural logarithm to solve exp({}) = value", variable),
                result.clone(),
            ));
            return Some(result);
        }

        // Pattern: a = exp(x)  →  x = ln(a)
        if let Some(result) = self.match_exp_pattern(&equation.right, &equation.left, var_name) {
            path.add_step(ResolutionStep::new(
                Operation::ApplyFunction("ln".to_string()),
                format!("Apply natural logarithm to solve exp({}) = value", variable),
                result.clone(),
            ));
            return Some(result);
        }

        // Pattern: a^x = b  →  x = ln(b) / ln(a)
        if let Some(result) = self.match_power_pattern(&equation.left, &equation.right, var_name) {
            path.add_step(ResolutionStep::new(
                Operation::ApplyLogProperty("change of base".to_string()),
                format!("Apply logarithm to solve for {} in exponent", variable),
                result.clone(),
            ));
            return Some(result);
        }

        // Pattern: b = a^x  →  x = ln(b) / ln(a)
        if let Some(result) = self.match_power_pattern(&equation.right, &equation.left, var_name) {
            path.add_step(ResolutionStep::new(
                Operation::ApplyLogProperty("change of base".to_string()),
                format!("Apply logarithm to solve for {} in exponent", variable),
                result.clone(),
            ));
            return Some(result);
        }

        None
    }

    /// Match pattern: exp(var) = value or coeff * exp(var) = value
    fn match_exp_pattern(
        &self,
        left: &Expression,
        right: &Expression,
        var: &str,
    ) -> Option<Expression> {
        if contains_variable(right, var) {
            return None;
        }

        // Pattern 1: exp(x) = a  →  x = ln(a)
        if let Expression::Function(crate::ast::Function::Exp, args) = left {
            if args.len() == 1 {
                if let Expression::Variable(v) = &args[0] {
                    if v.name == var {
                        let result =
                            Expression::Function(crate::ast::Function::Ln, vec![right.clone()]);
                        return Some(result.simplify());
                    }
                }

                // Pattern: exp(a*x) = b  →  a*x = ln(b)  →  x = ln(b)/a
                if let Some(coeff) = extract_coefficient(&args[0], var) {
                    let ln_applied =
                        Expression::Function(crate::ast::Function::Ln, vec![right.clone()]);
                    let result =
                        Expression::Binary(BinaryOp::Div, Box::new(ln_applied), Box::new(coeff));
                    return Some(result.simplify());
                }
            }
        }

        // Pattern 2: a * exp(x) = b  →  exp(x) = b/a  →  x = ln(b/a)
        if let Expression::Binary(BinaryOp::Mul, mul_left, mul_right) = left {
            if let Expression::Function(crate::ast::Function::Exp, args) = mul_left.as_ref() {
                if args.len() == 1 && !contains_variable(mul_right, var) {
                    if let Expression::Variable(v) = &args[0] {
                        if v.name == var {
                            let divided = Expression::Binary(
                                BinaryOp::Div,
                                Box::new(right.clone()),
                                Box::new(mul_right.as_ref().clone()),
                            );
                            let result =
                                Expression::Function(crate::ast::Function::Ln, vec![divided]);
                            return Some(result.simplify());
                        }
                    }
                }
            }

            if let Expression::Function(crate::ast::Function::Exp, args) = mul_right.as_ref() {
                if args.len() == 1 && !contains_variable(mul_left, var) {
                    if let Expression::Variable(v) = &args[0] {
                        if v.name == var {
                            let divided = Expression::Binary(
                                BinaryOp::Div,
                                Box::new(right.clone()),
                                Box::new(mul_left.as_ref().clone()),
                            );
                            let result =
                                Expression::Function(crate::ast::Function::Ln, vec![divided]);
                            return Some(result.simplify());
                        }
                    }
                }
            }
        }

        None
    }

    /// Match pattern: base^var = value
    fn match_power_pattern(
        &self,
        left: &Expression,
        right: &Expression,
        var: &str,
    ) -> Option<Expression> {
        if contains_variable(right, var) {
            return None;
        }

        // Pattern: a^x = b  →  x = ln(b) / ln(a)
        if let Expression::Power(base, exp) = left {
            if !contains_variable(base, var) && contains_variable(exp, var) {
                // Simple case: a^x = b
                if let Expression::Variable(v) = exp.as_ref() {
                    if v.name == var {
                        let ln_right =
                            Expression::Function(crate::ast::Function::Ln, vec![right.clone()]);
                        let ln_base = Expression::Function(
                            crate::ast::Function::Ln,
                            vec![base.as_ref().clone()],
                        );
                        let result = Expression::Binary(
                            BinaryOp::Div,
                            Box::new(ln_right),
                            Box::new(ln_base),
                        );
                        return Some(result.simplify());
                    }
                }

                // Pattern: a^(b*x) = c  →  b*x = ln(c)/ln(a)  →  x = ln(c)/(b*ln(a))
                if let Some(coeff) = extract_coefficient(exp, var) {
                    let ln_right =
                        Expression::Function(crate::ast::Function::Ln, vec![right.clone()]);
                    let ln_base =
                        Expression::Function(crate::ast::Function::Ln, vec![base.as_ref().clone()]);
                    let divided =
                        Expression::Binary(BinaryOp::Div, Box::new(ln_right), Box::new(ln_base));
                    let result =
                        Expression::Binary(BinaryOp::Div, Box::new(divided), Box::new(coeff));
                    return Some(result.simplify());
                }
            }
        }

        None
    }

    /// Check if an expression contains transcendental functions.
    ///
    /// A transcendental function is one that cannot be expressed as a finite
    /// combination of algebraic operations (addition, subtraction, multiplication,
    /// division, and root extraction). This includes trigonometric, exponential,
    /// and logarithmic functions.
    ///
    /// This helper is used by [`can_solve`](TranscendentalSolver::can_solve) to determine
    /// if an equation is suitable for the TranscendentalSolver.
    ///
    /// # Recognized Transcendental Functions
    ///
    /// - **Trigonometric**: sin, cos, tan, asin, acos, atan
    /// - **Hyperbolic**: sinh, cosh, tanh
    /// - **Exponential**: exp, a^x (when x contains a variable)
    /// - **Logarithmic**: ln, log, log2, log10
    ///
    /// # Examples
    ///
    /// ```
    /// use thales::solver::TranscendentalSolver;
    /// use thales::ast::{Expression, Variable, Function};
    ///
    /// // sin(x) contains transcendental function
    /// let x = Expression::Variable(Variable::new("x"));
    /// let sin_x = Expression::Function(Function::Sin, vec![x.clone()]);
    /// // Note: has_transcendental_function is private, tested via can_solve
    ///
    /// // x^2 does not contain transcendental function (algebraic)
    /// let x_squared = Expression::Power(
    ///     Box::new(x.clone()),
    ///     Box::new(Expression::Integer(2)),
    /// );
    ///
    /// // 2^x contains transcendental function (variable in exponent)
    /// let two_pow_x = Expression::Power(
    ///     Box::new(Expression::Integer(2)),
    ///     Box::new(x.clone()),
    /// );
    /// ```
    fn has_transcendental_function(expr: &Expression) -> bool {
        match expr {
            Expression::Function(func, _) => {
                matches!(
                    func,
                    crate::ast::Function::Sin
                        | crate::ast::Function::Cos
                        | crate::ast::Function::Tan
                        | crate::ast::Function::Asin
                        | crate::ast::Function::Acos
                        | crate::ast::Function::Atan
                        | crate::ast::Function::Sinh
                        | crate::ast::Function::Cosh
                        | crate::ast::Function::Tanh
                        | crate::ast::Function::Exp
                        | crate::ast::Function::Ln
                        | crate::ast::Function::Log
                        | crate::ast::Function::Log2
                        | crate::ast::Function::Log10
                )
            }
            Expression::Unary(_, inner) => Self::has_transcendental_function(inner),
            Expression::Binary(_, left, right) => {
                Self::has_transcendental_function(left) || Self::has_transcendental_function(right)
            }
            Expression::Power(base, exp) => {
                // Check if variable appears in exponent (exponential form)
                has_any_variable(exp)
                    || Self::has_transcendental_function(base)
                    || Self::has_transcendental_function(exp)
            }
            _ => false,
        }
    }

    /// Validate domain restrictions for inverse trigonometric functions.
    ///
    /// Ensures that input values to inverse trigonometric functions satisfy
    /// their domain restrictions. This prevents mathematical errors like
    /// attempting to compute asin(2).
    ///
    /// # Domain Restrictions
    ///
    /// - **asin(x)**: Requires `-1 ≤ x ≤ 1` (domain of arcsine)
    /// - **acos(x)**: Requires `-1 ≤ x ≤ 1` (domain of arccosine)
    /// - **atan(x)**: No restriction (all real numbers)
    ///
    /// # Parameters
    ///
    /// - `value`: The numeric value to validate
    /// - `func`: The inverse trigonometric function being applied
    ///
    /// # Returns
    ///
    /// - `Ok(())` if the value is within the valid domain
    /// - `Err(SolverError::Other)` if the value violates domain restrictions
    ///
    /// # Examples
    ///
    /// This is a private helper method that validates domains automatically during solving.
    /// See the domain error example in the [`TranscendentalSolver`] struct documentation
    /// for how domain validation errors are surfaced through the public API.
    fn validate_trig_domain(value: f64, func: &crate::ast::Function) -> Result<(), SolverError> {
        match func {
            crate::ast::Function::Asin | crate::ast::Function::Acos => {
                if value.abs() > 1.0 {
                    return Err(SolverError::Other(format!(
                        "Domain error: {:?} requires |value| ≤ 1, got {}",
                        func, value
                    )));
                }
            }
            _ => {}
        }
        Ok(())
    }
}

impl Solver for TranscendentalSolver {
    fn solve(
        &self,
        equation: &Equation,
        variable: &Variable,
    ) -> SolverResult<(Solution, ResolutionPath)> {
        let var_name = &variable.name;

        // Check if variable appears in equation
        let left_has_var = contains_variable(&equation.left, var_name);
        let right_has_var = contains_variable(&equation.right, var_name);

        if !left_has_var && !right_has_var {
            return Err(SolverError::CannotSolve(format!(
                "Variable '{}' not found in equation",
                var_name
            )));
        }

        // Initialize resolution path
        let initial_expr = Expression::Binary(
            BinaryOp::Sub,
            Box::new(equation.left.clone()),
            Box::new(equation.right.clone()),
        );
        let mut path = ResolutionPath::new(initial_expr);

        // Try trigonometric equation patterns
        if let Some(result) = self.solve_trig_equation(equation, variable, &mut path) {
            // Validate domain if result is a constant
            if let Expression::Function(func, args) = &result {
                if args.len() == 1 {
                    if let Some(val) = args[0].evaluate(&HashMap::new()) {
                        Self::validate_trig_domain(val, func)?;
                    }
                }
            }

            let evaluated = evaluate_constants(&result);
            path.set_result(evaluated.clone());
            return Ok((Solution::Unique(evaluated), path));
        }

        // Try logarithmic equation patterns
        if let Some(result) = self.solve_log_equation(equation, variable, &mut path) {
            let evaluated = evaluate_constants(&result);
            path.set_result(evaluated.clone());
            return Ok((Solution::Unique(evaluated), path));
        }

        // Try exponential equation patterns
        if let Some(result) = self.solve_exp_equation(equation, variable, &mut path) {
            let evaluated = evaluate_constants(&result);
            path.set_result(evaluated.clone());
            return Ok((Solution::Unique(evaluated), path));
        }

        // If no pattern matched, cannot solve
        Err(SolverError::CannotSolve(
            "Transcendental equation pattern not recognized or too complex".to_string(),
        ))
    }

    fn can_solve(&self, equation: &Equation) -> bool {
        // Check if equation contains transcendental functions
        Self::has_transcendental_function(&equation.left)
            || Self::has_transcendental_function(&equation.right)
    }
}

/// Result type for system solutions.
#[derive(Debug, Clone)]
pub enum SystemSolution {
    /// Unique solution: each variable has exactly one value.
    Unique(HashMap<Variable, Expression>),
    /// Infinite solutions: variables are expressed in terms of free parameters.
    Infinite {
        /// Variables that have specific values.
        bound: HashMap<Variable, Expression>,
        /// Variables that are free parameters (can take any value).
        free: Vec<Variable>,
    },
    /// No solution: the system is inconsistent.
    NoSolution,
}

/// A linear system of equations in matrix form Ax = b.
#[derive(Debug, Clone)]
pub struct LinearSystem {
    /// Coefficient matrix A (rows × cols where cols = number of variables).
    coefficients: Vec<Vec<f64>>,
    /// Constant vector b.
    constants: Vec<f64>,
    /// Variable names corresponding to columns.
    variables: Vec<Variable>,
}

impl LinearSystem {
    /// Create a linear system from equations and variables.
    ///
    /// Extracts coefficients from linear equations of the form:
    /// a₁x₁ + a₂x₂ + ... + aₙxₙ = b
    ///
    /// # Errors
    ///
    /// Returns an error if any equation is not linear in the given variables.
    pub fn from_equations(equations: &[Equation], variables: &[Variable]) -> SolverResult<Self> {
        let n_eqs = equations.len();
        let n_vars = variables.len();

        if n_eqs == 0 || n_vars == 0 {
            return Err(SolverError::Other("Empty system".to_string()));
        }

        let mut coefficients = Vec::with_capacity(n_eqs);
        let mut constants = Vec::with_capacity(n_eqs);

        for eq in equations {
            // Move everything to the left: left - right = 0
            let combined = Expression::Binary(
                BinaryOp::Sub,
                Box::new(eq.left.clone()),
                Box::new(eq.right.clone()),
            )
            .simplify();

            let (row, constant) = Self::extract_linear_coefficients(&combined, variables)?;
            coefficients.push(row);
            constants.push(-constant); // Move constant to RHS
        }

        Ok(Self {
            coefficients,
            constants,
            variables: variables.to_vec(),
        })
    }

    /// Extract linear coefficients from an expression.
    /// Returns (coefficients for each variable, constant term).
    fn extract_linear_coefficients(
        expr: &Expression,
        variables: &[Variable],
    ) -> SolverResult<(Vec<f64>, f64)> {
        let mut coeffs = vec![0.0; variables.len()];
        let mut constant = 0.0;

        // Collect all additive terms
        let terms = Self::collect_additive_terms(expr);

        for term in terms {
            // Check which variable this term contains
            let mut found_var = false;
            for (i, var) in variables.iter().enumerate() {
                if term.contains_variable(&var.name) {
                    // Extract coefficient
                    let coeff = Self::extract_coefficient(&term, var)?;
                    coeffs[i] += coeff;
                    found_var = true;
                    break;
                }
            }

            if !found_var {
                // This is a constant term
                let empty_vars: HashMap<String, f64> = HashMap::new();
                match term.evaluate(&empty_vars) {
                    Some(val) => constant += val,
                    None => {
                        return Err(SolverError::Other(format!(
                            "Cannot evaluate constant term: {}",
                            term
                        )));
                    }
                }
            }
        }

        Ok((coeffs, constant))
    }

    /// Collect all terms that are added or subtracted.
    fn collect_additive_terms(expr: &Expression) -> Vec<Expression> {
        match expr {
            Expression::Binary(BinaryOp::Add, left, right) => {
                let mut terms = Self::collect_additive_terms(left);
                terms.extend(Self::collect_additive_terms(right));
                terms
            }
            Expression::Binary(BinaryOp::Sub, left, right) => {
                let mut terms = Self::collect_additive_terms(left);
                // Negate right side terms
                for term in Self::collect_additive_terms(right) {
                    terms.push(Expression::Unary(UnaryOp::Neg, Box::new(term)));
                }
                terms
            }
            _ => vec![expr.clone()],
        }
    }

    /// Extract the coefficient of a variable from a term.
    fn extract_coefficient(term: &Expression, var: &Variable) -> SolverResult<f64> {
        match term {
            // Just the variable: coefficient is 1
            Expression::Variable(v) if v.name == var.name => Ok(1.0),

            // Negated variable: coefficient is -1
            Expression::Unary(UnaryOp::Neg, inner) => {
                let inner_coeff = Self::extract_coefficient(inner, var)?;
                Ok(-inner_coeff)
            }

            // Multiplication: c * x or x * c
            Expression::Binary(BinaryOp::Mul, left, right) => {
                let left_has_var = left.contains_variable(&var.name);
                let right_has_var = right.contains_variable(&var.name);

                if left_has_var && right_has_var {
                    return Err(SolverError::Other(format!(
                        "Non-linear term: {} * {} both contain {}",
                        left, right, var.name
                    )));
                }

                if left_has_var {
                    // Right should be the coefficient
                    let empty: HashMap<String, f64> = HashMap::new();
                    let coeff = right.evaluate(&empty).ok_or_else(|| {
                        SolverError::Other(format!("Cannot evaluate coefficient: {}", right))
                    })?;
                    let var_coeff = Self::extract_coefficient(left, var)?;
                    Ok(coeff * var_coeff)
                } else {
                    // Left should be the coefficient
                    let empty: HashMap<String, f64> = HashMap::new();
                    let coeff = left.evaluate(&empty).ok_or_else(|| {
                        SolverError::Other(format!("Cannot evaluate coefficient: {}", left))
                    })?;
                    let var_coeff = Self::extract_coefficient(right, var)?;
                    Ok(coeff * var_coeff)
                }
            }

            // Division: x / c
            Expression::Binary(BinaryOp::Div, left, right) => {
                if right.contains_variable(&var.name) {
                    return Err(SolverError::Other(format!(
                        "Non-linear: variable {} in denominator",
                        var.name
                    )));
                }
                let empty: HashMap<String, f64> = HashMap::new();
                let divisor = right.evaluate(&empty).ok_or_else(|| {
                    SolverError::Other(format!("Cannot evaluate divisor: {}", right))
                })?;
                if divisor.abs() < 1e-15 {
                    return Err(SolverError::DivisionByZero);
                }
                let var_coeff = Self::extract_coefficient(left, var)?;
                Ok(var_coeff / divisor)
            }

            _ => {
                // Check if this term contains the variable at all
                if term.contains_variable(&var.name) {
                    Err(SolverError::Other(format!(
                        "Cannot extract coefficient from: {}",
                        term
                    )))
                } else {
                    Ok(0.0)
                }
            }
        }
    }

    /// Solve the linear system using Gaussian elimination with partial pivoting.
    pub fn solve(&self) -> SolverResult<SystemSolution> {
        let n_eqs = self.coefficients.len();
        let n_vars = self.variables.len();

        // Create augmented matrix [A|b]
        let mut augmented: Vec<Vec<f64>> = self
            .coefficients
            .iter()
            .zip(self.constants.iter())
            .map(|(row, &c)| {
                let mut new_row = row.clone();
                new_row.push(c);
                new_row
            })
            .collect();

        // Forward elimination with partial pivoting
        let mut pivot_row = 0;
        let mut pivot_cols = Vec::new(); // Track which columns had pivots

        for col in 0..n_vars {
            if pivot_row >= n_eqs {
                break;
            }

            // Find pivot (largest absolute value in column)
            let mut max_row = pivot_row;
            let mut max_val = augmented[pivot_row][col].abs();
            for row in (pivot_row + 1)..n_eqs {
                if augmented[row][col].abs() > max_val {
                    max_val = augmented[row][col].abs();
                    max_row = row;
                }
            }

            // Skip column if all values are (near) zero
            if max_val < 1e-15 {
                continue;
            }

            // Swap rows if necessary
            if max_row != pivot_row {
                augmented.swap(pivot_row, max_row);
            }

            // Record pivot column
            pivot_cols.push(col);

            // Eliminate below
            let pivot_val = augmented[pivot_row][col];
            for row in (pivot_row + 1)..n_eqs {
                let factor = augmented[row][col] / pivot_val;
                augmented[row][col] = 0.0;
                for c in (col + 1)..=n_vars {
                    augmented[row][c] -= factor * augmented[pivot_row][c];
                }
            }

            pivot_row += 1;
        }

        let rank = pivot_cols.len();

        // Check for inconsistency: a row like [0 0 0 ... 0 | c] where c != 0
        for row in rank..n_eqs {
            let rhs = augmented[row][n_vars];
            let all_zeros = augmented[row][0..n_vars].iter().all(|&x| x.abs() < 1e-15);
            if all_zeros && rhs.abs() > 1e-15 {
                return Ok(SystemSolution::NoSolution);
            }
        }

        // Back substitution
        if rank == n_vars {
            // Unique solution
            let mut solution_values = vec![0.0; n_vars];

            // Back substitute from bottom to top
            for i in (0..rank).rev() {
                let col = pivot_cols[i];
                let mut sum = augmented[i][n_vars];
                for j in (col + 1)..n_vars {
                    sum -= augmented[i][j] * solution_values[j];
                }
                solution_values[col] = sum / augmented[i][col];
            }

            let mut result = HashMap::new();
            for (i, var) in self.variables.iter().enumerate() {
                let val = solution_values[i];
                // Convert to integer if close to an integer
                let expr = if (val - val.round()).abs() < 1e-10 {
                    Expression::Integer(val.round() as i64)
                } else {
                    Expression::Float(val)
                };
                result.insert(var.clone(), expr);
            }

            Ok(SystemSolution::Unique(result))
        } else {
            // Infinite solutions: some variables are free
            // Identify free variables (columns without pivots)
            let pivot_set: std::collections::HashSet<_> = pivot_cols.iter().cloned().collect();
            let free_cols: Vec<_> = (0..n_vars).filter(|c| !pivot_set.contains(c)).collect();

            let free_vars: Vec<_> = free_cols
                .iter()
                .map(|&c| self.variables[c].clone())
                .collect();

            // Back substitute to express bound variables in terms of free variables
            // For now, return a simplified result
            let mut bound = HashMap::new();

            // Back substitute for pivot variables
            for i in (0..rank).rev() {
                let col = pivot_cols[i];
                let rhs = augmented[i][n_vars];

                // Build expression: rhs - sum of (coefficient * free_var)
                let mut terms: Vec<Expression> = vec![];

                // Constant term
                if rhs.abs() > 1e-15 {
                    terms.push(if (rhs - rhs.round()).abs() < 1e-10 {
                        Expression::Integer(rhs.round() as i64)
                    } else {
                        Expression::Float(rhs)
                    });
                }

                // Free variable terms (negated because moved to RHS)
                for &free_col in &free_cols {
                    let coeff = -augmented[i][free_col] / augmented[i][col];
                    if coeff.abs() > 1e-15 {
                        let free_var = Expression::Variable(self.variables[free_col].clone());
                        let term = if (coeff - coeff.round()).abs() < 1e-10 {
                            let int_coeff = coeff.round() as i64;
                            if int_coeff == 1 {
                                free_var
                            } else if int_coeff == -1 {
                                Expression::Unary(UnaryOp::Neg, Box::new(free_var))
                            } else {
                                Expression::Binary(
                                    BinaryOp::Mul,
                                    Box::new(Expression::Integer(int_coeff)),
                                    Box::new(free_var),
                                )
                            }
                        } else {
                            Expression::Binary(
                                BinaryOp::Mul,
                                Box::new(Expression::Float(coeff)),
                                Box::new(free_var),
                            )
                        };
                        terms.push(term);
                    }
                }

                // Combine terms
                let expr = if terms.is_empty() {
                    Expression::Integer(0)
                } else if terms.len() == 1 {
                    terms.remove(0)
                } else {
                    let mut result = terms.remove(0);
                    for term in terms {
                        result =
                            Expression::Binary(BinaryOp::Add, Box::new(result), Box::new(term));
                    }
                    result
                };

                // Divide by pivot coefficient if not 1
                let pivot_coeff = augmented[i][col];
                let final_expr = if (pivot_coeff - 1.0).abs() < 1e-15 {
                    expr
                } else {
                    Expression::Binary(
                        BinaryOp::Div,
                        Box::new(expr),
                        Box::new(if (pivot_coeff - pivot_coeff.round()).abs() < 1e-10 {
                            Expression::Integer(pivot_coeff.round() as i64)
                        } else {
                            Expression::Float(pivot_coeff)
                        }),
                    )
                };

                bound.insert(self.variables[col].clone(), final_expr);
            }

            Ok(SystemSolution::Infinite {
                bound,
                free: free_vars,
            })
        }
    }

    /// Solve a 2x2 or 3x3 system using Cramer's rule.
    pub fn solve_cramers(&self) -> SolverResult<SystemSolution> {
        let n = self.variables.len();
        if self.coefficients.len() != n {
            return Err(SolverError::Other(
                "Cramer's rule requires square system".to_string(),
            ));
        }

        if n != 2 && n != 3 {
            return Err(SolverError::Other(
                "Cramer's rule only implemented for 2x2 and 3x3 systems".to_string(),
            ));
        }

        let det_a = if n == 2 {
            Self::det_2x2(&self.coefficients)
        } else {
            Self::det_3x3(&self.coefficients)
        };

        if det_a.abs() < 1e-15 {
            // Determinant is zero - system may have no or infinite solutions
            // Fall back to Gaussian elimination
            return self.solve();
        }

        let mut result = HashMap::new();

        for i in 0..n {
            // Create matrix with column i replaced by constants
            let mut modified: Vec<Vec<f64>> = self.coefficients.clone();
            for (row, &c) in self.constants.iter().enumerate() {
                modified[row][i] = c;
            }

            let det_i = if n == 2 {
                Self::det_2x2(&modified)
            } else {
                Self::det_3x3(&modified)
            };

            let val = det_i / det_a;
            let expr = if (val - val.round()).abs() < 1e-10 {
                Expression::Integer(val.round() as i64)
            } else {
                Expression::Float(val)
            };

            result.insert(self.variables[i].clone(), expr);
        }

        Ok(SystemSolution::Unique(result))
    }

    /// Compute 2x2 determinant.
    fn det_2x2(m: &[Vec<f64>]) -> f64 {
        m[0][0] * m[1][1] - m[0][1] * m[1][0]
    }

    /// Compute 3x3 determinant using expansion by first row.
    fn det_3x3(m: &[Vec<f64>]) -> f64 {
        let a = m[0][0];
        let b = m[0][1];
        let c = m[0][2];

        let minor1 = m[1][1] * m[2][2] - m[1][2] * m[2][1];
        let minor2 = m[1][0] * m[2][2] - m[1][2] * m[2][0];
        let minor3 = m[1][0] * m[2][1] - m[1][1] * m[2][0];

        a * minor1 - b * minor2 + c * minor3
    }
}

/// System of equations solver.
#[derive(Debug, Default)]
pub struct SystemSolver;

impl SystemSolver {
    /// Creates a new system of equations solver.
    pub fn new() -> Self {
        Self
    }

    /// Solve a system of linear equations for multiple variables.
    ///
    /// Uses Gaussian elimination with partial pivoting for general systems.
    /// For 2x2 and 3x3 systems, Cramer's rule is also available.
    ///
    /// # Arguments
    ///
    /// * `equations` - The equations to solve
    /// * `variables` - The variables to solve for
    ///
    /// # Returns
    ///
    /// * `SystemSolution::Unique` - If there is exactly one solution
    /// * `SystemSolution::Infinite` - If there are infinitely many solutions
    /// * `SystemSolution::NoSolution` - If the system is inconsistent
    ///
    /// # Examples
    ///
    /// ```
    /// use thales::solver::{SystemSolver, SystemSolution};
    /// use thales::ast::{Equation, Expression, Variable, BinaryOp};
    ///
    /// let solver = SystemSolver::new();
    ///
    /// // Solve: x + y = 5, x - y = 1
    /// let x = Variable::new("x");
    /// let y = Variable::new("y");
    ///
    /// let eq1 = Equation::new(
    ///     "eq1",
    ///     Expression::Binary(
    ///         BinaryOp::Add,
    ///         Box::new(Expression::Variable(x.clone())),
    ///         Box::new(Expression::Variable(y.clone())),
    ///     ),
    ///     Expression::Integer(5),
    /// );
    ///
    /// let eq2 = Equation::new(
    ///     "eq2",
    ///     Expression::Binary(
    ///         BinaryOp::Sub,
    ///         Box::new(Expression::Variable(x.clone())),
    ///         Box::new(Expression::Variable(y.clone())),
    ///     ),
    ///     Expression::Integer(1),
    /// );
    ///
    /// let result = solver.solve_linear_system(&[eq1, eq2], &[x.clone(), y.clone()]).unwrap();
    /// match result {
    ///     SystemSolution::Unique(sol) => {
    ///         // x = 3, y = 2
    ///         assert!(sol.contains_key(&x));
    ///         assert!(sol.contains_key(&y));
    ///     }
    ///     _ => panic!("Expected unique solution"),
    /// }
    /// ```
    pub fn solve_linear_system(
        &self,
        equations: &[Equation],
        variables: &[Variable],
    ) -> SolverResult<SystemSolution> {
        let system = LinearSystem::from_equations(equations, variables)?;
        system.solve()
    }

    /// Solve using Cramer's rule (2x2 and 3x3 systems only).
    pub fn solve_cramers(
        &self,
        equations: &[Equation],
        variables: &[Variable],
    ) -> SolverResult<SystemSolution> {
        let system = LinearSystem::from_equations(equations, variables)?;
        system.solve_cramers()
    }

    /// Solve a system of equations for multiple variables.
    ///
    /// This is a legacy method that delegates to solve_linear_system.
    pub fn solve_system(
        &self,
        equations: &[Equation],
        variables: &[Variable],
    ) -> SolverResult<HashMap<Variable, Solution>> {
        let result = self.solve_linear_system(equations, variables)?;

        match result {
            SystemSolution::Unique(sol) => {
                let mut out = HashMap::new();
                for (var, expr) in sol {
                    out.insert(var, Solution::Unique(expr));
                }
                Ok(out)
            }
            SystemSolution::Infinite { bound, free: _ } => {
                let mut out = HashMap::new();
                for (var, expr) in bound {
                    out.insert(
                        var,
                        Solution::Parametric {
                            expression: expr,
                            constraints: vec![],
                        },
                    );
                }
                Ok(out)
            }
            SystemSolution::NoSolution => Err(SolverError::NoSolution),
        }
    }
}

#[derive(Debug)]
/// Automatic solver dispatcher that selects the appropriate solving method.
///
/// `SmartSolver` examines the equation structure and dispatches to the most
/// suitable specialized solver. This eliminates the need to manually choose
/// between linear, quadratic, polynomial, or transcendental solvers.
///
/// # Priority Order
///
/// The solver tries methods in this priority order:
/// 1. **Linear** ([`LinearSolver`]): Fastest, handles equations like `ax + b = c`
/// 2. **Quadratic** ([`QuadraticSolver`]): Equations with x² terms
/// 3. **Polynomial** ([`PolynomialSolver`]): General polynomial equations
/// 4. **Transcendental** ([`TranscendentalSolver`]): Equations with sin, cos, exp, ln, log
///
/// This priority ensures simpler methods are tried first, falling back to more
/// complex methods only when needed.
///
/// # Examples
///
/// ## Linear Equation
///
/// ```
/// use thales::solver::{SmartSolver, Solver, Solution};
/// use thales::ast::{Equation, Expression, Variable, BinaryOp};
///
/// let solver = SmartSolver::new();
///
/// // Solve: 2x + 3 = 11
/// let x = Expression::Variable(Variable::new("x"));
/// let two_x = Expression::Binary(
///     BinaryOp::Mul,
///     Box::new(Expression::Integer(2)),
///     Box::new(x.clone()),
/// );
/// let left = Expression::Binary(
///     BinaryOp::Add,
///     Box::new(two_x),
///     Box::new(Expression::Integer(3)),
/// );
/// let equation = Equation::new("example", left, Expression::Integer(11));
///
/// let (solution, path) = solver.solve(&equation, &Variable::new("x")).unwrap();
///
/// // SmartSolver automatically selected LinearSolver
/// match solution {
///     Solution::Unique(expr) => {
///         assert_eq!(expr.evaluate(&std::collections::HashMap::new()), Some(4.0));
///     }
///     _ => panic!("Expected unique solution"),
/// }
/// ```
///
/// ## Transcendental Equation
///
/// ```no_run
/// use thales::solver::{SmartSolver, Solver, Solution};
/// use thales::ast::{Equation, Expression, Variable, Function};
///
/// let solver = SmartSolver::new();
///
/// // Solve: sin(x) = 0.5
/// let x = Expression::Variable(Variable::new("x"));
/// let sin_x = Expression::Function(Function::Sin, vec![x]);
/// let equation = Equation::new("trig", sin_x, Expression::Float(0.5));
///
/// let (solution, path) = solver.solve(&equation, &Variable::new("x")).unwrap();
///
/// // SmartSolver automatically selected TranscendentalSolver
/// match solution {
///     Solution::Unique(expr) => {
///         // expr contains asin(0.5)
///         let result = expr.evaluate(&std::collections::HashMap::new()).unwrap();
///         assert!((result - 0.5236).abs() < 0.001); // π/6 radians
///     }
///     _ => panic!("Expected unique solution"),
/// }
/// ```
///
/// ## Error Handling
///
/// ```
/// use thales::solver::{SmartSolver, Solver, SolverError};
/// use thales::ast::{Equation, Expression, Variable};
///
/// let solver = SmartSolver::new();
///
/// // Variable not in equation
/// let equation = Equation::new("bad", Expression::Integer(0), Expression::Integer(5));
/// let result = solver.solve(&equation, &Variable::new("x"));
///
/// // Since x doesn't appear in the equation, solver cannot handle it
/// assert!(result.is_err());
/// match result {
///     Err(SolverError::CannotSolve(_)) | Err(SolverError::UnsupportedEquationType) => {
///         // Expected - x not in equation or equation not supported
///     }
///     _ => panic!("Expected CannotSolve or UnsupportedEquationType error"),
/// }
/// ```
///
/// # Implementation Notes
///
/// - `SmartSolver` maintains instances of all specialized solvers internally
/// - The `can_solve()` check is fast and only examines equation structure
/// - If no solver can handle the equation, returns [`SolverError::UnsupportedEquationType`]
/// - The solving process is deterministic - same equation always uses same solver
///
/// # See Also
///
/// - [`solve_for`]: High-level API that uses `SmartSolver` and handles value substitution
/// - [`Solver`]: Base trait implemented by all solvers
/// - [`LinearSolver`], [`QuadraticSolver`], [`PolynomialSolver`], [`TranscendentalSolver`]: Specialized solvers
pub struct SmartSolver {
    linear: LinearSolver,
    quadratic: QuadraticSolver,
    polynomial: PolynomialSolver,
    transcendental: TranscendentalSolver,
}

impl SmartSolver {
    /// Creates a new smart solver with all specialized solvers initialized.
    pub fn new() -> Self {
        Self {
            linear: LinearSolver::new(),
            quadratic: QuadraticSolver::new(),
            polynomial: PolynomialSolver::new(),
            transcendental: TranscendentalSolver::new(),
        }
    }
}

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

impl Solver for SmartSolver {
    fn solve(
        &self,
        equation: &Equation,
        variable: &Variable,
    ) -> SolverResult<(Solution, ResolutionPath)> {
        // TODO: Analyze equation and dispatch to appropriate solver
        // Priority order: linear -> quadratic -> polynomial -> transcendental
        if self.linear.can_solve(equation) {
            self.linear.solve(equation, variable)
        } else if self.quadratic.can_solve(equation) {
            self.quadratic.solve(equation, variable)
        } else if self.polynomial.can_solve(equation) {
            self.polynomial.solve(equation, variable)
        } else if self.transcendental.can_solve(equation) {
            self.transcendental.solve(equation, variable)
        } else {
            Err(SolverError::UnsupportedEquationType)
        }
    }

    fn can_solve(&self, equation: &Equation) -> bool {
        self.linear.can_solve(equation)
            || self.quadratic.can_solve(equation)
            || self.polynomial.can_solve(equation)
            || self.transcendental.can_solve(equation)
    }
}

// ============================================================================
// High-Level API
// ============================================================================

/// Solve an equation for a specific variable with known values substituted.
///
/// This is the primary high-level API for equation solving. It combines symbolic
/// solving with numeric evaluation in three steps:
///
/// 1. **Symbolic solving**: Uses [`SmartSolver`] to solve for the target variable
/// 2. **Value substitution**: Replaces known variables with their numeric values
/// 3. **Simplification**: Evaluates constants and simplifies the result
///
/// This function is ideal when you have an equation with multiple variables and
/// want to solve for one variable given values for the others.
///
/// # Arguments
///
/// * `equation` - The equation to solve (e.g., `ax + b = c`)
/// * `target` - Name of the variable to solve for (e.g., `"x"`)
/// * `known_values` - HashMap mapping variable names to their numeric values
///
/// # Returns
///
/// A [`ResolutionPath`] containing:
/// - All solving steps performed (isolating variable, applying operations)
/// - The final result expression (fully evaluated if all values known)
/// - Operation history for display/debugging
///
/// # Errors
///
/// Returns [`SolverError`] if:
/// - [`SolverError::UnsupportedEquationType`]: No solver can handle this equation type
/// - [`SolverError::NoSolution`]: Equation is inconsistent (e.g., `0 = 5`)
/// - [`SolverError::InfiniteSolutions`]: Equation is an identity (e.g., `x = x`)
/// - [`SolverError::CannotSolve`]: Target variable not found or in unsolvable form
/// - [`SolverError::Other`]: Evaluation failed or unsupported solution type
///
/// # Examples
///
/// ## Basic Linear Equation with Substitution
///
/// ```
/// use thales::solver::solve_for;
/// use thales::ast::{Equation, Expression, Variable, BinaryOp};
/// use std::collections::HashMap;
///
/// // Solve: ax + b = c for x, given a=2, b=3, c=11
/// let a = Expression::Variable(Variable::new("a"));
/// let x = Expression::Variable(Variable::new("x"));
/// let b = Expression::Variable(Variable::new("b"));
/// let c = Expression::Variable(Variable::new("c"));
///
/// let ax = Expression::Binary(BinaryOp::Mul, Box::new(a), Box::new(x));
/// let left = Expression::Binary(BinaryOp::Add, Box::new(ax), Box::new(b));
/// let equation = Equation::new("linear", left, c);
///
/// let mut known = HashMap::new();
/// known.insert("a".to_string(), 2.0);
/// known.insert("b".to_string(), 3.0);
/// known.insert("c".to_string(), 11.0);
///
/// let path = solve_for(&equation, "x", &known).unwrap();
///
/// // Result is x = 4.0
/// assert_eq!(path.result.evaluate(&HashMap::new()), Some(4.0));
///
/// // Path contains solving steps
/// assert!(!path.steps.is_empty());
/// ```
///
/// ## Physics Formula: Ohm's Law
///
/// ```
/// use thales::solver::solve_for;
/// use thales::ast::{Equation, Expression, Variable, BinaryOp};
/// use std::collections::HashMap;
///
/// // V = I * R, solve for I given V=12V, R=4Ω
/// let v = Expression::Variable(Variable::new("V"));
/// let i = Expression::Variable(Variable::new("I"));
/// let r = Expression::Variable(Variable::new("R"));
///
/// let ir = Expression::Binary(BinaryOp::Mul, Box::new(i), Box::new(r));
/// let equation = Equation::new("ohms_law", v, ir);
///
/// let mut known = HashMap::new();
/// known.insert("V".to_string(), 12.0);
/// known.insert("R".to_string(), 4.0);
///
/// let path = solve_for(&equation, "I", &known).unwrap();
///
/// // I = V/R = 12/4 = 3A
/// assert_eq!(path.result.evaluate(&HashMap::new()), Some(3.0));
/// ```
///
/// ## Symbolic Solution (No Known Values)
///
/// ```
/// use thales::solver::solve_for;
/// use thales::ast::{Equation, Expression, Variable, BinaryOp};
/// use std::collections::HashMap;
///
/// // Solve: 2x + 3 = y for x (no known values)
/// let x = Expression::Variable(Variable::new("x"));
/// let y = Expression::Variable(Variable::new("y"));
/// let two_x = Expression::Binary(
///     BinaryOp::Mul,
///     Box::new(Expression::Integer(2)),
///     Box::new(x),
/// );
/// let left = Expression::Binary(
///     BinaryOp::Add,
///     Box::new(two_x),
///     Box::new(Expression::Integer(3)),
/// );
/// let equation = Equation::new("symbolic", left, y);
///
/// let known = HashMap::new(); // No known values
///
/// let path = solve_for(&equation, "x", &known).unwrap();
///
/// // Result is symbolic: x = (y - 3) / 2
/// // Can evaluate when y is provided later
/// let mut eval_context = HashMap::new();
/// eval_context.insert("y".to_string(), 11.0);
/// assert_eq!(path.result.evaluate(&eval_context), Some(4.0));
/// ```
///
/// ## Error Handling
///
/// ```
/// use thales::solver::{solve_for, SolverError};
/// use thales::ast::{Equation, Expression};
/// use std::collections::HashMap;
///
/// // Variable not in equation: 0 = 5
/// let equation = Equation::new("bad", Expression::Integer(0), Expression::Integer(5));
/// let known = HashMap::new();
///
/// let result = solve_for(&equation, "x", &known);
///
/// // Should fail because x doesn't appear in equation
/// assert!(result.is_err());
/// match result {
///     Err(SolverError::CannotSolve(_)) | Err(SolverError::UnsupportedEquationType) => {
///         // Expected - x not in equation
///     }
///     _ => panic!("Should have failed with CannotSolve or UnsupportedEquationType"),
/// }
/// ```
///
/// # Parameter Details
///
/// ## Variable Substitution
///
/// The `known_values` HashMap substitutes variables before final evaluation:
///
/// ```text
/// Equation: ax + b = c
/// Target: x
/// Known: {a: 2, b: 3, c: 11}
///
/// Step 1 (Symbolic): x = (c - b) / a
/// Step 2 (Substitute): x = (11 - 3) / 2
/// Step 3 (Evaluate): x = 4
/// ```
///
/// Variables not in `known_values` remain symbolic in the result.
///
/// ## Result Expression
///
/// The `ResolutionPath.result` field contains the final expression:
/// - If all variables known: evaluates to a single number
/// - If some variables unknown: contains symbolic expression
/// - Use [`Expression::evaluate`] with additional values to compute final result
///
/// # Implementation Notes
///
/// - Uses [`SmartSolver`] internally for automatic method selection
/// - Only supports unique solutions currently (not multiple/parametric)
/// - Simplification is automatic - no manual step required
/// - Thread-safe: can be called concurrently on different equations
///
/// # See Also
///
/// - [`SmartSolver`]: The underlying solver used for symbolic solving
/// - [`ResolutionPath`]: Return type containing steps and result
/// - [`compute_partial_derivative`]: For derivative computation (uncertainty propagation)
/// - [`Expression::evaluate`]: For evaluating symbolic results with values
pub fn solve_for(
    equation: &Equation,
    target: &str,
    known_values: &HashMap<String, f64>,
) -> Result<ResolutionPath, SolverError> {
    // Create Variable from target string
    let target_var = Variable::new(target);

    // Try solving with SmartSolver
    let solver = SmartSolver::new();
    let (solution, mut path) = solver.solve(equation, &target_var)?;

    // Extract the solution expression
    let solution_expr = match solution {
        Solution::Unique(expr) => expr,
        Solution::Multiple(_) => {
            return Err(SolverError::Other(
                "Multiple solutions not yet supported in solve_for".to_string(),
            ))
        }
        Solution::None => return Err(SolverError::NoSolution),
        Solution::Infinite => return Err(SolverError::InfiniteSolutions),
        Solution::Parametric { .. } => {
            return Err(SolverError::Other(
                "Parametric solutions not yet supported in solve_for".to_string(),
            ))
        }
    };

    // Substitute known values
    if !known_values.is_empty() {
        let substituted = substitute_values(&solution_expr, known_values);
        let simplified = substituted.simplify();
        let evaluated = evaluate_constants(&simplified);

        path.add_step(ResolutionStep::new(
            Operation::Substitute {
                variable: Variable::new("known_values"),
                value: Expression::Integer(0), // Placeholder
            },
            "Substitute known values and evaluate".to_string(),
            evaluated.clone(),
        ));

        path.set_result(evaluated);
    } else {
        path.set_result(solution_expr);
    }

    Ok(path)
}

/// Substitute known variable values into an expression.
fn substitute_values(expr: &Expression, values: &HashMap<String, f64>) -> Expression {
    match expr {
        Expression::Variable(v) => {
            if let Some(&value) = values.get(&v.name) {
                Expression::Float(value)
            } else {
                expr.clone()
            }
        }
        Expression::Unary(op, inner) => {
            Expression::Unary(*op, Box::new(substitute_values(inner, values)))
        }
        Expression::Binary(op, left, right) => Expression::Binary(
            *op,
            Box::new(substitute_values(left, values)),
            Box::new(substitute_values(right, values)),
        ),
        Expression::Function(func, args) => Expression::Function(
            func.clone(),
            args.iter()
                .map(|arg| substitute_values(arg, values))
                .collect(),
        ),
        Expression::Power(base, exp) => Expression::Power(
            Box::new(substitute_values(base, values)),
            Box::new(substitute_values(exp, values)),
        ),
        _ => expr.clone(),
    }
}

/// Compute a partial derivative for uncertainty propagation and sensitivity analysis.
///
/// Given an equation defining an output variable in terms of input variables,
/// this function computes the partial derivative ∂output/∂input and evaluates
/// it at the given point. This is essential for:
///
/// - **Uncertainty propagation**: Computing how measurement errors affect results
/// - **Sensitivity analysis**: Understanding which inputs most affect the output
/// - **Error bars**: Calculating confidence intervals for computed values
///
/// # Mathematical Background
///
/// For a function V(l, w, h), the partial derivative ∂V/∂l represents the rate
/// of change of V with respect to l while holding w and h constant.
///
/// **Example**: Box volume V = l × w × h
/// ```text
/// ∂V/∂l = w × h    (derivative with respect to length)
/// ∂V/∂w = l × h    (derivative with respect to width)
/// ∂V/∂h = l × w    (derivative with respect to height)
/// ```
///
/// # Uncertainty Propagation Formula
///
/// For independent input variables with uncertainties δx₁, δx₂, ..., the
/// uncertainty in the output δy is:
///
/// ```text
/// δy = √[(∂y/∂x₁ · δx₁)² + (∂y/∂x₂ · δx₂)² + ...]
/// ```
///
/// Use [`compute_all_partial_derivatives`] to get all derivatives at once.
///
/// # Arguments
///
/// * `equation` - Equation with output variable isolated (e.g., `V = l * w * h`)
/// * `output_var` - Name of the output variable (left or right side, e.g., `"V"`)
/// * `input_var` - Name of the input variable to differentiate with respect to (e.g., `"l"`)
/// * `values` - HashMap of all variable values at the evaluation point
///
/// # Returns
///
/// The numerical value of ∂output/∂input evaluated at the given point.
///
/// # Errors
///
/// Returns [`SolverError`] if:
/// - [`SolverError::CannotSolve`]: Output variable not found or not isolated in equation
/// - [`SolverError::Other`]: Evaluation failed due to missing values in the HashMap
///
/// # Examples
///
/// ## Box Volume Derivative
///
/// ```
/// use thales::ast::{Equation, Expression, Variable, BinaryOp};
/// use thales::solver::compute_partial_derivative;
/// use std::collections::HashMap;
///
/// // Equation: V = l * w * h
/// let l = Expression::Variable(Variable::new("l"));
/// let w = Expression::Variable(Variable::new("w"));
/// let h = Expression::Variable(Variable::new("h"));
/// let v = Expression::Variable(Variable::new("V"));
///
/// let lw = Expression::Binary(BinaryOp::Mul, Box::new(l), Box::new(w));
/// let lwh = Expression::Binary(BinaryOp::Mul, Box::new(lw), Box::new(h));
/// let equation = Equation::new("box_volume", v, lwh);
///
/// let mut values = HashMap::new();
/// values.insert("l".to_string(), 2.0);
/// values.insert("w".to_string(), 3.0);
/// values.insert("h".to_string(), 4.0);
///
/// // ∂V/∂l = w * h = 3 * 4 = 12.0
/// let dv_dl = compute_partial_derivative(&equation, "V", "l", &values).unwrap();
/// assert_eq!(dv_dl, 12.0);
///
/// // ∂V/∂w = l * h = 2 * 4 = 8.0
/// let dv_dw = compute_partial_derivative(&equation, "V", "w", &values).unwrap();
/// assert_eq!(dv_dw, 8.0);
///
/// // ∂V/∂h = l * w = 2 * 3 = 6.0
/// let dv_dh = compute_partial_derivative(&equation, "V", "h", &values).unwrap();
/// assert_eq!(dv_dh, 6.0);
/// ```
///
/// ## Ohm's Law Sensitivity
///
/// ```
/// use thales::ast::{Equation, Expression, Variable, BinaryOp};
/// use thales::solver::compute_partial_derivative;
/// use std::collections::HashMap;
///
/// // P = V² / R (power dissipation)
/// let v = Expression::Variable(Variable::new("V"));
/// let r = Expression::Variable(Variable::new("R"));
/// let p = Expression::Variable(Variable::new("P"));
///
/// let v_squared = Expression::Power(Box::new(v), Box::new(Expression::Integer(2)));
/// let power = Expression::Binary(BinaryOp::Div, Box::new(v_squared), Box::new(r));
/// let equation = Equation::new("power", p, power);
///
/// let mut values = HashMap::new();
/// values.insert("V".to_string(), 12.0);
/// values.insert("R".to_string(), 4.0);
///
/// // ∂P/∂V = 2V/R = 2(12)/4 = 6.0 W/V
/// let dp_dv = compute_partial_derivative(&equation, "P", "V", &values).unwrap();
/// assert_eq!(dp_dv, 6.0);
///
/// // This means a 1V change in voltage causes ~6W change in power
/// ```
///
/// ## Uncertainty Propagation Example
///
/// ```
/// use thales::ast::{Equation, Expression, Variable, BinaryOp};
/// use thales::solver::compute_partial_derivative;
/// use std::collections::HashMap;
///
/// // Area = π * r² with r = 5.0 ± 0.1 cm
/// let pi = std::f64::consts::PI;
/// let r = Expression::Variable(Variable::new("r"));
/// let a = Expression::Variable(Variable::new("A"));
///
/// let r_squared = Expression::Power(Box::new(r), Box::new(Expression::Integer(2)));
/// let pi_r_sq = Expression::Binary(
///     BinaryOp::Mul,
///     Box::new(Expression::Float(pi)),
///     Box::new(r_squared),
/// );
/// let equation = Equation::new("circle_area", a, pi_r_sq);
///
/// let mut values = HashMap::new();
/// values.insert("r".to_string(), 5.0);
///
/// // ∂A/∂r = 2πr
/// let da_dr = compute_partial_derivative(&equation, "A", "r", &values).unwrap();
/// assert!((da_dr - 2.0 * pi * 5.0).abs() < 0.001);
///
/// // Uncertainty: δA = |∂A/∂r| × δr = 31.416 × 0.1 ≈ 3.14 cm²
/// let delta_r = 0.1;
/// let delta_a = da_dr * delta_r;
/// assert!((delta_a - 3.14159).abs() < 0.001);
/// ```
///
/// # Implementation Notes
///
/// The function:
/// 1. Extracts the expression for the output variable from the equation
/// 2. Computes the symbolic derivative using [`Expression::differentiate`]
/// 3. Simplifies the derivative expression
/// 4. Evaluates it numerically with the provided values
///
/// **Requirement**: The equation must have the output variable isolated on one side.
/// For example, `V = l*w*h` is valid, but `V - l*w*h = 0` is not (yet).
///
/// # See Also
///
/// - [`compute_all_partial_derivatives`]: Compute all partial derivatives at once
/// - [`Expression::differentiate`]: Underlying symbolic differentiation
/// - [`solve_for`]: Solve equations before computing derivatives
pub fn compute_partial_derivative(
    equation: &Equation,
    output_var: &str,
    input_var: &str,
    values: &HashMap<String, f64>,
) -> Result<f64, SolverError> {
    // Get the expression for the output variable
    // If equation is "output = expr", use expr
    // If equation is "expr = output", use expr
    let output_expr = if let Expression::Variable(v) = &equation.left {
        if v.name == output_var {
            &equation.right
        } else if let Expression::Variable(v2) = &equation.right {
            if v2.name == output_var {
                &equation.left
            } else {
                return Err(SolverError::CannotSolve(format!(
                    "Output variable '{}' not found in equation",
                    output_var
                )));
            }
        } else {
            return Err(SolverError::CannotSolve(format!(
                "Output variable '{}' not found in equation",
                output_var
            )));
        }
    } else if let Expression::Variable(v) = &equation.right {
        if v.name == output_var {
            &equation.left
        } else {
            return Err(SolverError::CannotSolve(format!(
                "Output variable '{}' not found in equation",
                output_var
            )));
        }
    } else {
        return Err(SolverError::CannotSolve(
            "Equation does not have output variable isolated".to_string(),
        ));
    };

    // Compute the derivative symbolically
    let derivative_expr = output_expr.differentiate(input_var);

    // Simplify the derivative
    let simplified = derivative_expr.simplify();

    // Evaluate the derivative at the given values
    simplified.evaluate(values).ok_or_else(|| {
        SolverError::Other(format!(
            "Failed to evaluate derivative - missing or invalid values"
        ))
    })
}

/// Compute all partial derivatives for complete uncertainty propagation.
///
/// Given an equation defining an output variable in terms of input variables,
/// this function computes all partial derivatives ∂output/∂input_i and evaluates
/// them at the given point. This is the recommended way to compute derivatives
/// when you need multiple partial derivatives, as it provides a clean interface
/// for uncertainty propagation.
///
/// # Uncertainty Propagation
///
/// This function provides all the partial derivatives needed to compute output
/// uncertainty from input uncertainties using the standard formula:
///
/// ```text
/// δy = √[(∂y/∂x₁ · δx₁)² + (∂y/∂x₂ · δx₂)² + ... + (∂y/∂xₙ · δxₙ)²]
/// ```
///
/// where:
/// - `δy` is the uncertainty in the output
/// - `∂y/∂xᵢ` are the partial derivatives (computed by this function)
/// - `δxᵢ` are the uncertainties in the input variables
///
/// # Arguments
///
/// * `equation` - Equation with output variable isolated (e.g., `V = l * w * h`)
/// * `output_var` - Name of the output variable (e.g., `"V"`)
/// * `input_vars` - Slice of input variable names to compute derivatives for
/// * `values` - HashMap of all variable values at the evaluation point
///
/// # Returns
///
/// A HashMap mapping each input variable name to its partial derivative value.
/// The returned HashMap has the same keys as provided in `input_vars`.
///
/// # Errors
///
/// Returns [`SolverError`] if any partial derivative computation fails:
/// - [`SolverError::CannotSolve`]: Output variable not found or not isolated
/// - [`SolverError::Other`]: Evaluation failed due to missing values
///
/// # Examples
///
/// ## Box Volume with Complete Uncertainty
///
/// ```
/// use thales::ast::{Equation, Expression, Variable, BinaryOp};
/// use thales::solver::compute_all_partial_derivatives;
/// use std::collections::HashMap;
///
/// // Equation: V = l * w * h
/// let l = Expression::Variable(Variable::new("l"));
/// let w = Expression::Variable(Variable::new("w"));
/// let h = Expression::Variable(Variable::new("h"));
/// let v = Expression::Variable(Variable::new("V"));
///
/// let lw = Expression::Binary(BinaryOp::Mul, Box::new(l), Box::new(w));
/// let lwh = Expression::Binary(BinaryOp::Mul, Box::new(lw), Box::new(h));
/// let equation = Equation::new("box_volume", v, lwh);
///
/// // Measurements with uncertainties
/// let mut values = HashMap::new();
/// values.insert("l".to_string(), 2.0);  // 2.0 ± 0.1 cm
/// values.insert("w".to_string(), 3.0);  // 3.0 ± 0.1 cm
/// values.insert("h".to_string(), 4.0);  // 4.0 ± 0.1 cm
///
/// let input_vars = vec!["l".to_string(), "w".to_string(), "h".to_string()];
/// let derivatives = compute_all_partial_derivatives(
///     &equation,
///     "V",
///     &input_vars,
///     &values
/// ).unwrap();
///
/// // Verify partial derivatives
/// assert_eq!(derivatives.get("l").unwrap(), &12.0); // ∂V/∂l = w*h = 3*4
/// assert_eq!(derivatives.get("w").unwrap(), &8.0);  // ∂V/∂w = l*h = 2*4
/// assert_eq!(derivatives.get("h").unwrap(), &6.0);  // ∂V/∂h = l*w = 2*3
///
/// // Compute uncertainty: δV = √[(∂V/∂l·δl)² + (∂V/∂w·δw)² + (∂V/∂h·δh)²]
/// let delta_l = 0.1;
/// let delta_w = 0.1;
/// let delta_h = 0.1;
///
/// let delta_v = (
///     (derivatives["l"] * delta_l).powi(2) +
///     (derivatives["w"] * delta_w).powi(2) +
///     (derivatives["h"] * delta_h).powi(2)
/// ).sqrt();
///
/// // V = 24.0 ± 1.56 cm³
/// assert!((delta_v - 1.562).abs() < 0.01);
/// ```
///
/// ## Slide Rule Calculation with Error Propagation
///
/// ```
/// use thales::ast::{Equation, Expression, Variable, BinaryOp};
/// use thales::solver::compute_all_partial_derivatives;
/// use std::collections::HashMap;
///
/// // Power: P = V * I with measurement errors
/// let v = Expression::Variable(Variable::new("V"));
/// let i = Expression::Variable(Variable::new("I"));
/// let p = Expression::Variable(Variable::new("P"));
///
/// let power = Expression::Binary(BinaryOp::Mul, Box::new(v), Box::new(i));
/// let equation = Equation::new("power", p, power);
///
/// // Slide rule readings (±2% typical accuracy)
/// let mut values = HashMap::new();
/// values.insert("V".to_string(), 12.0);  // 12.0V ± 2%
/// values.insert("I".to_string(), 2.5);   // 2.5A ± 2%
///
/// let input_vars = vec!["V".to_string(), "I".to_string()];
/// let derivatives = compute_all_partial_derivatives(
///     &equation,
///     "P",
///     &input_vars,
///     &values
/// ).unwrap();
///
/// // ∂P/∂V = I = 2.5 W/V
/// assert_eq!(derivatives["V"], 2.5);
///
/// // ∂P/∂I = V = 12.0 W/A
/// assert_eq!(derivatives["I"], 12.0);
///
/// // Combined uncertainty (2% error on each measurement)
/// let delta_v = 12.0 * 0.02;  // 0.24V
/// let delta_i = 2.5 * 0.02;   // 0.05A
///
/// let delta_p = (
///     (derivatives["V"] * delta_v).powi(2) +
///     (derivatives["I"] * delta_i).powi(2)
/// ).sqrt();
///
/// // P = 30W ± 0.84W (2.8% combined error)
/// assert!((delta_p - 0.84).abs() < 0.01);
/// ```
///
/// ## Sensitivity Analysis
///
/// ```
/// use thales::ast::{Equation, Expression, Variable, BinaryOp};
/// use thales::solver::compute_all_partial_derivatives;
/// use std::collections::HashMap;
///
/// // Area = length * width (simple formula for sensitivity analysis)
/// let l = Expression::Variable(Variable::new("L"));
/// let w = Expression::Variable(Variable::new("W"));
/// let a = Expression::Variable(Variable::new("A"));
///
/// let area = Expression::Binary(BinaryOp::Mul, Box::new(l), Box::new(w));
/// let equation = Equation::new("area", a, area);
///
/// let mut values = HashMap::new();
/// values.insert("L".to_string(), 10.0);  // Length (m)
/// values.insert("W".to_string(), 5.0);   // Width (m)
///
/// let input_vars = vec!["L".to_string(), "W".to_string()];
/// let derivatives = compute_all_partial_derivatives(
///     &equation,
///     "A",
///     &input_vars,
///     &values
/// ).unwrap();
///
/// // ∂A/∂L = W = 5.0
/// assert_eq!(derivatives["L"], 5.0);
///
/// // ∂A/∂W = L = 10.0
/// assert_eq!(derivatives["W"], 10.0);
///
/// // Length has stronger effect (larger derivative value)
/// // ∂A/∂W = 10.0 is twice ∂A/∂L = 5.0
/// assert!(derivatives["W"] > derivatives["L"]);
/// ```
///
/// # Performance Notes
///
/// This function calls [`compute_partial_derivative`] for each input variable.
/// If you only need one or two derivatives, calling [`compute_partial_derivative`]
/// directly may be more efficient.
///
/// # See Also
///
/// - [`compute_partial_derivative`]: Compute a single partial derivative
/// - [`Expression::differentiate`]: Underlying symbolic differentiation
/// - [`solve_for`]: Solve equations before computing derivatives
pub fn compute_all_partial_derivatives(
    equation: &Equation,
    output_var: &str,
    input_vars: &[String],
    values: &HashMap<String, f64>,
) -> Result<HashMap<String, f64>, SolverError> {
    let mut derivatives = HashMap::new();

    for input_var in input_vars {
        let derivative = compute_partial_derivative(equation, output_var, input_var, values)?;
        derivatives.insert(input_var.clone(), derivative);
    }

    Ok(derivatives)
}

// TODO: Add equation simplification before solving
// TODO: Add symbolic manipulation utilities
// TODO: Add support for inequalities
// TODO: Add support for absolute value equations
// TODO: Add support for piecewise functions
// TODO: Add step-by-step explanation generation

#[cfg(test)]
mod system_solver_tests {
    use super::*;
    use crate::ast::{BinaryOp, Equation, Expression, Variable};

    fn var(name: &str) -> Expression {
        Expression::Variable(Variable::new(name))
    }

    fn int(n: i64) -> Expression {
        Expression::Integer(n)
    }

    fn add(left: Expression, right: Expression) -> Expression {
        Expression::Binary(BinaryOp::Add, Box::new(left), Box::new(right))
    }

    fn sub(left: Expression, right: Expression) -> Expression {
        Expression::Binary(BinaryOp::Sub, Box::new(left), Box::new(right))
    }

    fn mul(left: Expression, right: Expression) -> Expression {
        Expression::Binary(BinaryOp::Mul, Box::new(left), Box::new(right))
    }

    #[test]
    fn test_2x2_unique_solution() {
        // Solve: x + y = 5, x - y = 1
        // Solution: x = 3, y = 2
        let x = Variable::new("x");
        let y = Variable::new("y");

        let eq1 = Equation::new("eq1", add(var("x"), var("y")), int(5));
        let eq2 = Equation::new("eq2", sub(var("x"), var("y")), int(1));

        let solver = SystemSolver::new();
        let result = solver
            .solve_linear_system(&[eq1, eq2], &[x.clone(), y.clone()])
            .unwrap();

        match result {
            SystemSolution::Unique(sol) => {
                let x_val = sol.get(&x).unwrap();
                let y_val = sol.get(&y).unwrap();

                let empty: HashMap<String, f64> = HashMap::new();
                assert_eq!(x_val.evaluate(&empty), Some(3.0));
                assert_eq!(y_val.evaluate(&empty), Some(2.0));
            }
            _ => panic!("Expected unique solution"),
        }
    }

    #[test]
    fn test_2x2_with_coefficients() {
        // Solve: 2x + 3y = 8, 4x - y = 2
        // Solution: x = 1, y = 2
        let x = Variable::new("x");
        let y = Variable::new("y");

        let eq1 = Equation::new(
            "eq1",
            add(mul(int(2), var("x")), mul(int(3), var("y"))),
            int(8),
        );
        let eq2 = Equation::new("eq2", sub(mul(int(4), var("x")), var("y")), int(2));

        let solver = SystemSolver::new();
        let result = solver
            .solve_linear_system(&[eq1, eq2], &[x.clone(), y.clone()])
            .unwrap();

        match result {
            SystemSolution::Unique(sol) => {
                let x_val = sol.get(&x).unwrap();
                let y_val = sol.get(&y).unwrap();

                let empty: HashMap<String, f64> = HashMap::new();
                assert_eq!(x_val.evaluate(&empty), Some(1.0));
                assert_eq!(y_val.evaluate(&empty), Some(2.0));
            }
            _ => panic!("Expected unique solution"),
        }
    }

    #[test]
    fn test_3x3_unique_solution() {
        // Solve: x + y + z = 6, 2x + y - z = 1, x - y + 2z = 5
        // Solution: x = 1, y = 2, z = 3
        let x = Variable::new("x");
        let y = Variable::new("y");
        let z = Variable::new("z");

        let eq1 = Equation::new("eq1", add(add(var("x"), var("y")), var("z")), int(6));
        let eq2 = Equation::new(
            "eq2",
            sub(add(mul(int(2), var("x")), var("y")), var("z")),
            int(1),
        );
        let eq3 = Equation::new(
            "eq3",
            add(sub(var("x"), var("y")), mul(int(2), var("z"))),
            int(5),
        );

        let solver = SystemSolver::new();
        let result = solver
            .solve_linear_system(&[eq1, eq2, eq3], &[x.clone(), y.clone(), z.clone()])
            .unwrap();

        match result {
            SystemSolution::Unique(sol) => {
                let empty: HashMap<String, f64> = HashMap::new();
                assert_eq!(sol.get(&x).unwrap().evaluate(&empty), Some(1.0));
                assert_eq!(sol.get(&y).unwrap().evaluate(&empty), Some(2.0));
                assert_eq!(sol.get(&z).unwrap().evaluate(&empty), Some(3.0));
            }
            _ => panic!("Expected unique solution"),
        }
    }

    #[test]
    fn test_underdetermined_system() {
        // Solve: x + y = 5 (one equation, two unknowns)
        // Should have infinite solutions
        let x = Variable::new("x");
        let y = Variable::new("y");

        let eq1 = Equation::new("eq1", add(var("x"), var("y")), int(5));

        let solver = SystemSolver::new();
        let result = solver
            .solve_linear_system(&[eq1], &[x.clone(), y.clone()])
            .unwrap();

        match result {
            SystemSolution::Infinite { bound, free } => {
                // One variable should be free
                assert!(!free.is_empty());
                // The other should be bound to an expression
                assert!(!bound.is_empty());
            }
            _ => panic!("Expected infinite solutions"),
        }
    }

    #[test]
    fn test_inconsistent_system() {
        // Solve: x + y = 5, x + y = 6 (no solution)
        let x = Variable::new("x");
        let y = Variable::new("y");

        let eq1 = Equation::new("eq1", add(var("x"), var("y")), int(5));
        let eq2 = Equation::new("eq2", add(var("x"), var("y")), int(6));

        let solver = SystemSolver::new();
        let result = solver
            .solve_linear_system(&[eq1, eq2], &[x.clone(), y.clone()])
            .unwrap();

        assert!(matches!(result, SystemSolution::NoSolution));
    }

    #[test]
    fn test_cramers_rule_2x2() {
        // Same as test_2x2_unique_solution but using Cramer's rule
        let x = Variable::new("x");
        let y = Variable::new("y");

        let eq1 = Equation::new("eq1", add(var("x"), var("y")), int(5));
        let eq2 = Equation::new("eq2", sub(var("x"), var("y")), int(1));

        let solver = SystemSolver::new();
        let result = solver
            .solve_cramers(&[eq1, eq2], &[x.clone(), y.clone()])
            .unwrap();

        match result {
            SystemSolution::Unique(sol) => {
                let empty: HashMap<String, f64> = HashMap::new();
                assert_eq!(sol.get(&x).unwrap().evaluate(&empty), Some(3.0));
                assert_eq!(sol.get(&y).unwrap().evaluate(&empty), Some(2.0));
            }
            _ => panic!("Expected unique solution"),
        }
    }

    #[test]
    fn test_linear_system_struct() {
        // Test LinearSystem::from_equations
        let x = Variable::new("x");
        let y = Variable::new("y");

        let eq1 = Equation::new("eq1", add(var("x"), var("y")), int(5));
        let eq2 = Equation::new("eq2", sub(var("x"), var("y")), int(1));

        let system = LinearSystem::from_equations(&[eq1, eq2], &[x.clone(), y.clone()]).unwrap();

        // Verify coefficients: x + y = 5 -> [1, 1 | 5], x - y = 1 -> [1, -1 | 1]
        assert_eq!(system.coefficients.len(), 2);
        assert_eq!(system.constants.len(), 2);
    }

    #[test]
    fn test_overdetermined_consistent() {
        // Solve: x + y = 5, x - y = 1, 2x = 6
        // All three equations are consistent with x = 3, y = 2
        let x = Variable::new("x");
        let y = Variable::new("y");

        let eq1 = Equation::new("eq1", add(var("x"), var("y")), int(5));
        let eq2 = Equation::new("eq2", sub(var("x"), var("y")), int(1));
        let eq3 = Equation::new("eq3", mul(int(2), var("x")), int(6));

        let solver = SystemSolver::new();
        let result = solver
            .solve_linear_system(&[eq1, eq2, eq3], &[x.clone(), y.clone()])
            .unwrap();

        match result {
            SystemSolution::Unique(sol) => {
                let empty: HashMap<String, f64> = HashMap::new();
                assert_eq!(sol.get(&x).unwrap().evaluate(&empty), Some(3.0));
                assert_eq!(sol.get(&y).unwrap().evaluate(&empty), Some(2.0));
            }
            _ => panic!("Expected unique solution"),
        }
    }
}