RustedSciThe 0.4.15

//! # Symbolic Engine Module
//!
//! This module provides a comprehensive symbolic mathematics engine for creating, manipulating,
//! and evaluating symbolic expressions. It serves as the core foundation for symbolic computation
//! in the RustedSciThe framework, enabling analytical differentiation, expression simplification,
//! and conversion to executable functions.
//!
//! ## Purpose
//!
//! The symbolic engine allows users to:
//! - Parse and create symbolic mathematical expressions
//! - Perform analytical differentiation and integration
//! - Simplify complex expressions algebraically
//! - Convert symbolic expressions to executable Rust functions (lambdification)
//! - Handle indexed variables for matrix/vector operations
//! - Support trigonometric, exponential, and logarithmic functions
//!
//! ## Main Structures and Methods
//!
//! ### `Expr` Enum
//! The core symbolic expression type supporting:
//! - **Variables**: `Var(String)` - symbolic variables like "x", "y"
//! - **Constants**: `Const(f64)` - numerical constants
//! - **Operations**: `Add`, `Sub`, `Mul`, `Div`, `Pow` - basic arithmetic
//! - **Functions**: `Exp`, `Ln`, `sin`, `cos`, etc. - mathematical functions
//!
//! ### Key Methods
//! - `Symbols(symbols: &str)` - Create multiple variables from comma-separated string
//! - `IndexedVar(index: usize, var_name: &str)` - Create indexed variables (x0, x1, etc.)
//! - `diff(var: &str)` - Analytical differentiation
//! - `lambdify()` - Convert to executable function
//! - `simplify_()` - Algebraic simplification
//! - `set_variable()` - Substitute variables with values
//!
//! ## Interesting Code Features
//!
//! 1. **Recursive Expression Tree**: Uses Box<Expr> for nested expressions, enabling
//!    arbitrarily complex mathematical structures
//!
//! 2. **Operator Overloading**: Implements std::ops traits (Add, Sub, Mul, Div) for
//!    natural mathematical syntax: `x + y * z`
//!
//! 3. **Pattern Matching Differentiation**: Uses exhaustive match statements to implement
//!    analytical differentiation rules (product rule, chain rule, etc.)
//!
//! 4. **Smart Simplification**: The `simplify_()` method applies algebraic rules like
//!    `x + 0 = x`, `x * 1 = x`, `0 * x = 0` recursively
//!
//! 5. **Indexed Variable System**: Supports both 1D (x0, x1) and 2D (A_2_3) indexing
//!    for matrix operations and large systems
//!
//! 6. **Macro System**: Provides convenient macros like `symbols!(x, y, z)` for
//!    ergonomic variable creation
//!
//! 7. **Non-standard Function Names**: Uses mathematical notation (tg, ctg) instead
//!    of programming conventions (tan, cot) for trigonometric functions

#![allow(non_camel_case_types)]

use crate::symbolic::symbolic_metadata::TraversalCache;
use std::collections::HashMap;
use std::f64;

use std::fmt;

/// Core symbolic expression enum representing mathematical expressions as an abstract syntax tree.
///
/// Each variant represents a different type of mathematical construct, from simple variables
/// and constants to complex nested operations. The enum uses Box<Expr> for recursive structures,
/// allowing arbitrarily deep expression trees.
///
/// # Examples
/// ```rust, ignore
/// use symbolic_engine::Expr;
/// let x = Expr::Var("x".to_string());
/// let expr = Expr::Add(Box::new(x), Box::new(Expr::Const(2.0)));
/// ```
#[derive(Clone, Debug, PartialEq)]
pub enum Expr {
    /// Symbolic variable with a name (e.g., "x", "y", "velocity")
    Var(String),
    /// Numerical constant value
    Const(f64),
    /// Addition operation: left + right
    Add(Box<Expr>, Box<Expr>),
    /// Subtraction operation: left - right
    Sub(Box<Expr>, Box<Expr>),
    /// Multiplication operation: left * right
    Mul(Box<Expr>, Box<Expr>),
    /// Division operation: left / right
    Div(Box<Expr>, Box<Expr>),
    /// Power operation: base ^ exponent
    Pow(Box<Expr>, Box<Expr>),
    /// Exponential function: e^x
    Exp(Box<Expr>),
    /// Natural logarithm: ln(x)
    Ln(Box<Expr>),
    /// Sine function: sin(x)
    sin(Box<Expr>),
    /// Cosine function: cos(x)
    cos(Box<Expr>),
    /// Tangent function: tan(x) - uses mathematical notation 'tg'
    tg(Box<Expr>),
    /// Cotangent function: cot(x) - uses mathematical notation 'ctg'
    ctg(Box<Expr>),
    /// Arcsine function: arcsin(x)
    arcsin(Box<Expr>),
    /// Arccosine function: arccos(x)
    arccos(Box<Expr>),
    /// Arctangent function: arctan(x) - uses mathematical notation 'arctg'
    arctg(Box<Expr>),
    /// Arccotangent function: arccot(x) - uses mathematical notation 'arcctg'
    arcctg(Box<Expr>),
}

#[derive(Clone, Debug)]
struct PrettyBlock {
    lines: Vec<String>,
    attach_row: usize,
}

impl PrettyBlock {
    fn single(line: String) -> Self {
        Self {
            lines: vec![line],
            attach_row: 0,
        }
    }

    fn width(&self) -> usize {
        self.lines
            .iter()
            .map(|line| line.chars().count())
            .max()
            .unwrap_or(0)
    }

    fn height(&self) -> usize {
        self.lines.len()
    }

    fn pad_line_right(line: &str, width: usize) -> String {
        let current = line.chars().count();
        if current >= width {
            line.to_string()
        } else {
            format!("{}{}", line, " ".repeat(width - current))
        }
    }

    fn centered_in(&self, width: usize) -> Vec<String> {
        let block_width = self.width();
        if block_width >= width {
            return self
                .lines
                .iter()
                .map(|line| Self::pad_line_right(line, block_width))
                .collect();
        }

        let pad = width - block_width;
        let left = pad / 2;
        let right = pad - left;
        self.lines
            .iter()
            .map(|line| {
                format!(
                    "{}{}{}",
                    " ".repeat(left),
                    Self::pad_line_right(line, block_width),
                    " ".repeat(right)
                )
            })
            .collect()
    }

    fn wrap_parentheses(self) -> Self {
        let Self { lines, attach_row } = self;
        if lines.is_empty() {
            return Self::single("()".to_string());
        }
        if lines.len() == 1 {
            return Self::single(format!("({})", lines[0]));
        }

        let width = lines.iter().map(|line| line.chars().count()).max().unwrap_or(0);
        let mut wrapped = Vec::with_capacity(lines.len());
        let last_index = lines.len() - 1;
        for (idx, line) in lines.into_iter().enumerate() {
            let padded = Self::pad_line_right(&line, width);
            if idx == 0 {
                wrapped.push(format!("({}", padded));
            } else if idx == last_index {
                wrapped.push(format!("{})", padded));
            } else {
                wrapped.push(format!("{} ", padded));
            }
        }
        Self {
            lines: wrapped,
            attach_row,
        }
    }
}

/// Display implementation for pretty printing symbolic expressions.
///
/// Converts expressions to human-readable mathematical notation with minimal
/// parentheses based on operator precedence rules.
///
/// ## Key Principles for Bracket Rules
///
/// ### Operator Precedence Hierarchy (highest to lowest):
/// - **Variables/Constants (100)**: Never need brackets (`x`, `5.0`)
/// - **Functions (90)**: Already have parentheses (`sin(x)`, `exp(y)`)
/// - **Power (80)**: Right-associative (`x^y^z` = `x^(y^z)`)
/// - **Multiplication/Division (60)**: Left-associative (`x*y/z` = `(x*y)/z`)
/// - **Addition/Subtraction (40)**: Left-associative (`x+y-z` = `(x+y)-z`)
///
/// ### When Brackets Are Added:
/// 1. **Lower precedence in higher precedence context**: `(x + y) * z`
/// 2. **Right-associative operations with equal/lower precedence**: `x - (y - z)`, `x / (y / z)`, `x ^ (y ^ z)`
/// 3. **Division with complex denominators**: `x / (y + z)`
/// 4. **Power with complex base**: `(x + y) ^ 2`
///
/// ### When Brackets Are Omitted:
/// 1. **Variables and constants**: Never need brackets
/// 2. **Functions**: Already have built-in parentheses
/// 3. **Same or higher precedence**: `x * y + z` (no brackets around `x * y`)
/// 4. **Left-associative with higher precedence**: `x + y * z` (no brackets around `y * z`)
///
/// ## Implementation Logic
///
/// The precedence system uses numerical values where higher numbers indicate
/// higher precedence. Two helper methods determine bracket necessity:
///
/// - `needs_brackets_left()`: Checks if left operand needs brackets (simple precedence comparison)
/// - `needs_brackets_right()`: Handles right-associativity for `-`, `/`, `^` operations
///
/// ### Examples of Output:
/// - `x + y * z` instead of `(x + (y * z))`
/// - `(x + y) * z` instead of `((x + y) * z)`
/// - `x - (y - z)` (preserves necessary brackets for right-associative subtraction)
/// - `sin(x) + y` instead of `(sin(x) + y)`
/// - `x ^ y ^ z` becomes `x ^ (y ^ z)` (right-associative power)
impl fmt::Display for Expr {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
        match self {
            Expr::Var(name) => write!(f, "{}", name),
            Expr::Const(val) => write!(f, "{}", val),

            Expr::Add(lhs, rhs) => {
                self.fmt_left_operand(f, lhs)?;
                write!(f, " + ")?;
                self.fmt_right_operand(f, rhs)
            }

            Expr::Sub(lhs, rhs) => {
                self.fmt_left_operand(f, lhs)?;
                write!(f, " - ")?;
                self.fmt_right_operand(f, rhs)
            }

            Expr::Mul(lhs, rhs) => {
                self.fmt_left_operand(f, lhs)?;
                write!(f, " * ")?;
                self.fmt_right_operand(f, rhs)
            }

            Expr::Div(lhs, rhs) => {
                self.fmt_left_operand(f, lhs)?;
                write!(f, " / ")?;
                self.fmt_right_operand(f, rhs)
            }

            Expr::Pow(base, exp) => {
                self.fmt_left_operand(f, base)?;
                write!(f, " ^ ")?;
                self.fmt_right_operand(f, exp)
            }

            Expr::Exp(expr) => write!(f, "exp({})", expr),
            Expr::Ln(expr) => write!(f, "ln({})", expr),
            Expr::sin(expr) => write!(f, "sin({})", expr),
            Expr::cos(expr) => write!(f, "cos({})", expr),
            Expr::tg(expr) => write!(f, "tg({})", expr),
            Expr::ctg(expr) => write!(f, "ctg({})", expr),
            Expr::arcsin(expr) => write!(f, "arcsin({})", expr),
            Expr::arccos(expr) => write!(f, "arccos({})", expr),
            Expr::arctg(expr) => write!(f, "arctg({})", expr),
            Expr::arcctg(expr) => write!(f, "arcctg({})", expr),
        }
    }
}

impl std::ops::Add for Expr {
    type Output = Self;

    fn add(self, rhs: Self) -> Self::Output {
        Expr::Add(self.boxed(), rhs.boxed())
    }
}

impl std::ops::Sub for Expr {
    type Output = Self;

    fn sub(self, rhs: Self) -> Self::Output {
        Expr::Sub(self.boxed(), rhs.boxed())
    }
}

impl std::ops::Mul for Expr {
    type Output = Self;

    fn mul(self, rhs: Self) -> Self::Output {
        Expr::Mul(self.boxed(), rhs.boxed())
    }
}

impl std::ops::Div for Expr {
    type Output = Self;

    fn div(self, rhs: Self) -> Self::Output {
        Expr::Div(self.boxed(), rhs.boxed())
    }
}
impl std::ops::AddAssign for Expr {
    fn add_assign(&mut self, rhs: Self) {
        *self = Expr::Add(Box::new(self.clone()), Box::new(rhs));
    }
}

impl std::ops::SubAssign for Expr {
    fn sub_assign(&mut self, rhs: Self) {
        *self = Expr::Sub(Box::new(self.clone()), Box::new(rhs));
    }
}

impl std::ops::MulAssign for Expr {
    fn mul_assign(&mut self, rhs: Self) {
        *self = Expr::Mul(Box::new(self.clone()), Box::new(rhs));
    }
}

impl std::ops::DivAssign for Expr {
    fn div_assign(&mut self, rhs: Self) {
        *self = Expr::Div(Box::new(self.clone()), Box::new(rhs));
    }
}

impl std::ops::Neg for Expr {
    type Output = Self;

    fn neg(self) -> Self::Output {
        Expr::Mul(Box::new(Expr::Const(-1.0)), Box::new(self))
    }
}

/*


impl  std::f64:: for Expr {
    type Output = Self;
    fn exp(self) -> Self::Output {
        Expr::Exp(self.boxed())
    }
}

    ipml std::f64::::pow for Expr {
    type Output = Self;

    fn pow(self, rhs: Self) -> Self::Output {
        Expr::Pow(self.boxed(), rhs.boxed())
    }
}

impl f64::ops::ln for Expr {
    type Output = Self;
    fn ln(self) -> Self::Output {
        Expr::Ln(self.boxed())
    }
}
    */
// Implement differentiation, based on the recursive definition

impl Expr {
    /// Get operator precedence (higher number = higher precedence)
    fn precedence(&self) -> u8 {
        match self {
            Expr::Var(_) | Expr::Const(_) => 100, // Highest - never need brackets
            Expr::Exp(_)
            | Expr::Ln(_)
            | Expr::sin(_)
            | Expr::cos(_)
            | Expr::tg(_)
            | Expr::ctg(_)
            | Expr::arcsin(_)
            | Expr::arccos(_)
            | Expr::arctg(_)
            | Expr::arcctg(_) => 90, // Functions
            Expr::Pow(_, _) => 80,                // Power
            Expr::Mul(_, _) | Expr::Div(_, _) => 60, // Multiplication/Division
            Expr::Add(_, _) | Expr::Sub(_, _) => 40, // Addition/Subtraction
        }
    }

    /// Check if operand needs brackets when used as left operand
    fn needs_brackets_left(&self, operand: &Expr) -> bool {
        operand.precedence() < self.precedence()
    }

    /// Check if operand needs brackets when used as right operand
    fn needs_brackets_right(&self, operand: &Expr) -> bool {
        match self {
            Expr::Sub(_, _) | Expr::Div(_, _) | Expr::Pow(_, _) => {
                operand.precedence() <= self.precedence()
            }
            _ => operand.precedence() < self.precedence(),
        }
    }

    /// Format left operand with conditional brackets
    fn fmt_left_operand(&self, f: &mut fmt::Formatter, operand: &Expr) -> fmt::Result {
        if self.needs_brackets_left(operand) {
            write!(f, "({})", operand)
        } else {
            write!(f, "{}", operand)
        }
    }

    /// Format right operand with conditional brackets
    fn fmt_right_operand(&self, f: &mut fmt::Formatter, operand: &Expr) -> fmt::Result {
        if self.needs_brackets_right(operand) {
            write!(f, "({})", operand)
        } else {
            write!(f, "{}", operand)
        }
    }

    /// Render an expression in a fully explicit inline form.
    ///
    /// This is used only as a fallback for nested fractions where the main
    /// layout engine would otherwise lose grouping information.
    fn to_explicit_inline_string(&self) -> String {
        match self {
            Expr::Var(_) | Expr::Const(_) => self.to_string(),
            Expr::Add(lhs, rhs) => format!(
                "({} + {})",
                lhs.to_explicit_inline_string(),
                rhs.to_explicit_inline_string()
            ),
            Expr::Sub(lhs, rhs) => format!(
                "({} - {})",
                lhs.to_explicit_inline_string(),
                rhs.to_explicit_inline_string()
            ),
            Expr::Mul(lhs, rhs) => format!(
                "({} * {})",
                lhs.to_explicit_inline_string(),
                rhs.to_explicit_inline_string()
            ),
            Expr::Div(lhs, rhs) => format!(
                "({} / {})",
                lhs.to_explicit_inline_string(),
                rhs.to_explicit_inline_string()
            ),
            Expr::Pow(base, exp) => format!(
                "({} ^ {})",
                base.to_explicit_inline_string(),
                exp.to_explicit_inline_string()
            ),
            Expr::Exp(expr) => format!("exp({})", expr.to_explicit_inline_string()),
            Expr::Ln(expr) => format!("ln({})", expr.to_explicit_inline_string()),
            Expr::sin(expr) => format!("sin({})", expr.to_explicit_inline_string()),
            Expr::cos(expr) => format!("cos({})", expr.to_explicit_inline_string()),
            Expr::tg(expr) => format!("tg({})", expr.to_explicit_inline_string()),
            Expr::ctg(expr) => format!("ctg({})", expr.to_explicit_inline_string()),
            Expr::arcsin(expr) => format!("arcsin({})", expr.to_explicit_inline_string()),
            Expr::arccos(expr) => format!("arccos({})", expr.to_explicit_inline_string()),
            Expr::arctg(expr) => format!("arctg({})", expr.to_explicit_inline_string()),
            Expr::arcctg(expr) => format!("arcctg({})", expr.to_explicit_inline_string()),
        }
    }

    /// Format function with multi-line content horizontally
    fn _format_function_horizontal(func_name: &str, inner_lines: &[String]) -> Vec<String> {
        if inner_lines.len() <= 1 {
            return vec![format!(
                "{}({})",
                func_name,
                inner_lines.get(0).unwrap_or(&String::new())
            )];
        }

        let mut result = Vec::new();

        // First line: function_name( first_line )
        result.push(format!("{}( {} )", func_name, inner_lines[0]));

        // Middle lines: align content properly with the content position, not the parenthesis
        let content_offset = func_name.len() + 2; // "func( " length
        for i in 1..inner_lines.len() {
            result.push(format!("{}{}", " ".repeat(content_offset), inner_lines[i]));
        }

        result
    }

    /// Convert digit to Unicode superscript
    fn to_superscript_digit(digit: u8) -> char {
        match digit {
            0 => '⁰',
            1 => '¹',
            2 => '²',
            3 => '³',
            4 => '⁴',
            5 => '⁵',
            6 => '⁶',
            7 => '⁷',
            8 => '⁸',
            9 => '⁹',
            _ => '?',
        }
    }

    /// Convert single character to Unicode superscript if possible
    fn to_superscript_char(c: char) -> Option<char> {
        match c {
            'a' => Some('ᵃ'),
            'b' => Some('ᵇ'),
            'c' => Some('ᶜ'),
            'd' => Some('ᵈ'),
            'e' => Some('ᵉ'),
            'f' => Some('ᶠ'),
            'g' => Some('ᵍ'),
            'h' => Some('ʰ'),
            'i' => Some('ⁱ'),
            'j' => Some('ʲ'),
            'k' => Some('ᵏ'),
            'l' => Some('ˡ'),
            'm' => Some('ᵐ'),
            'n' => Some('ⁿ'),
            'o' => Some('ᵒ'),
            'p' => Some('ᵖ'),
            'r' => Some('ʳ'),
            's' => Some('ˢ'),
            't' => Some('ᵗ'),
            'u' => Some('ᵘ'),
            'v' => Some('ᵛ'),
            'w' => Some('ʷ'),
            'x' => Some('ˣ'),
            'y' => Some('ʸ'),
            'z' => Some('ᶻ'),
            '+' => Some('⁺'),
            '-' => Some('⁻'),
            '=' => Some('⁼'),
            '(' => Some('⁽'),
            ')' => Some('⁾'),
            _ => None,
        }
    }

    /// Check if expression can be rendered as simple Unicode superscript
    fn can_use_unicode_superscript(&self) -> bool {
        match self {
            Expr::Const(n) => *n >= 0.0 && *n <= 9.0 && n.fract() == 0.0,
            Expr::Var(name) => {
                name.len() == 1 && Self::to_superscript_char(name.chars().next().unwrap()).is_some()
            }
            _ => false,
        }
    }

    /// Convert expression to Unicode superscript string
    fn to_unicode_superscript(&self) -> String {
        match self {
            Expr::Const(n) => {
                let digit = *n as u8;
                Self::to_superscript_digit(digit).to_string()
            }
            Expr::Var(name) => {
                let c = name.chars().next().unwrap();
                Self::to_superscript_char(c).unwrap().to_string()
            }
            _ => format!("^({})", self),
        }
    }

    /// Convert Greek letter names to Unicode symbols
    fn convert_greek_letters(&self) -> Expr {
        match self {
            Expr::Var(name) => {
                let lower_name = name.to_lowercase();

                // Check for Greek letters with or without indices
                for (greek_name, greek_symbol) in [
                    ("alpha", "α"),
                    ("beta", "β"),
                    ("gamma", "γ"),
                    ("delta", "δ"),
                    ("epsilon", "ε"),
                    ("zeta", "ζ"),
                    ("eta", "η"),
                    ("theta", "θ"),
                    ("iota", "ι"),
                    ("kappa", "κ"),
                    ("lambda", "λ"),
                    ("mu", "μ"),
                    ("nu", "ν"),
                    ("xi", "ξ"),
                    ("omicron", "ο"),
                    ("pi", "π"),
                    ("rho", "ρ"),
                    ("sigma", "σ"),
                    ("tau", "τ"),
                    ("upsilon", "υ"),
                    ("phi", "φ"),
                    ("chi", "χ"),
                    ("psi", "ψ"),
                    ("omega", "ω"),
                ] {
                    if lower_name.starts_with(greek_name) {
                        let suffix = &lower_name[greek_name.len()..];
                        return Expr::Var(format!("{}{}", greek_symbol, suffix));
                    }
                }

                self.clone()
            }
            Expr::Add(lhs, rhs) => Expr::Add(
                Box::new(lhs.convert_greek_letters()),
                Box::new(rhs.convert_greek_letters()),
            ),
            Expr::Sub(lhs, rhs) => Expr::Sub(
                Box::new(lhs.convert_greek_letters()),
                Box::new(rhs.convert_greek_letters()),
            ),
            Expr::Mul(lhs, rhs) => Expr::Mul(
                Box::new(lhs.convert_greek_letters()),
                Box::new(rhs.convert_greek_letters()),
            ),
            Expr::Div(lhs, rhs) => Expr::Div(
                Box::new(lhs.convert_greek_letters()),
                Box::new(rhs.convert_greek_letters()),
            ),
            Expr::Pow(base, exp) => Expr::Pow(
                Box::new(base.convert_greek_letters()),
                Box::new(exp.convert_greek_letters()),
            ),
            Expr::Exp(expr) => Expr::Exp(Box::new(expr.convert_greek_letters())),
            Expr::Ln(expr) => Expr::Ln(Box::new(expr.convert_greek_letters())),
            Expr::sin(expr) => Expr::sin(Box::new(expr.convert_greek_letters())),
            Expr::cos(expr) => Expr::cos(Box::new(expr.convert_greek_letters())),
            Expr::tg(expr) => Expr::tg(Box::new(expr.convert_greek_letters())),
            Expr::ctg(expr) => Expr::ctg(Box::new(expr.convert_greek_letters())),
            Expr::arcsin(expr) => Expr::arcsin(Box::new(expr.convert_greek_letters())),
            Expr::arccos(expr) => Expr::arccos(Box::new(expr.convert_greek_letters())),
            Expr::arctg(expr) => Expr::arctg(Box::new(expr.convert_greek_letters())),
            Expr::arcctg(expr) => Expr::arcctg(Box::new(expr.convert_greek_letters())),
            _ => self.clone(),
        }
    }

    /// Pretty print with mathematical formatting
    pub fn pretty_print(&self) -> String {
        self.convert_greek_letters()
            .pretty_print_internal()
            .join("\n")
    }

    /// Calculate the actual visual width of an expression when pretty printed
    #[allow(dead_code)]
    fn visual_width(&self) -> usize {
        match self {
            Expr::Var(name) => name.chars().count(),
            Expr::Const(val) => format!("{}", val).chars().count(),
            Expr::Div(num, den) => {
                let num_width = num.visual_width();
                let den_width = den.visual_width();
                num_width.max(den_width).max(3)
            }
            Expr::Pow(base, exp) => {
                if exp.can_use_unicode_superscript() {
                    base.visual_width() + exp.to_unicode_superscript().chars().count()
                } else {
                    base.visual_width().max(exp.visual_width())
                }
            }
            _ => {
                // For other operations, calculate based on operands
                match self {
                    Expr::Add(lhs, rhs) | Expr::Sub(lhs, rhs) | Expr::Mul(lhs, rhs) => {
                        lhs.visual_width() + 3 + rhs.visual_width() // operator + spaces
                    }
                    Expr::Exp(_)
                    | Expr::Ln(_)
                    | Expr::sin(_)
                    | Expr::cos(_)
                    | Expr::tg(_)
                    | Expr::ctg(_)
                    | Expr::arcsin(_)
                    | Expr::arccos(_)
                    | Expr::arctg(_)
                    | Expr::arcctg(_) => {
                        // Function name + parentheses + content
                        let func_name = match self {
                            Expr::Exp(_) => "exp",
                            Expr::Ln(_) => "ln",
                            Expr::sin(_) => "sin",
                            Expr::cos(_) => "cos",
                            Expr::tg(_) => "tg",
                            Expr::ctg(_) => "ctg",
                            Expr::arcsin(_) => "arcsin",
                            Expr::arccos(_) => "arccos",
                            Expr::arctg(_) => "arctg",
                            Expr::arcctg(_) => "arcctg",
                            _ => "",
                        };
                        let inner = match self {
                            Expr::Exp(e)
                            | Expr::Ln(e)
                            | Expr::sin(e)
                            | Expr::cos(e)
                            | Expr::tg(e)
                            | Expr::ctg(e)
                            | Expr::arcsin(e)
                            | Expr::arccos(e)
                            | Expr::arctg(e)
                            | Expr::arcctg(e) => e.visual_width(),
                            _ => 0,
                        };
                        func_name.chars().count() + 2 + inner // name + "()" + content
                    }
                    _ => 10, // fallback
                }
            }
        }
    }

    /// Internal pretty print using level-based layout system
    fn pretty_print_internal(&self) -> Vec<String> {
        self.pretty_block().lines
    }

    fn pretty_block(&self) -> PrettyBlock {
        match self {
            Expr::Var(name) => PrettyBlock::single(name.clone()),
            Expr::Const(val) => PrettyBlock::single(format!("{}", val)),
            Expr::Add(lhs, rhs) => Self::combine_infix_blocks(
                lhs.pretty_block(),
                " + ",
                rhs.pretty_block(),
                self.needs_brackets_left(lhs),
                self.needs_brackets_right(rhs),
            ),
            Expr::Sub(lhs, rhs) => Self::combine_infix_blocks(
                lhs.pretty_block(),
                " - ",
                rhs.pretty_block(),
                self.needs_brackets_left(lhs),
                self.needs_brackets_right(rhs),
            ),
            Expr::Mul(lhs, rhs) => Self::combine_infix_blocks(
                lhs.pretty_block(),
                " * ",
                rhs.pretty_block(),
                self.needs_brackets_left(lhs),
                self.needs_brackets_right(rhs),
            ),
            Expr::Div(num, den) => Self::combine_fraction_blocks(num.pretty_block(), den.pretty_block()),
            Expr::Pow(base, exp) => {
                let exp_block = exp.pretty_block();
                let mut base_block = base.pretty_block();

                if exp.can_use_unicode_superscript()
                    && exp_block.lines.len() == 1
                    && base_block.lines.len() == 1
                {
                    let base_str = if self.needs_brackets_left(base) {
                        format!("({})", base_block.lines[0].clone())
                    } else {
                        base_block.lines[0].clone()
                    };
                    PrettyBlock::single(format!("{}{}", base_str, exp.to_unicode_superscript()))
                } else {
                    let base_is_fraction = matches!(**base, Expr::Div(_, _));
                    if self.needs_brackets_left(base) && !base_is_fraction {
                        base_block = base_block.wrap_parentheses();
                    }
                    Self::render_power_block(base_block, exp_block)
                }
            }
            Expr::Exp(expr) => {
                let exp_block = expr.pretty_block();
                if expr.can_use_unicode_superscript() && exp_block.lines.len() == 1 {
                    PrettyBlock::single(format!("e{}", expr.to_unicode_superscript()))
                } else {
                    Self::render_exp_block(exp_block)
                }
            }
            Expr::Ln(expr) => Self::render_function_block("ln", expr.pretty_block()),
            Expr::sin(expr) => Self::render_function_block("sin", expr.pretty_block()),
            Expr::cos(expr) => Self::render_function_block("cos", expr.pretty_block()),
            Expr::tg(expr) => Self::render_function_block("tg", expr.pretty_block()),
            Expr::ctg(expr) => Self::render_function_block("ctg", expr.pretty_block()),
            Expr::arcsin(expr) => Self::render_function_block("arcsin", expr.pretty_block()),
            Expr::arccos(expr) => Self::render_function_block("arccos", expr.pretty_block()),
            Expr::arctg(expr) => Self::render_function_block("arctg", expr.pretty_block()),
            Expr::arcctg(expr) => Self::render_function_block("arcctg", expr.pretty_block()),
        }
    }

    fn combine_infix_blocks(
        mut left: PrettyBlock,
        operator: &str,
        mut right: PrettyBlock,
        left_needs_brackets: bool,
        right_needs_brackets: bool,
    ) -> PrettyBlock {
        if left_needs_brackets {
            left = left.wrap_parentheses();
        }
        if right_needs_brackets {
            right = right.wrap_parentheses();
        }

        if left.lines.len() == 1 && right.lines.len() == 1 {
            return PrettyBlock::single(format!(
                "{}{}{}",
                left.lines[0], operator, right.lines[0]
            ));
        }

        let left_width = left.width();
        let right_width = right.width();
        let op_width = operator.chars().count();
        let attach_row = left.attach_row.max(right.attach_row);
        let left_top_pad = attach_row.saturating_sub(left.attach_row);
        let right_top_pad = attach_row.saturating_sub(right.attach_row);
        let total_height = (left_top_pad + left.height()).max(right_top_pad + right.height());

        let mut lines = Vec::with_capacity(total_height);
        for row in 0..total_height {
            let left_line = if row >= left_top_pad && row < left_top_pad + left.height() {
                PrettyBlock::pad_line_right(&left.lines[row - left_top_pad], left_width)
            } else {
                " ".repeat(left_width)
            };

            let right_line = if row >= right_top_pad && row < right_top_pad + right.height() {
                PrettyBlock::pad_line_right(&right.lines[row - right_top_pad], right_width)
            } else {
                " ".repeat(right_width)
            };

            let middle = if row == attach_row {
                operator.to_string()
            } else {
                " ".repeat(op_width)
            };
            lines.push(format!("{}{}{}", left_line, middle, right_line));
        }

        PrettyBlock { lines, attach_row }
    }

    fn combine_fraction_blocks(num: PrettyBlock, den: PrettyBlock) -> PrettyBlock {
        let width = num.width().max(den.width()).max(3);
        let mut lines = Vec::new();
        lines.extend(num.centered_in(width));
        lines.push("─".repeat(width));
        lines.extend(den.centered_in(width));
        PrettyBlock {
            lines,
            attach_row: num.height(),
        }
    }

    fn render_block_at(
        canvas: &mut [Vec<char>],
        block: &PrettyBlock,
        top: usize,
        left: usize,
    ) {
        for (row_idx, line) in block.lines.iter().enumerate() {
            let canvas_row = top + row_idx;
            if canvas_row >= canvas.len() {
                continue;
            }
            for (col_idx, ch) in line.chars().enumerate() {
                let canvas_col = left + col_idx;
                if canvas_col < canvas[canvas_row].len() && ch != ' ' {
                    canvas[canvas_row][canvas_col] = ch;
                }
            }
        }
    }

    fn render_power_block(base: PrettyBlock, exp: PrettyBlock) -> PrettyBlock {
        let base_width = base.width();
        let exp_width = exp.width();
        let exp_left = base_width + 1;
        let total_width = exp_left + exp_width;
        let base_is_exp_marker = base.height() == 1
            && base
                .lines
                .first()
                .map(|line| line.as_str())
                == Some("e");
        let base_top = if exp.height() == 1 {
            1
        } else if base.height() > 1 || base_is_exp_marker {
            exp.height() + 1
        } else {
            exp.height() - 1
        };
        let total_height = (base_top + base.height()).max(exp.height());
        let attach_row = base_top + base.attach_row;

        let mut canvas = vec![vec![' '; total_width]; total_height];
        Self::render_block_at(&mut canvas, &exp, 0, exp_left);
        Self::render_block_at(&mut canvas, &base, base_top, 0);

        let lines = canvas
            .into_iter()
            .map(|row| row.into_iter().collect::<String>().trim_end().to_string())
            .collect::<Vec<_>>();

        PrettyBlock { lines, attach_row }
    }

    fn render_exp_block(exp: PrettyBlock) -> PrettyBlock {
        // Natural exponentials should follow the same power geometry as other
        // powers, but with `e` as the base instead of a textual `e^(...)`.
        Self::render_power_block(PrettyBlock::single("e".to_string()), exp)
    }

    fn render_function_block(func_name: &str, inner: PrettyBlock) -> PrettyBlock {
        if inner.lines.len() == 1 {
            return PrettyBlock::single(format!("{}({})", func_name, inner.lines[0]));
        }

        let name_width = func_name.chars().count();
        let prefix_width = name_width + 1;
        let inner_width = inner.width();
        let total_width = prefix_width + inner_width + 2;
        let total_height = inner.height();
        let attach_row = inner.attach_row.min(total_height.saturating_sub(1));

        // A multi-line argument is wrapped in scalable mathematical
        // parentheses. The function name sits on the argument's attachment
        // axis, so an outer infix operator composes with the complete call.
        let mut canvas = vec![vec![' '; total_width]; total_height];
        for (idx, ch) in func_name.chars().enumerate() {
            canvas[attach_row][idx] = ch;
        }

        Self::render_block_at(&mut canvas, &inner, 0, prefix_width + 1);
        let right_col = prefix_width + inner_width + 1;
        for row in 0..total_height {
            let (left_bracket, right_bracket) = if row == 0 {
                ('⎛', '⎞')
            } else if row + 1 == total_height {
                ('⎝', '⎠')
            } else {
                ('⎜', '⎟')
            };
            canvas[row][prefix_width] = left_bracket;
            canvas[row][right_col] = right_bracket;
        }

        let lines = canvas
            .into_iter()
            .map(|row| row.into_iter().collect::<String>().trim_end().to_string())
            .collect::<Vec<_>>();

        PrettyBlock { lines, attach_row }
    }

    /// BASIC FEATURES

    /// Check if operand needs brackets (continued)
    fn _needs_brackets_continued(&self, operand: &Expr) -> bool {
        operand.precedence() < self.precedence()
    }

    /// BASIC FEATURES

    /// Creates multiple symbolic variables from a comma-separated string.
    ///
    /// Parses a string containing variable names separated by commas and returns
    /// a vector of Expr::Var instances. Whitespace is automatically trimmed.
    ///
    /// # Arguments
    /// * `symbols` - Comma-separated string of variable names (e.g., "x, y, z")
    ///
    /// # Returns
    /// Vector of Expr::Var instances for each variable name
    ///
    /// # Examples
    /// ```rust, ignore
    /// let vars = Expr::Symbols("x, y, z");
    /// assert_eq!(vars.len(), 3);
    /// ```
    pub fn Symbols(symbols: &str) -> Vec<Expr> {
        let symbols = symbols.to_string();
        let vec_trimmed: Vec<String> = symbols.split(',').map(|s| s.trim().to_string()).collect();
        let vector_of_symbolic_vars: Vec<Expr> = vec_trimmed
            .iter()
            .filter(|s| !s.is_empty())
            .map(|s| Expr::Var(s.to_string()))
            .collect();
        vector_of_symbolic_vars
    }
    /// Substitutes a variable with a constant value throughout the expression.
    ///
    /// Recursively traverses the expression tree and replaces all occurrences
    /// of the specified variable with the given constant value.
    ///
    /// # Arguments
    /// * `var` - Name of the variable to substitute
    /// * `value` - Numerical value to substitute for the variable
    ///
    /// # Returns
    /// New expression with the variable substituted
    pub fn set_variable(&self, var: &str, value: f64) -> Expr {
        match self {
            Expr::Var(name) if name == var => Expr::Const(value),
            Expr::Add(lhs, rhs) => Expr::Add(
                Box::new(lhs.set_variable(var, value)),
                Box::new(rhs.set_variable(var, value)),
            ),
            Expr::Sub(lhs, rhs) => Expr::Sub(
                Box::new(lhs.set_variable(var, value)),
                Box::new(rhs.set_variable(var, value)),
            ),
            Expr::Mul(lhs, rhs) => Expr::Mul(
                Box::new(lhs.set_variable(var, value)),
                Box::new(rhs.set_variable(var, value)),
            ),
            Expr::Div(lhs, rhs) => Expr::Div(
                Box::new(lhs.set_variable(var, value)),
                Box::new(rhs.set_variable(var, value)),
            ),
            Expr::Pow(base, exp) => Expr::Pow(
                Box::new(base.set_variable(var, value)),
                Box::new(exp.set_variable(var, value)),
            ),
            Expr::Exp(expr) => Expr::Exp(Box::new(expr.set_variable(var, value))),
            Expr::Ln(expr) => Expr::Ln(Box::new(expr.set_variable(var, value))),
            Expr::sin(expr) => Expr::sin(Box::new(expr.set_variable(var, value))),
            Expr::cos(expr) => Expr::cos(Box::new(expr.set_variable(var, value))),
            Expr::tg(expr) => Expr::tg(Box::new(expr.set_variable(var, value))),
            Expr::ctg(expr) => Expr::ctg(Box::new(expr.set_variable(var, value))),
            Expr::arcsin(expr) => Expr::arcsin(Box::new(expr.set_variable(var, value))),
            Expr::arccos(expr) => Expr::arccos(Box::new(expr.set_variable(var, value))),
            Expr::arctg(expr) => Expr::arctg(Box::new(expr.set_variable(var, value))),
            Expr::arcctg(expr) => Expr::arcctg(Box::new(expr.set_variable(var, value))),
            _ => self.clone(),
        }
    }

    /// Substitutes multiple variables with constant values using a HashMap.
    ///
    /// More efficient than multiple set_variable calls when substituting many variables.
    /// Only variables present in the map are substituted.
    ///
    /// # Arguments
    /// * `var_map` - HashMap mapping variable names to their replacement values
    ///
    /// # Returns
    /// New expression with all mapped variables substituted
    pub fn set_variable_from_map(&self, var_map: &HashMap<String, f64>) -> Expr {
        match self {
            Expr::Var(name) if var_map.contains_key(name) => Expr::Const(var_map[name]),
            Expr::Add(lhs, rhs) => Expr::Add(
                Box::new(lhs.set_variable_from_map(var_map)),
                Box::new(rhs.set_variable_from_map(var_map)),
            ),
            Expr::Sub(lhs, rhs) => Expr::Sub(
                Box::new(lhs.set_variable_from_map(var_map)),
                Box::new(rhs.set_variable_from_map(var_map)),
            ),
            Expr::Mul(lhs, rhs) => Expr::Mul(
                Box::new(lhs.set_variable_from_map(var_map)),
                Box::new(rhs.set_variable_from_map(var_map)),
            ),
            Expr::Div(lhs, rhs) => Expr::Div(
                Box::new(lhs.set_variable_from_map(var_map)),
                Box::new(rhs.set_variable_from_map(var_map)),
            ),
            Expr::Pow(base, exp) => Expr::Pow(
                Box::new(base.set_variable_from_map(var_map)),
                Box::new(exp.set_variable_from_map(var_map)),
            ),
            Expr::Exp(expr) => Expr::Exp(Box::new(expr.set_variable_from_map(var_map))),
            Expr::Ln(expr) => Expr::Ln(Box::new(expr.set_variable_from_map(var_map))),

            Expr::sin(expr) => Expr::sin(Box::new(expr.set_variable_from_map(var_map))),
            Expr::cos(expr) => Expr::cos(Box::new(expr.set_variable_from_map(var_map))),
            Expr::tg(expr) => Expr::tg(Box::new(expr.set_variable_from_map(var_map))),
            Expr::ctg(expr) => Expr::ctg(Box::new(expr.set_variable_from_map(var_map))),
            Expr::arcsin(expr) => Expr::arcsin(Box::new(expr.set_variable_from_map(var_map))),
            Expr::arccos(expr) => Expr::arccos(Box::new(expr.set_variable_from_map(var_map))),
            Expr::arctg(expr) => Expr::arctg(Box::new(expr.set_variable_from_map(var_map))),
            Expr::arcctg(expr) => Expr::arcctg(Box::new(expr.set_variable_from_map(var_map))),
            _ => self.clone(),
        }
    }
    /// Renames a variable throughout the expression.
    ///
    /// Recursively replaces all occurrences of old_var with new_var.
    /// Useful for variable substitution and expression manipulation.
    ///
    /// # Arguments
    /// * `old_var` - Current variable name to replace
    /// * `new_var` - New variable name
    ///
    /// # Returns
    /// New expression with the variable renamed
    pub fn rename_variable(&self, old_var: &str, new_var: &str) -> Expr {
        match self {
            Expr::Var(name) if name == old_var => Expr::Var(new_var.to_string()),
            Expr::Add(lhs, rhs) => Expr::Add(
                Box::new(lhs.rename_variable(old_var, new_var)),
                Box::new(rhs.rename_variable(old_var, new_var)),
            ),
            Expr::Sub(lhs, rhs) => Expr::Sub(
                Box::new(lhs.rename_variable(old_var, new_var)),
                Box::new(rhs.rename_variable(old_var, new_var)),
            ),
            Expr::Mul(lhs, rhs) => Expr::Mul(
                Box::new(lhs.rename_variable(old_var, new_var)),
                Box::new(rhs.rename_variable(old_var, new_var)),
            ),
            Expr::Div(lhs, rhs) => Expr::Div(
                Box::new(lhs.rename_variable(old_var, new_var)),
                Box::new(rhs.rename_variable(old_var, new_var)),
            ),
            Expr::Pow(base, exp) => Expr::Pow(
                Box::new(base.rename_variable(old_var, new_var)),
                Box::new(exp.rename_variable(old_var, new_var)),
            ),
            Expr::Exp(expr) => Expr::Exp(Box::new(expr.rename_variable(old_var, new_var))),
            Expr::Ln(expr) => Expr::Ln(Box::new(expr.rename_variable(old_var, new_var))),

            Expr::sin(expr) => Expr::sin(Box::new(expr.rename_variable(old_var, new_var))),
            Expr::cos(expr) => Expr::cos(Box::new(expr.rename_variable(old_var, new_var))),
            Expr::tg(expr) => Expr::tg(Box::new(expr.rename_variable(old_var, new_var))),
            Expr::ctg(expr) => Expr::ctg(Box::new(expr.rename_variable(old_var, new_var))),
            Expr::arcsin(expr) => Expr::arcsin(Box::new(expr.rename_variable(old_var, new_var))),
            Expr::arccos(expr) => Expr::arccos(Box::new(expr.rename_variable(old_var, new_var))),
            Expr::arctg(expr) => Expr::arctg(Box::new(expr.rename_variable(old_var, new_var))),
            Expr::arcctg(expr) => Expr::arcctg(Box::new(expr.rename_variable(old_var, new_var))),

            _ => self.clone(),
        }
    }
    /// rename variables from a map
    pub fn rename_variables(&self, var_map: &HashMap<String, String>) -> Expr {
        match self {
            Expr::Var(name) if var_map.contains_key(name) => Expr::Var(var_map[name].to_string()),
            Expr::Add(lhs, rhs) => Expr::Add(
                Box::new(lhs.rename_variables(var_map)),
                Box::new(rhs.rename_variables(var_map)),
            ),
            Expr::Sub(lhs, rhs) => Expr::Sub(
                Box::new(lhs.rename_variables(var_map)),
                Box::new(rhs.rename_variables(var_map)),
            ),
            Expr::Mul(lhs, rhs) => Expr::Mul(
                Box::new(lhs.rename_variables(var_map)),
                Box::new(rhs.rename_variables(var_map)),
            ),
            Expr::Div(lhs, rhs) => Expr::Div(
                Box::new(lhs.rename_variables(var_map)),
                Box::new(rhs.rename_variables(var_map)),
            ),
            Expr::Pow(base, exp) => Expr::Pow(
                Box::new(base.rename_variables(var_map)),
                Box::new(exp.rename_variables(var_map)),
            ),
            Expr::Exp(expr) => Expr::Exp(Box::new(expr.rename_variables(var_map))),
            Expr::Ln(expr) => Expr::Ln(Box::new(expr.rename_variables(var_map))),

            Expr::sin(expr) => Expr::sin(Box::new(expr.rename_variables(var_map))),
            Expr::cos(expr) => Expr::cos(Box::new(expr.rename_variables(var_map))),
            Expr::tg(expr) => Expr::tg(Box::new(expr.rename_variables(var_map))),
            Expr::ctg(expr) => Expr::ctg(Box::new(expr.rename_variables(var_map))),
            Expr::arcsin(expr) => Expr::arcsin(Box::new(expr.rename_variables(var_map))),
            Expr::arccos(expr) => Expr::arccos(Box::new(expr.rename_variables(var_map))),
            Expr::arctg(expr) => Expr::arctg(Box::new(expr.rename_variables(var_map))),

            _ => self.clone(),
        }
    }
    /// substitute a variable with an expression
    pub fn substitute_variable(&self, var: &str, expr: &Expr) -> Expr {
        match self {
            Expr::Var(name) if name == var => expr.clone(),
            Expr::Add(lhs, rhs) => Expr::Add(
                Box::new(lhs.substitute_variable(var, expr)),
                Box::new(rhs.substitute_variable(var, expr)),
            ),
            Expr::Sub(lhs, rhs) => Expr::Sub(
                Box::new(lhs.substitute_variable(var, expr)),
                Box::new(rhs.substitute_variable(var, expr)),
            ),
            Expr::Mul(lhs, rhs) => Expr::Mul(
                Box::new(lhs.substitute_variable(var, expr)),
                Box::new(rhs.substitute_variable(var, expr)),
            ),
            Expr::Div(lhs, rhs) => Expr::Div(
                Box::new(lhs.substitute_variable(var, expr)),
                Box::new(rhs.substitute_variable(var, expr)),
            ),
            Expr::Pow(base, exp) => Expr::Pow(
                Box::new(base.substitute_variable(var, expr)),
                Box::new(exp.substitute_variable(var, expr)),
            ),
            Expr::Exp(expr_old) => Expr::Exp(Box::new(expr_old.substitute_variable(var, expr))),
            Expr::Ln(expr_old) => Expr::Ln(Box::new(expr_old.substitute_variable(var, expr))),

            Expr::sin(expr) => Expr::sin(Box::new(expr.substitute_variable(var, expr))),
            Expr::cos(expr) => Expr::cos(Box::new(expr.substitute_variable(var, expr))),
            Expr::tg(expr) => Expr::tg(Box::new(expr.substitute_variable(var, expr))),
            Expr::ctg(expr) => Expr::ctg(Box::new(expr.substitute_variable(var, expr))),
            Expr::arcsin(expr) => Expr::arcsin(Box::new(expr.substitute_variable(var, expr))),
            Expr::arccos(expr) => Expr::arccos(Box::new(expr.substitute_variable(var, expr))),
            Expr::arctg(expr) => Expr::arctg(Box::new(expr.substitute_variable(var, expr))),
            Expr::arcctg(expr) => Expr::arcctg(Box::new(expr.substitute_variable(var, expr))),
            _ => self.clone(),
        }
    }
    /// check if the expression contains a variable
    pub fn contains_variable(&self, var_name: &str) -> bool {
        let mut cache = TraversalCache::new();
        cache.contains_variable(self, var_name)
    }

    /// Convenience method to wrap expression in Box for recursive structures.
    ///
    /// Essential for creating nested expressions since Expr variants use Box<Expr>.
    pub fn boxed(self) -> Box<Self> {
        Box::new(self)
    }

    /// Mutates current expression to a variable and returns the new variable expression.
    ///
    /// # Arguments
    /// * `var` - Variable name to create
    pub fn var_expr(&mut self, var: &str) -> Expr {
        let expr = Expr::Var(var.to_string());
        *self = expr.clone();
        expr
    }

    /// Mutates current expression to a constant value.
    ///
    /// # Arguments
    /// * `val` - Constant value to set
    pub fn const_expr(&mut self, val: f64) {
        *self = Expr::Const(val);
    }

    /// Creates exponential function e^(self).
    ///
    /// # Returns
    /// New Expr::Exp containing this expression
    pub fn exp(mut self) -> Expr {
        self = Expr::Exp(self.boxed());
        self
    }

    /// Creates natural logarithm ln(self).
    ///
    /// # Returns
    /// New Expr::Ln containing this expression
    pub fn ln(mut self) -> Expr {
        self = Expr::Ln(self.boxed());
        self
    }

    /// Creates base-10 logarithm using ln(x)/ln(10) identity.
    ///
    /// # Returns
    /// Expression equivalent to log10(self)
    pub fn log10(mut self) -> Expr {
        self = Expr::Ln(self.boxed()) / Expr::Const(2.30258509);
        self
    }

    /// Creates power expression self^rhs.
    ///
    /// # Arguments
    /// * `rhs` - Exponent expression
    ///
    /// # Returns
    /// New Expr::Pow with self as base and rhs as exponent
    pub fn pow(mut self, rhs: Expr) -> Expr {
        self = Expr::Pow(self.boxed(), rhs.boxed());
        self
    }

    /// Checks if expression is exactly zero (constant 0.0).
    ///
    /// # Returns
    /// true if expression is Const(0.0), false otherwise
    pub fn is_zero(&self) -> bool {
        match self {
            Expr::Const(val) => val == &0.0,
            _ => false,
        }
    }

    //__________________________________INDEXED VARIABLES____________________________________

    /// Creates a single indexed variable with format "name + index" (e.g., "x5").
    ///
    /// Useful for creating sequences of related variables in mathematical systems.
    ///
    /// # Arguments
    /// * `index` - Numerical index to append
    /// * `var_name` - Base variable name
    ///
    /// # Returns
    /// Expr::Var with indexed name
    pub fn IndexedVar(index: usize, var_name: &str) -> Expr {
        let indexed_var_name = format!("{}{}", var_name, index);
        Expr::Var(indexed_var_name)
    }

    /// Creates multiple indexed variables and their string representations.
    ///
    /// Generates a sequence of variables like x0, x1, x2, ... up to num_vars-1.
    /// Returns both the Expr objects and their string names for convenience.
    ///
    /// # Arguments
    /// * `num_vars` - Number of variables to create
    /// * `var_name` - Base name for variables
    ///
    /// # Returns
    /// Tuple of (Vec<Expr>, Vec<String>) containing expressions and names
    pub fn IndexedVars(num_vars: usize, var_name: &str) -> (Vec<Expr>, Vec<String>) {
        let vec_of_expr = (0..num_vars)
            .map(|i| Expr::IndexedVar(i, var_name))
            .collect();
        let vec_of_names = (0..num_vars)
            .map(|i| format!("{}_{}", var_name, i))
            .collect();
        (vec_of_expr, vec_of_names)
    }
    pub fn IndexedVarsMatrix(
        num_vars: usize,
        var_names: Vec<String>,
    ) -> (Vec<Vec<Expr>>, Vec<Vec<String>>) {
        let mut matrix = Vec::new();
        let mut matrix_of_expr = Vec::new();
        for i in 0..num_vars {
            let mut matrix_i = Vec::new();
            let mut matrix_of_expr_i = Vec::new();
            for j in 0..var_names.len() {
                let indexed_var_name = format!("{}_{}", var_names[j], i);
                matrix_i.push(indexed_var_name.clone());
                matrix_of_expr_i.push(Expr::Var(indexed_var_name));
            }
            matrix.push(matrix_i);
            matrix_of_expr.push(matrix_of_expr_i);
        }
        (matrix_of_expr, matrix)
    }
    /// Creates a 2D indexed variable with format "name_row_col" (e.g., "A_2_3").
    ///
    /// Essential for matrix operations and 2D grid systems in numerical methods.
    ///
    /// # Arguments
    /// * `index_row` - Row index
    /// * `index_col` - Column index  
    /// * `var_name` - Base variable name
    ///
    /// # Returns
    /// Expr::Var with 2D indexed name
    pub fn IndexedVar2D(index_row: usize, index_col: usize, var_name: &str) -> Expr {
        let indexed_var_name = format!("{}_{}_{}", var_name, index_row, index_col);
        Expr::Var(indexed_var_name)
    }

    /// Creates a 2D matrix of indexed variables.
    ///
    /// Generates variables in matrix form like A_0_0, A_0_1, A_1_0, A_1_1, etc.
    /// Returns both the 2D structure and flattened name list.
    ///
    /// # Arguments
    /// * `num_rows` - Number of rows in matrix
    /// * `num_cols` - Number of columns in matrix
    /// * `var_name` - Base name for matrix elements
    ///
    /// # Returns
    /// Tuple of (Vec<Vec<Expr>>, Vec<String>) - matrix structure and flat name list
    pub fn IndexedVars2D(
        num_rows: usize,
        num_cols: usize,
        var_name: &str,
    ) -> (Vec<Vec<Expr>>, Vec<String>) {
        let mut vec_of_names: Vec<String> = Vec::new();
        let matrix = (0..num_rows)
            .map(|i| {
                (0..num_cols)
                    .map(|j| {
                        let indexed_var_name = format!("{}_{}_{}", var_name, i, j);

                        Expr::Var(indexed_var_name)
                    })
                    .collect::<Vec<_>>()
            })
            .collect::<Vec<_>>();
        for i in 0..num_rows {
            for j in 0..num_cols {
                let indexed_var_name = format!("{}_{}_{}", var_name, i, j);
                vec_of_names.push(indexed_var_name);
            }
        }
        (matrix, vec_of_names)
    }

    /// Creates a flattened vector of 2D indexed variables.
    ///
    /// More memory-efficient than IndexedVars2D when matrix structure isn't needed.
    ///
    /// # Arguments
    /// * `num_rows` - Number of rows
    /// * `num_cols` - Number of columns
    /// * `var_name` - Base variable name
    ///
    /// # Returns
    /// Flat Vec<Expr> of all matrix elements
    pub fn IndexedVars2Dflat(num_rows: usize, num_cols: usize, var_name: &str) -> Vec<Expr> {
        (0..num_rows)
            .flat_map(|i| (0..num_cols).map(move |j| Expr::IndexedVar2D(i, j, var_name)))
            .collect()
    }

    /// Creates a symbolic polynomial of specified degree.
    ///
    /// Generates polynomial c0 + c1*x + c2*x^2 + ... + cn*x^n with symbolic coefficients.
    /// Useful for curve fitting and polynomial approximation.
    ///
    /// # Arguments
    /// * `degree` - Highest power in polynomial
    /// * `var_name` - Variable name for polynomial argument
    ///
    /// # Returns
    /// Tuple of (polynomial expression, coefficient names)
    pub fn polyval(degree: usize, var_name: &str) -> (Expr, Vec<String>) {
        let mut eq: Expr = Self::Const(0.0);
        let mut unknowns = Vec::new();
        let arg_expr = Self::Var(var_name.to_string());
        for i in 0..degree + 1 {
            let coeff = Self::Var(format!("c{}", i));
            unknowns.push(format!("c{}", i));
            eq = eq + coeff * (arg_expr.clone().pow(Self::Const(i as f64))).simplify_();
        }
        (eq, unknowns)
    }
}

//___________________________________MACROS____________________________________

/// Macro to create symbolic variables from a comma-separated list
/// Usage: symbols!(x, y, z) -> creates variables x, y, z
#[macro_export]
macro_rules! symbols {
    ($($var:ident),+ $(,)?) => {
        {
            let var_names = stringify!($($var),+);
            let vars = Expr::Symbols(var_names);
            let mut iter = vars.into_iter();
            ($(
                {
                    let $var = iter.next().unwrap();
                    $var
                }
            ),+)
        }
    };
}

/// Macro to create indexed variables
/// Usage: indexed_vars!(5, "x") -> creates x0, x1, x2, x3, x4
#[macro_export]
macro_rules! indexed_vars {
    ($count:expr, $name:expr) => {
        Expr::IndexedVars($count, $name)
    };
}

/// Macro to create 2D indexed variables
/// Usage: indexed_vars_2d!(3, 3, "A") -> creates A_0_0, A_0_1, etc.
#[macro_export]
macro_rules! indexed_vars_2d {
    ($rows:expr, $cols:expr, $name:expr) => {
        Expr::IndexedVars2D($rows, $cols, $name)
    };
}

/// Macro to create a single indexed variable
/// Usage: indexed_var!(5, "x") -> creates x5
#[macro_export]
macro_rules! indexed_var {
    ($index:expr, $name:expr) => {
        Expr::IndexedVar($index, $name)
    };
}

/// Macro to create a 2D indexed variable
/// Usage: indexed_var_2d!(2, 3, "A") -> creates A_2_3
#[macro_export]
macro_rules! indexed_var_2d {
    ($row:expr, $col:expr, $name:expr) => {
        Expr::IndexedVar2D($row, $col, $name)
    };
}