use crate::environment::Environment;
use crate::exact::ExactNum;
use crate::multipoly::MultiPoly;
use crate::node::Node;
use crate::polynomial::Polynomial;
use num_integer::Integer;
use num_traits::{Signed, ToPrimitive};
use std::collections::HashMap;

/// Extract square factors from `n` so that `√n = outside · √inside` with `inside` square-free.
fn extract_square_factors(mut n: u64) -> (u64, u64) {
    if n == 0 {
        return (0, 0);
    }
    let mut outside = 1u64;
    let mut d = 2u64;
    while d * d <= n {
        while n.is_multiple_of(d * d) {
            outside *= d;
            n /= d * d;
        }
        d += 1;
    }
    (outside, n)
}

pub trait Simplifiable {
    fn simplify(&self, env: &Environment) -> Result<Node, String>;
}

impl Simplifiable for Node {
    fn simplify(&self, env: &Environment) -> Result<Node, String> {
        match self {
            Node::Add(left, right) => {
                let left_simplified = left.simplify(env)?;
                let right_simplified = right.simplify(env)?;

                if let (Node::Num(ref l), Node::Num(ref r)) = (&left_simplified, &right_simplified)
                {
                    return Ok(Node::Num(l + r));
                }

                if let Node::Num(ref n) = left_simplified {
                    if n.is_zero() {
                        return Ok(right_simplified);
                    }
                }
                if let Node::Num(ref n) = right_simplified {
                    if n.is_zero() {
                        return Ok(left_simplified);
                    }
                }

                // sin²(x) + cos²(x) → 1
                if let Some(result) = try_pythagorean(&left_simplified, &right_simplified) {
                    return Ok(result);
                }

                // a/d + b/d → (a+b)/d
                if let Some(combined) =
                    try_combine_fractions(&left_simplified, &right_simplified, false, env)
                {
                    return Ok(combined);
                }

                // a√X + b√X → (a+b)√X
                if let Some(combined) =
                    try_combine_like_radicals(&left_simplified, &right_simplified, false, env)
                {
                    return Ok(combined);
                }

                let result = Node::Add(Box::new(left_simplified), Box::new(right_simplified));
                let mut term_map: HashMap<String, ExactNum> = HashMap::new();
                if collect_terms(&result, &mut term_map, env).is_ok() {
                    Ok(rebuild_expression(term_map))
                } else if let Some(normalized) = try_polynomial_normalize(&result) {
                    Ok(normalized)
                } else if let Some(normalized) = try_rational_normalize(&result, env) {
                    Ok(normalized)
                } else {
                    Ok(result)
                }
            }
            Node::Num(_) => {
                // ExactNum::Rational is already in lowest terms (BigRational handles reduction).
                // Nothing to simplify for plain numbers.
                Ok(self.clone())
            }
            Node::Multiply(left, right) => {
                let left_simplified = left.simplify(env)?;
                let right_simplified = right.simplify(env)?;

                // Handle multiplication by zero
                if let Node::Num(ref n) = left_simplified {
                    if n.is_zero() {
                        return Ok(Node::Num(ExactNum::zero()));
                    }
                }
                if let Node::Num(ref n) = right_simplified {
                    if n.is_zero() {
                        return Ok(Node::Num(ExactNum::zero()));
                    }
                }

                // Multiplying by one
                if let Node::Num(ref n) = left_simplified {
                    if n.is_one() {
                        return Ok(right_simplified);
                    }
                }
                if let Node::Num(ref n) = right_simplified {
                    if n.is_one() {
                        return Ok(left_simplified);
                    }
                }

                // If both are numbers, multiply them directly
                if let (Node::Num(ref l), Node::Num(ref r)) = (&left_simplified, &right_simplified)
                {
                    return Ok(Node::Num(l * r));
                }

                // k * (-f) → (-k) * f — absorb negation into coefficient
                if let Node::Num(ref k) = left_simplified {
                    if let Node::Negate(inner) = right_simplified {
                        return Node::Multiply(Box::new(Node::Num(-k.clone())), inner)
                            .simplify(env);
                    }
                }
                // (-f) * k → (-k) * f
                if let Node::Negate(inner) = &left_simplified {
                    if let Node::Num(ref k) = right_simplified {
                        return Node::Multiply(Box::new(Node::Num(-k.clone())), inner.clone())
                            .simplify(env);
                    }
                }
                // f * (-g) → -(f * g), (-f) * g → -(f * g) — general negation extraction
                if let Node::Negate(inner) = &right_simplified {
                    return Node::Negate(Box::new(Node::Multiply(
                        Box::new(left_simplified),
                        inner.clone(),
                    )))
                    .simplify(env);
                }
                if let Node::Negate(inner) = &left_simplified {
                    return Node::Negate(Box::new(Node::Multiply(
                        inner.clone(),
                        Box::new(right_simplified),
                    )))
                    .simplify(env);
                }

                // **Handle implicit multiplication of number and variable (e.g., 5 * x -> 5x)**
                if let (Node::Num(ref l_coef), Node::Variable(ref var)) =
                    (&left_simplified, &right_simplified)
                {
                    return Ok(Node::Multiply(
                        Box::new(Node::Num(l_coef.clone())),
                        Box::new(Node::Variable(var.clone())),
                    ));
                }
                if let (Node::Variable(ref var), Node::Num(ref r_coef)) =
                    (&left_simplified, &right_simplified)
                {
                    return Ok(Node::Multiply(
                        Box::new(Node::Num(r_coef.clone())),
                        Box::new(Node::Variable(var.clone())),
                    ));
                }

                // x^a * x^b → x^(a+b)
                if let (Node::Power(ref base1, ref exp1), Node::Power(ref base2, ref exp2)) =
                    (&left_simplified, &right_simplified)
                {
                    if base1 == base2 {
                        if let (Node::Num(ref a), Node::Num(ref b)) = (exp1.as_ref(), exp2.as_ref())
                        {
                            return Ok(Node::Power(base1.clone(), Box::new(Node::Num(a + b))));
                        }
                    }
                }

                // x * x^a → x^(a+1)
                if let Node::Power(ref base, ref exp) = right_simplified {
                    if *base.as_ref() == left_simplified {
                        if let Node::Num(ref a) = exp.as_ref() {
                            return Ok(Node::Power(
                                base.clone(),
                                Box::new(Node::Num(a + &ExactNum::one())),
                            ));
                        }
                    }
                }
                // x^a * x → x^(a+1)
                if let Node::Power(ref base, ref exp) = left_simplified {
                    if *base.as_ref() == right_simplified {
                        if let Node::Num(ref a) = exp.as_ref() {
                            return Ok(Node::Power(
                                base.clone(),
                                Box::new(Node::Num(a + &ExactNum::one())),
                            ));
                        }
                    }
                }

                // x * x → x^2
                if left_simplified == right_simplified && !matches!(left_simplified, Node::Num(_)) {
                    return Ok(Node::Power(
                        Box::new(left_simplified),
                        Box::new(Node::Num(ExactNum::two())),
                    ));
                }

                let result = Node::Multiply(Box::new(left_simplified), Box::new(right_simplified));
                if let Some(normalized) = try_polynomial_normalize(&result) {
                    Ok(normalized)
                } else if let Some(normalized) = try_rational_normalize(&result, env) {
                    Ok(normalized)
                } else {
                    Ok(result)
                }
            }
            Node::Power(base, exponent) => {
                let base_simplified = base.simplify(env)?;
                let exponent_simplified = exponent.simplify(env)?;

                // 0^n → 0 for n > 0, 1^n → 1
                if let Node::Num(ref b) = base_simplified {
                    if b.is_zero() {
                        if let Node::Num(ref e) = exponent_simplified {
                            if !e.is_negative() {
                                return Ok(Node::Num(ExactNum::zero()));
                            }
                        }
                    }
                    if b.is_one() {
                        return Ok(Node::Num(ExactNum::one()));
                    }
                }

                // x^0 → 1
                if let Node::Num(ref n) = exponent_simplified {
                    if n.is_zero() {
                        return Ok(Node::Num(ExactNum::one()));
                    }
                }

                // x^1 → x
                if let Node::Num(ref n) = exponent_simplified {
                    if n.is_one() {
                        return Ok(base_simplified);
                    }
                }

                // If both the base and exponent are numbers, evaluate the power
                if let (Node::Num(ref b), Node::Num(ref e)) =
                    (&base_simplified, &exponent_simplified)
                {
                    return Ok(Node::Num(b.powf(e)));
                }

                // (x^a)^b → x^(a*b) when both exponents are numeric
                if let Node::Power(inner_base, inner_exp) = &base_simplified {
                    if let (Node::Num(ref a), Node::Num(ref b)) =
                        (&**inner_exp, &exponent_simplified)
                    {
                        return Ok(Node::Power(inner_base.clone(), Box::new(Node::Num(a * b))));
                    }
                }

                // (-1)^(2n) → 1 when n is integer (even exponent of -1)
                if is_neg_one(&base_simplified) && is_even_integer_expr(&exponent_simplified, env) {
                    return Ok(Node::Num(ExactNum::one()));
                }

                Ok(Node::Power(
                    Box::new(base_simplified),
                    Box::new(exponent_simplified),
                ))
            }
            Node::Subtract(left, right) => {
                let left_simplified = left.simplify(env)?;
                let right_simplified = right.simplify(env)?;

                if let (Node::Num(ref l), Node::Num(ref r)) = (&left_simplified, &right_simplified)
                {
                    return Ok(Node::Num(l - r));
                }

                if let Node::Num(ref n) = right_simplified {
                    if n.is_zero() {
                        return Ok(left_simplified);
                    }
                }
                if let Node::Num(ref n) = left_simplified {
                    if n.is_zero() {
                        return Ok(Node::Negate(Box::new(right_simplified)));
                    }
                }

                // 1 - sin²(x) → cos²(x), 1 - cos²(x) → sin²(x)
                if let Node::Num(ref n) = left_simplified {
                    if n.is_one() {
                        if let Some(args) = is_trig_squared(&right_simplified, "sin") {
                            return Ok(Node::Power(
                                Box::new(Node::Function("cos".to_string(), args)),
                                Box::new(Node::Num(ExactNum::two())),
                            ));
                        }
                        if let Some(args) = is_trig_squared(&right_simplified, "cos") {
                            return Ok(Node::Power(
                                Box::new(Node::Function("sin".to_string(), args)),
                                Box::new(Node::Num(ExactNum::two())),
                            ));
                        }
                    }
                }

                // sin²(x) - 1 → -cos²(x), cos²(x) - 1 → -sin²(x)
                if let Node::Num(ref n) = right_simplified {
                    if n.is_one() {
                        if let Some(args) = is_trig_squared(&left_simplified, "sin") {
                            return Ok(Node::Negate(Box::new(Node::Power(
                                Box::new(Node::Function("cos".to_string(), args)),
                                Box::new(Node::Num(ExactNum::two())),
                            ))));
                        }
                        if let Some(args) = is_trig_squared(&left_simplified, "cos") {
                            return Ok(Node::Negate(Box::new(Node::Power(
                                Box::new(Node::Function("sin".to_string(), args)),
                                Box::new(Node::Num(ExactNum::two())),
                            ))));
                        }
                    }
                }

                // a/d - b/d → (a-b)/d
                if let Some(combined) =
                    try_combine_fractions(&left_simplified, &right_simplified, true, env)
                {
                    return Ok(combined);
                }

                // a√X - b√X → (a-b)√X
                if let Some(combined) =
                    try_combine_like_radicals(&left_simplified, &right_simplified, true, env)
                {
                    return Ok(combined);
                }

                let result = Node::Subtract(Box::new(left_simplified), Box::new(right_simplified));
                let mut term_map: HashMap<String, ExactNum> = HashMap::new();
                if collect_terms(&result, &mut term_map, env).is_ok() {
                    Ok(rebuild_expression(term_map))
                } else if let Some(normalized) = try_polynomial_normalize(&result) {
                    Ok(normalized)
                } else if let Some(normalized) = try_rational_normalize(&result, env) {
                    Ok(normalized)
                } else {
                    Ok(result)
                }
            }
            Node::Negate(operand) => {
                let simplified = operand.simplify(env)?;
                if let Node::Num(ref n) = simplified {
                    return Ok(Node::Num(-n.clone()));
                }
                if let Node::Negate(inner) = simplified {
                    return Ok(*inner);
                }
                // -(a + b) → (-a) - b, -(a - b) → b - a
                if let Node::Add(a, b) = simplified {
                    return Node::Subtract(Box::new(Node::Negate(a)), b).simplify(env);
                }
                if let Node::Subtract(a, b) = simplified {
                    return Node::Subtract(b, a).simplify(env);
                }
                Ok(Node::Negate(Box::new(simplified)))
            }
            Node::Divide(left, right) => {
                let left_simplified = left.simplify(env)?;
                let right_simplified = right.simplify(env)?;

                // Cancel common leading negative signs in fractions
                if has_leading_negative(&left_simplified) && has_leading_negative(&right_simplified)
                {
                    let pos_num = negate_leading(&left_simplified);
                    let pos_den = negate_leading(&right_simplified);
                    return Node::Divide(Box::new(pos_num), Box::new(pos_den)).simplify(env);
                }

                if let Node::Num(ref n) = right_simplified {
                    if n.is_one() {
                        return Ok(left_simplified);
                    }
                }

                if let (Node::Num(ref l), Node::Num(ref r)) = (&left_simplified, &right_simplified)
                {
                    return Ok(Node::Num(l / r));
                }

                // (n/expr) / m → (n/m) / expr — collapse nested numeric divisions
                if let Node::Num(ref r) = right_simplified {
                    if let Node::Divide(ref inner_num, ref inner_den) = left_simplified {
                        if let Node::Num(ref l) = **inner_num {
                            let combined = l / r;
                            return Ok(Node::Divide(
                                Box::new(Node::Num(combined)),
                                inner_den.clone(),
                            ));
                        }
                    }
                }

                // (k·expr) / m → (k/m)·expr — cancel common numeric factors
                // Also: -(k·expr) / m → -((k/m)·expr)
                if let Node::Num(ref denom_num) = right_simplified {
                    if !denom_num.is_zero() {
                        // Handle -(k·expr) / m
                        if let Node::Negate(ref inner) = left_simplified {
                            if let Node::Multiply(_, _) = **inner {
                                let inner_div =
                                    Node::Divide(inner.clone(), Box::new(right_simplified.clone()));
                                let inner_result = inner_div.simplify(env)?;
                                return Ok(Node::Negate(Box::new(inner_result)));
                            }
                        }
                        if let Node::Multiply(ref ml, ref mr) = left_simplified {
                            if let Node::Num(ref numer_coeff) = **ml {
                                let reduced = numer_coeff / denom_num;
                                if reduced.is_one() {
                                    return Ok(*mr.clone());
                                }
                                let neg_one = ExactNum::integer(-1);
                                if reduced == neg_one {
                                    return Ok(Node::Negate(mr.clone()));
                                }
                                return Node::Multiply(Box::new(Node::Num(reduced)), mr.clone())
                                    .simplify(env);
                            }
                            if let Node::Num(ref numer_coeff) = **mr {
                                let reduced = numer_coeff / denom_num;
                                if reduced.is_one() {
                                    return Ok(*ml.clone());
                                }
                                let neg_one = ExactNum::integer(-1);
                                if reduced == neg_one {
                                    return Ok(Node::Negate(ml.clone()));
                                }
                                return Node::Multiply(ml.clone(), Box::new(Node::Num(reduced)))
                                    .simplify(env);
                            }
                        }
                    }
                }

                // k / (m·expr) → (k/m) / expr — cancel numeric factor in denominator
                if let Node::Num(ref numer_num) = left_simplified {
                    if let Node::Multiply(ref dl, ref dr) = right_simplified {
                        if let Node::Num(ref denom_coeff) = **dl {
                            if !denom_coeff.is_zero() {
                                let reduced = numer_num / denom_coeff;
                                return Node::Divide(Box::new(Node::Num(reduced)), dr.clone())
                                    .simplify(env);
                            }
                        }
                        if let Node::Num(ref denom_coeff) = **dr {
                            if !denom_coeff.is_zero() {
                                let reduced = numer_num / denom_coeff;
                                return Node::Divide(Box::new(Node::Num(reduced)), dl.clone())
                                    .simplify(env);
                            }
                        }
                    }
                }

                // (k·expr1) / (m·expr2) → cancel common integer factor
                if let Node::Multiply(ref nl, ref nr) = left_simplified {
                    if let Node::Multiply(ref dl, ref dr) = right_simplified {
                        let num_coeff = if let Node::Num(ref k) = **nl {
                            Some((k, nr.as_ref()))
                        } else if let Node::Num(ref k) = **nr {
                            Some((k, nl.as_ref()))
                        } else {
                            None
                        };
                        let den_coeff = if let Node::Num(ref m) = **dl {
                            Some((m, dr.as_ref()))
                        } else if let Node::Num(ref m) = **dr {
                            Some((m, dl.as_ref()))
                        } else {
                            None
                        };
                        if let (Some((k, e1)), Some((m, e2))) = (num_coeff, den_coeff) {
                            if let (ExactNum::Rational(ref kr), ExactNum::Rational(ref mr)) = (k, m)
                            {
                                if kr.is_integer() && mr.is_integer() {
                                    let ki = kr.to_integer();
                                    let mi = mr.to_integer();
                                    let g = ki.gcd(&mi);
                                    let abs_g = if g.is_negative() { -&g } else { g.clone() };
                                    if abs_g > num_bigint::BigInt::from(1) {
                                        let new_k = &ki / &g;
                                        let new_m = &mi / &g;
                                        let one = num_bigint::BigInt::from(1);
                                        let neg_one = num_bigint::BigInt::from(-1);
                                        let make_num = |n: num_bigint::BigInt| {
                                            use num_rational::BigRational;
                                            ExactNum::Rational(BigRational::from_integer(n))
                                        };
                                        let new_num = if new_k == one {
                                            e1.clone()
                                        } else if new_k == neg_one {
                                            Node::Negate(Box::new(e1.clone()))
                                        } else {
                                            Node::Multiply(
                                                Box::new(Node::Num(make_num(new_k))),
                                                Box::new(e1.clone()),
                                            )
                                        };
                                        let new_den = if new_m == one {
                                            e2.clone()
                                        } else if new_m == neg_one {
                                            Node::Negate(Box::new(e2.clone()))
                                        } else {
                                            Node::Multiply(
                                                Box::new(Node::Num(make_num(new_m))),
                                                Box::new(e2.clone()),
                                            )
                                        };
                                        return Node::Divide(Box::new(new_num), Box::new(new_den))
                                            .simplify(env);
                                    }
                                }
                            }
                        }
                    }
                }

                // x / x → 1
                if left_simplified == right_simplified && !matches!(left_simplified, Node::Num(_)) {
                    return Ok(Node::Num(ExactNum::one()));
                }

                // sin(x) / cos(x) → tan(x), cos(x) / sin(x) → cot(x)
                if let (
                    Node::Function(ref fname1, ref args1),
                    Node::Function(ref fname2, ref args2),
                ) = (&left_simplified, &right_simplified)
                {
                    if fname1 == "sin" && fname2 == "cos" && args1 == args2 {
                        return Ok(Node::Function("tan".to_string(), args1.clone()));
                    }
                    if fname1 == "cos" && fname2 == "sin" && args1 == args2 {
                        return Ok(Node::Function("cot".to_string(), args1.clone()));
                    }
                }

                // k·sin(x) / cos(x) → k·tan(x), k·cos(x) / sin(x) → k·cot(x)
                if let Node::Multiply(ref ml, ref mr) = left_simplified {
                    if let Node::Function(ref fname2, ref args2) = right_simplified {
                        let (coeff, func) = if matches!(**ml, Node::Function(_, _)) {
                            (mr, ml)
                        } else {
                            (ml, mr)
                        };
                        if let Node::Function(ref fname1, ref args1) = **func {
                            if args1 == args2 {
                                if fname1 == "sin" && fname2 == "cos" {
                                    return Node::Multiply(
                                        coeff.clone(),
                                        Box::new(Node::Function("tan".to_string(), args1.clone())),
                                    )
                                    .simplify(env);
                                }
                                if fname1 == "cos" && fname2 == "sin" {
                                    return Node::Multiply(
                                        coeff.clone(),
                                        Box::new(Node::Function("cot".to_string(), args1.clone())),
                                    )
                                    .simplify(env);
                                }
                            }
                        }
                    }
                }

                // 1 / sin(x) → csc(x), 1 / cos(x) → sec(x), 1 / tan(x) → cot(x)
                if let Node::Num(ref n) = left_simplified {
                    if n.is_one() {
                        if let Node::Function(ref fname, ref args) = right_simplified {
                            let recip = match fname.as_str() {
                                "sin" => Some("csc"),
                                "cos" => Some("sec"),
                                "tan" => Some("cot"),
                                _ => None,
                            };
                            if let Some(recip_name) = recip {
                                return Ok(Node::Function(recip_name.to_string(), args.clone()));
                            }
                        }
                    }
                }

                // x^a / x^b → x^(a-b)
                if let (Node::Power(ref base1, ref exp1), Node::Power(ref base2, ref exp2)) =
                    (&left_simplified, &right_simplified)
                {
                    if base1 == base2 {
                        if let (Node::Num(ref a), Node::Num(ref b)) = (exp1.as_ref(), exp2.as_ref())
                        {
                            let diff = a - b;
                            if diff.is_zero() {
                                return Ok(Node::Num(ExactNum::one()));
                            } else if diff.is_one() {
                                return Ok(*base1.clone());
                            }
                            return Ok(Node::Power(base1.clone(), Box::new(Node::Num(diff))));
                        }
                    }
                }

                // x^a / x → x^(a-1)
                if let Node::Power(ref base, ref exp) = left_simplified {
                    if *base.as_ref() == right_simplified {
                        if let Node::Num(ref a) = exp.as_ref() {
                            let diff = a - &ExactNum::one();
                            if diff.is_zero() {
                                return Ok(Node::Num(ExactNum::one()));
                            } else if diff.is_one() {
                                return Ok(*base.clone());
                            }
                            return Ok(Node::Power(base.clone(), Box::new(Node::Num(diff))));
                        }
                    }
                }

                // x / x^a → x^(1-a)
                if let Node::Power(ref base, ref exp) = right_simplified {
                    if *base.as_ref() == left_simplified {
                        if let Node::Num(ref a) = exp.as_ref() {
                            let diff = &ExactNum::one() - a;
                            if diff.is_zero() {
                                return Ok(Node::Num(ExactNum::one()));
                            } else if diff.is_one() {
                                return Ok(*base.clone());
                            }
                            return Ok(Node::Power(base.clone(), Box::new(Node::Num(diff))));
                        }
                    }
                }

                if let Some(simplified) = try_polynomial_divide(&left_simplified, &right_simplified)
                {
                    return Ok(simplified);
                }

                Ok(Node::Divide(
                    Box::new(left_simplified),
                    Box::new(right_simplified),
                ))
            }

            Node::Summation(index_var, start, end, body) => {
                let start_simplified = start.simplify(env)?;
                let end_simplified = end.simplify(env)?;
                let body_simplified = body.simplify(env)?;

                // Try to evaluate if bounds are constant values
                if let (Node::Num(ref start_n), Node::Num(ref end_n)) =
                    (&start_simplified, &end_simplified)
                {
                    if start_n.is_integer() && end_n.is_integer() {
                        let start_val = start_n.to_f64();
                        let end_val = end_n.to_f64();

                        // For small ranges (fewer than 100 terms), we can expand inline
                        let range_size = (end_val - start_val + 1.0) as usize;
                        if range_size <= 10 {
                            let mut sum_node = Node::Num(ExactNum::zero());

                            // Create a temporary environment for each iteration
                            let mut sum_env = env.clone();

                            let start_i = start_val as i64;
                            let end_i = end_val as i64;

                            // Evaluate each term and add them together
                            for i in start_i..=end_i {
                                sum_env.set_exact(index_var, ExactNum::integer(i));

                                // Create a substituted body for this iteration
                                let substituted_body = crate::substitute::substitute_variable(
                                    &body_simplified,
                                    index_var,
                                    &Node::Num(ExactNum::integer(i)),
                                )?;

                                // Add this term to our running sum
                                sum_node =
                                    Node::Add(Box::new(sum_node), Box::new(substituted_body));
                            }

                            return Ok(sum_node);
                        }
                    }
                }

                // If we can't or shouldn't evaluate the summation, return it with simplified components
                Ok(Node::Summation(
                    index_var.clone(),
                    Box::new(start_simplified),
                    Box::new(end_simplified),
                    Box::new(body_simplified),
                ))
            }
            Node::Abs(operand) => {
                let simplified = operand.simplify(env)?;
                if let Node::Num(ref n) = simplified {
                    return Ok(Node::Num(n.abs()));
                }
                // |x| → x when x is nonnegative
                if let Node::Variable(ref v) = simplified {
                    if env.assumptions().is_nonneg(v) {
                        return Ok(simplified);
                    }
                    if env.assumptions().is_negative(v) {
                        return Ok(Node::Negate(Box::new(simplified)));
                    }
                }
                // |-x| → |x|
                if let Node::Negate(inner) = simplified {
                    return Ok(Node::Abs(inner));
                }
                // ||x|| → |x|
                if let Node::Abs(_) = simplified {
                    return Ok(simplified);
                }
                Ok(Node::Abs(Box::new(simplified)))
            }
            Node::Sqrt(operand) => {
                let simplified = operand.simplify(env)?;
                if let Node::Num(ref n) = simplified {
                    let s = n.sqrt();
                    if matches!(s, ExactNum::Rational(_)) {
                        return Ok(Node::Num(s));
                    }
                    // Non-perfect-square integer: extract square factors
                    if let Some(val) = n.to_i64() {
                        if val > 0 {
                            let (outside, inside) = extract_square_factors(val as u64);
                            if inside == 1 {
                                return Ok(Node::Num(ExactNum::integer(outside as i64)));
                            }
                            let sqrt_inside =
                                Node::Sqrt(Box::new(Node::Num(ExactNum::integer(inside as i64))));
                            if outside == 1 {
                                return Ok(sqrt_inside);
                            }
                            return Ok(Node::Multiply(
                                Box::new(Node::Num(ExactNum::integer(outside as i64))),
                                Box::new(sqrt_inside),
                            ));
                        }
                    }
                    // Non-perfect-square rational: try numerator and denominator separately
                    if let ExactNum::Rational(ref r) = n {
                        if !r.is_negative() {
                            if let (Some(num), Some(den)) = (r.numer().to_i64(), r.denom().to_i64())
                            {
                                if num > 0 && den > 0 {
                                    let (num_out, num_in) = extract_square_factors(num as u64);
                                    let (den_out, den_in) = extract_square_factors(den as u64);
                                    if num_in == 1 && den_in == 1 {
                                        return Ok(Node::Num(ExactNum::rational(
                                            num_out as i64,
                                            den_out as i64,
                                        )));
                                    }
                                    let sqrt_part = Node::Sqrt(Box::new(Node::Num(
                                        ExactNum::rational(num_in as i64, den_in as i64),
                                    )));
                                    if num_out == 1 && den_out == 1 {
                                        return Ok(sqrt_part);
                                    }
                                    return Ok(Node::Multiply(
                                        Box::new(Node::Num(ExactNum::rational(
                                            num_out as i64,
                                            den_out as i64,
                                        ))),
                                        Box::new(sqrt_part),
                                    ));
                                }
                            }
                        }
                    }
                    // Keep symbolic — do NOT fall back to float
                    return Ok(Node::Sqrt(Box::new(simplified)));
                }
                // Product radicand: decompose and extract
                if let Some(result) = try_simplify_sqrt_product(&simplified, env, true) {
                    return Ok(result);
                }
                // sqrt(x²) → x when x positive, |x| otherwise
                if let Node::Power(ref base, ref exp) = simplified {
                    if let Node::Num(ref e) = **exp {
                        if e == &ExactNum::two() {
                            if let Node::Variable(ref v) = **base {
                                if env.assumptions().is_nonneg(v) {
                                    return Ok(*base.clone());
                                }
                            }
                            return Ok(Node::Abs(base.clone()));
                        }
                    }
                }
                Ok(Node::Sqrt(Box::new(simplified)))
            }
            Node::Function(name, args) => {
                let simplified_args: Vec<Node> = args
                    .iter()
                    .map(|a| a.simplify(env))
                    .collect::<Result<Vec<_>, _>>()?;

                if simplified_args.len() == 1 {
                    let arg = &simplified_args[0];
                    match name.as_str() {
                        "ln" => {
                            // ln(e^x) → x
                            if let Node::Power(base, exp) = arg {
                                if let Node::Num(ref b) = **base {
                                    if (b.to_f64() - std::f64::consts::E).abs() < 1e-14 {
                                        return Ok(*exp.clone());
                                    }
                                }
                                // ln(a^b) → b·ln(a), then re-simplify since ln(a)
                                // may itself expand (e.g. a = x·y → ln(x)+ln(y))
                                let inner_ln =
                                    Node::Function("ln".to_string(), vec![*base.clone()])
                                        .simplify(env)?;
                                return Node::Multiply(exp.clone(), Box::new(inner_ln))
                                    .simplify(env);
                            }
                            // ln(a·b) → ln(a) + ln(b), re-simplify each ln
                            if let Node::Multiply(a, b) = arg {
                                let ln_a = Node::Function("ln".to_string(), vec![*a.clone()])
                                    .simplify(env)?;
                                let ln_b = Node::Function("ln".to_string(), vec![*b.clone()])
                                    .simplify(env)?;
                                return Node::Add(Box::new(ln_a), Box::new(ln_b)).simplify(env);
                            }
                            // ln(a/b) → ln(a) - ln(b), re-simplify each ln
                            if let Node::Divide(a, b) = arg {
                                let ln_a = Node::Function("ln".to_string(), vec![*a.clone()])
                                    .simplify(env)?;
                                let ln_b = Node::Function("ln".to_string(), vec![*b.clone()])
                                    .simplify(env)?;
                                return Node::Subtract(Box::new(ln_a), Box::new(ln_b))
                                    .simplify(env);
                            }
                        }
                        "exp" => {
                            // exp(ln(x)) → x
                            if let Node::Function(inner_name, inner_args) = arg {
                                if inner_name == "ln" && inner_args.len() == 1 {
                                    return Ok(inner_args[0].clone());
                                }
                            }
                        }
                        "sqrt" => {
                            // Numeric: extract square factors to avoid float fallback
                            if let Node::Num(ref n) = arg {
                                // Perfect rational square → exact result
                                let s = n.sqrt();
                                if matches!(s, ExactNum::Rational(_)) {
                                    return Ok(Node::Num(s));
                                }
                                // Non-perfect-square integer: √n = outside · √inside
                                if let Some(val) = n.to_i64() {
                                    if val > 0 {
                                        let (outside, inside) = extract_square_factors(val as u64);
                                        if inside == 1 {
                                            return Ok(Node::Num(ExactNum::integer(
                                                outside as i64,
                                            )));
                                        }
                                        let sqrt_inside = Node::Sqrt(Box::new(Node::Num(
                                            ExactNum::integer(inside as i64),
                                        )));
                                        if outside == 1 {
                                            return Ok(sqrt_inside);
                                        }
                                        return Ok(Node::Multiply(
                                            Box::new(Node::Num(ExactNum::integer(outside as i64))),
                                            Box::new(sqrt_inside),
                                        ));
                                    }
                                }
                                // Non-perfect-square rational
                                if let ExactNum::Rational(ref r) = n {
                                    if !r.is_negative() {
                                        if let (Some(num), Some(den)) =
                                            (r.numer().to_i64(), r.denom().to_i64())
                                        {
                                            if num > 0 && den > 0 {
                                                let (num_out, num_in) =
                                                    extract_square_factors(num as u64);
                                                let (den_out, den_in) =
                                                    extract_square_factors(den as u64);
                                                if num_in == 1 && den_in == 1 {
                                                    return Ok(Node::Num(ExactNum::rational(
                                                        num_out as i64,
                                                        den_out as i64,
                                                    )));
                                                }
                                                let sqrt_part = Node::Sqrt(Box::new(Node::Num(
                                                    ExactNum::rational(
                                                        num_in as i64,
                                                        den_in as i64,
                                                    ),
                                                )));
                                                if num_out == 1 && den_out == 1 {
                                                    return Ok(sqrt_part);
                                                }
                                                return Ok(Node::Multiply(
                                                    Box::new(Node::Num(ExactNum::rational(
                                                        num_out as i64,
                                                        den_out as i64,
                                                    ))),
                                                    Box::new(sqrt_part),
                                                ));
                                            }
                                        }
                                    }
                                }
                                // Fallback: keep symbolic as Sqrt node
                                return Ok(Node::Sqrt(Box::new(arg.clone())));
                            }
                            // Product radicand: decompose and extract
                            if let Some(result) = try_simplify_sqrt_product(arg, env, false) {
                                return Ok(result);
                            }
                            // sqrt(x²) → x when x nonneg, |x| otherwise
                            if let Node::Power(base, exp) = arg {
                                if let Node::Num(ref e) = **exp {
                                    if e == &ExactNum::two() {
                                        if let Node::Variable(ref v) = **base {
                                            if env.assumptions().is_nonneg(v) {
                                                return Ok(*base.clone());
                                            }
                                        }
                                        return Ok(Node::Abs(base.clone()));
                                    }
                                }
                            }
                        }
                        // sin(-x) → -sin(x)
                        "sin" | "tan" | "sinh" | "tanh" => {
                            if let Node::Negate(inner) = arg {
                                return Ok(Node::Negate(Box::new(Node::Function(
                                    name.clone(),
                                    vec![*inner.clone()],
                                ))));
                            }
                        }
                        // cos(-x) → cos(x)
                        "cos" | "cosh" => {
                            if let Node::Negate(inner) = arg {
                                return Ok(Node::Function(name.clone(), vec![*inner.clone()]));
                            }
                        }
                        _ => {}
                    }
                }

                let all_numeric = simplified_args.iter().all(|a| matches!(a, Node::Num(_)));
                if all_numeric {
                    let f64_args: Vec<f64> = simplified_args
                        .iter()
                        .map(|a| {
                            if let Node::Num(n) = a {
                                n.to_f64()
                            } else {
                                unreachable!()
                            }
                        })
                        .collect();
                    if let Ok(result) = crate::functions::call_function(name, f64_args) {
                        if result.is_finite() {
                            return Ok(Node::Num(ExactNum::from_f64(result)));
                        }
                    }
                }

                Ok(Node::Function(name.clone(), simplified_args))
            }
            _ => Ok(self.clone()),
        }
    }
}

fn has_leading_negative(node: &Node) -> bool {
    match node {
        Node::Negate(_) => true,
        Node::Num(n) => n.is_negative(),
        Node::Multiply(a, _) => has_leading_negative(a),
        Node::Add(a, _) | Node::Subtract(a, _) => has_leading_negative(a),
        _ => false,
    }
}

fn negate_leading(node: &Node) -> Node {
    match node {
        Node::Negate(inner) => *inner.clone(),
        Node::Num(n) => Node::Num(-n.clone()),
        Node::Multiply(a, b) => Node::Multiply(Box::new(negate_leading(a)), b.clone()),
        Node::Subtract(a, b) => Node::Subtract(b.clone(), a.clone()),
        _ => Node::Negate(Box::new(node.clone())),
    }
}

fn collect_terms_inner(
    node: &Node,
    term_map: &mut HashMap<String, ExactNum>,
    sign: &ExactNum,
) -> Result<(), String> {
    match node {
        Node::Add(left, right) => {
            collect_terms_inner(left, term_map, sign)?;
            collect_terms_inner(right, term_map, sign)?;
        }
        Node::Subtract(left, right) => {
            collect_terms_inner(left, term_map, sign)?;
            let neg_sign = sign.clone() * ExactNum::integer(-1);
            collect_terms_inner(right, term_map, &neg_sign)?;
        }
        Node::Negate(inner) => {
            let neg_sign = sign.clone() * ExactNum::integer(-1);
            collect_terms_inner(inner, term_map, &neg_sign)?;
        }
        Node::Multiply(left, right) => {
            if let (Node::Num(ref coef), Node::Variable(ref var)) = (&**left, &**right) {
                let entry = term_map.entry(var.clone()).or_insert_with(ExactNum::zero);
                *entry = entry.clone() + coef.clone() * sign.clone();
            } else {
                return Err("Unsupported multiply form in collect_terms".to_string());
            }
        }
        Node::Variable(var) => {
            let entry = term_map.entry(var.clone()).or_insert_with(ExactNum::zero);
            *entry = entry.clone() + sign.clone();
        }
        Node::Num(num) => {
            let entry = term_map
                .entry("".to_string())
                .or_insert_with(ExactNum::zero);
            *entry = entry.clone() + num.clone() * sign.clone();
        }
        _ => return Err("Unsupported node type in collect_terms".to_string()),
    }
    Ok(())
}

fn collect_terms(
    node: &Node,
    term_map: &mut HashMap<String, ExactNum>,
    _env: &Environment,
) -> Result<(), String> {
    collect_terms_inner(node, term_map, &ExactNum::one())
}

fn rebuild_expression(term_map: HashMap<String, ExactNum>) -> Node {
    let mut terms: Vec<(String, ExactNum)> = term_map.into_iter().collect();

    // Sort: variables alphabetically first, constant term last
    terms.sort_by(|a, b| match (a.0.is_empty(), b.0.is_empty()) {
        (true, false) => std::cmp::Ordering::Greater,
        (false, true) => std::cmp::Ordering::Less,
        _ => a.0.cmp(&b.0),
    });

    // Build (abs_node, is_negative) pairs for non-zero terms
    let mut signed_terms: Vec<(Node, bool)> = vec![];

    for (var, coef) in terms {
        if coef.is_zero() {
            continue;
        }
        let negative = coef.is_negative();
        let abs_coef = if negative {
            -coef.clone()
        } else {
            coef.clone()
        };

        let node = if var.is_empty() {
            Node::Num(abs_coef)
        } else if abs_coef.is_one() {
            Node::Variable(var)
        } else {
            Node::Multiply(Box::new(Node::Num(abs_coef)), Box::new(Node::Variable(var)))
        };
        signed_terms.push((node, negative));
    }

    if signed_terms.is_empty() {
        return Node::Num(ExactNum::zero());
    }

    let (first_node, first_neg) = signed_terms.remove(0);
    let mut result = if first_neg {
        Node::Negate(Box::new(first_node))
    } else {
        first_node
    };

    for (node, negative) in signed_terms {
        result = if negative {
            Node::Subtract(Box::new(result), Box::new(node))
        } else {
            Node::Add(Box::new(result), Box::new(node))
        };
    }

    result
}

fn find_single_variable(node: &Node) -> Option<String> {
    let mut vars = std::collections::HashSet::new();
    collect_variables(node, &mut vars);
    if vars.len() == 1 {
        vars.into_iter().next()
    } else {
        None
    }
}

fn collect_variables(node: &Node, vars: &mut std::collections::HashSet<String>) {
    match node {
        Node::Variable(v) => {
            vars.insert(v.clone());
        }
        Node::Add(l, r)
        | Node::Subtract(l, r)
        | Node::Multiply(l, r)
        | Node::Divide(l, r)
        | Node::Power(l, r) => {
            collect_variables(l, vars);
            collect_variables(r, vars);
        }
        Node::Negate(inner) | Node::Sqrt(inner) | Node::Abs(inner) => {
            collect_variables(inner, vars);
        }
        _ => {}
    }
}

fn contains_divide(node: &Node) -> bool {
    match node {
        Node::Divide(_, _) => true,
        Node::Add(a, b) | Node::Subtract(a, b) | Node::Multiply(a, b) => {
            contains_divide(a) || contains_divide(b)
        }
        Node::Negate(a) => contains_divide(a),
        _ => false,
    }
}

fn try_rational_normalize(node: &Node, env: &Environment) -> Option<Node> {
    if !contains_divide(node) {
        return None;
    }
    let (num, den) = crate::expression::to_rational_form(node)?;
    let num_simplified = num.simplify(env).ok()?;
    let den_simplified = den.simplify(env).ok()?;
    if let Node::Num(ref n) = den_simplified {
        if n.is_one() {
            return Some(num_simplified);
        }
    }
    Some(Node::Divide(
        Box::new(num_simplified),
        Box::new(den_simplified),
    ))
}

fn try_polynomial_normalize(node: &Node) -> Option<Node> {
    if let Some(var) = find_single_variable(node) {
        let poly = Polynomial::from_node(node, &var).ok()?;
        return Some(poly.to_node());
    }
    // Multivariate fallback
    let mp = MultiPoly::from_node(node).ok()?;
    Some(mp.to_node())
}

fn try_factored_display(poly: &Polynomial) -> Option<Node> {
    use crate::mod_poly::factor_over_q;
    use crate::polynomial::{lcm_bigint, rational_to_node};
    use num_bigint::BigInt;
    use num_traits::One;

    // Only factor degree > 1 polynomials
    let deg = poly.degree()?;
    if deg <= 1 {
        return None;
    }

    let (content, factors) = factor_over_q(poly);
    // Only use factored form when there are multiple factors (distinct or repeated)
    if factors.len() <= 1 {
        return None;
    }

    // Group identical monic factors and convert to integer-coefficient form
    let mut grouped: Vec<(Polynomial, usize)> = Vec::new();
    let mut adjusted_content = content;

    for f in &factors {
        // Clear denominators: find LCM of coefficient denominators
        let mut lcm = BigInt::one();
        for i in 0..=f.degree().unwrap_or(0) {
            let c = f.coeff(i);
            lcm = lcm_bigint(&lcm, c.denom());
        }
        let scale = num_rational::BigRational::from_integer(lcm);
        let f_int = f.scalar_mul(&scale);
        adjusted_content /= &scale;

        // Group identical factors
        let found = grouped.iter_mut().find(|(p, _)| {
            p.degree() == f_int.degree()
                && (0..=p.degree().unwrap_or(0)).all(|i| p.coeff(i) == f_int.coeff(i))
        });
        if let Some(entry) = found {
            entry.1 += 1;
        } else {
            grouped.push((f_int, 1));
        }
    }

    // Build the factored node
    let mut parts: Vec<Node> = Vec::new();

    // Add content coefficient if not 1
    if !adjusted_content.is_one() {
        parts.push(rational_to_node(&adjusted_content));
    }

    for (factor, mult) in &grouped {
        let factor_node = factor.to_node();
        let term = if *mult > 1 {
            Node::Power(
                Box::new(factor_node),
                Box::new(Node::Num(ExactNum::integer(*mult as i64))),
            )
        } else {
            factor_node
        };
        parts.push(term);
    }

    if parts.is_empty() {
        return Some(Node::Num(ExactNum::one()));
    }

    let mut result = parts.remove(0);
    for part in parts {
        result = Node::Multiply(Box::new(result), Box::new(part));
    }

    Some(result)
}

fn try_polynomial_divide(numer: &Node, denom: &Node) -> Option<Node> {
    let mut vars = std::collections::HashSet::new();
    collect_variables(numer, &mut vars);
    collect_variables(denom, &mut vars);

    if vars.len() == 1 {
        let var = vars.into_iter().next()?;
        return try_univariate_divide(numer, denom, &var);
    }

    if vars.len() >= 2 {
        return try_multivariate_divide(numer, denom);
    }

    None
}

fn try_univariate_divide(numer: &Node, denom: &Node, var: &str) -> Option<Node> {
    use crate::polynomial::rational_gcd;
    let n = Polynomial::from_node(numer, var).ok()?;
    let d = Polynomial::from_node(denom, var).ok()?;

    if d.is_zero() {
        return None;
    }

    // Step 1: Polynomial GCD cancellation
    let g = n.gcd(&d);
    let (n_red, d_red) = if g.degree().unwrap_or(0) > 0 {
        let (nr, nr_rem) = n.div_rem(&g).ok()?;
        let (dr, dr_rem) = d.div_rem(&g).ok()?;
        if !nr_rem.is_zero() || !dr_rem.is_zero() {
            return None;
        }
        (nr, dr)
    } else {
        (n.clone(), d.clone())
    };

    // Step 2: Content simplification — cancel GCD of rational coefficients
    let c_n = n_red.content();
    let c_d = d_red.content();
    let c_gcd = rational_gcd(&c_n, &c_d);
    let one = num_rational::BigRational::from_integer(num_bigint::BigInt::from(1));
    let (n_final, d_final) = if c_gcd != one {
        let inv = num_rational::BigRational::from_integer(num_bigint::BigInt::from(1)) / &c_gcd;
        (n_red.scalar_mul(&inv), d_red.scalar_mul(&inv))
    } else {
        (n_red, d_red)
    };

    // Step 3: Try factored display for numerator and denominator
    let denom_factored = try_factored_display(&d_final);
    let numer_factored = try_factored_display(&n_final);

    // Check if anything changed (GCD/content simplification or factored display)
    let changed =
        n_final != n || d_final != d || denom_factored.is_some() || numer_factored.is_some();
    if !changed {
        return None;
    }

    // Step 4: Build result
    if d_final.is_constant() {
        let d_val = d_final.coeff(0);
        if d_val == one {
            return Some(numer_factored.unwrap_or_else(|| n_final.to_node()));
        }
    }

    let denom_node = denom_factored.unwrap_or_else(|| d_final.to_node());
    let numer_node = numer_factored.unwrap_or_else(|| n_final.to_node());

    Some(Node::Divide(Box::new(numer_node), Box::new(denom_node)))
}

fn try_multivariate_divide(numer: &Node, denom: &Node) -> Option<Node> {
    let n = MultiPoly::from_node(numer).ok()?;
    let d = MultiPoly::from_node(denom).ok()?;

    if d.is_zero() {
        return None;
    }

    // Step 1: Polynomial GCD cancellation
    let g = MultiPoly::gcd(&n, &d);
    let (n_red, d_red) = if !g.is_constant() {
        (n.exact_div(&g), d.exact_div(&g))
    } else {
        (n.clone(), d.clone())
    };

    // Step 2: Rational content simplification
    use crate::polynomial::rational_gcd;
    use num_traits::One;
    let c_n = n_red.rational_content();
    let c_d = d_red.rational_content();
    let c_gcd = rational_gcd(&c_n, &c_d);
    let (n_final, d_final) = if !c_gcd.is_one() {
        (
            n_red.scalar_div_rational(&c_gcd),
            d_red.scalar_div_rational(&c_gcd),
        )
    } else {
        (n_red, d_red)
    };

    // Step 3: Check if anything changed
    if n_final == n && d_final == d {
        return None;
    }

    if d_final.is_one() {
        return Some(n_final.to_node());
    }
    if let Some(d_val) = d_final.as_constant() {
        if !num_traits::Zero::is_zero(d_val) {
            let inv = num_rational::BigRational::from_integer(num_bigint::BigInt::from(1)) / d_val;
            return Some(n_final.scalar_mul(&inv).to_node());
        }
    }

    Some(Node::Divide(
        Box::new(n_final.to_node()),
        Box::new(d_final.to_node()),
    ))
}

fn is_trig_squared(node: &Node, func_name: &str) -> Option<Vec<Node>> {
    if let Node::Power(base, exp) = node {
        if let Node::Num(ref e) = **exp {
            if e == &ExactNum::two() {
                if let Node::Function(name, args) = base.as_ref() {
                    if name == func_name {
                        return Some(args.clone());
                    }
                }
            }
        }
    }
    None
}

fn try_pythagorean(left: &Node, right: &Node) -> Option<Node> {
    // sin²(x) + cos²(x) → 1
    if let (Some(sin_args), Some(cos_args)) =
        (is_trig_squared(left, "sin"), is_trig_squared(right, "cos"))
    {
        if sin_args == cos_args {
            return Some(Node::Num(ExactNum::one()));
        }
    }
    // cos²(x) + sin²(x) → 1
    if let (Some(cos_args), Some(sin_args)) =
        (is_trig_squared(left, "cos"), is_trig_squared(right, "sin"))
    {
        if cos_args == sin_args {
            return Some(Node::Num(ExactNum::one()));
        }
    }

    // a·sin²(x) + a·cos²(x) → a (with coefficient)
    if let (Some((coeff_l, sin_args)), Some((coeff_r, cos_args))) = (
        extract_coeff_trig_sq(left, "sin"),
        extract_coeff_trig_sq(right, "cos"),
    ) {
        if sin_args == cos_args && coeff_l == coeff_r {
            return Some(Node::Num(coeff_l));
        }
    }
    if let (Some((coeff_l, cos_args)), Some((coeff_r, sin_args))) = (
        extract_coeff_trig_sq(left, "cos"),
        extract_coeff_trig_sq(right, "sin"),
    ) {
        if cos_args == sin_args && coeff_l == coeff_r {
            return Some(Node::Num(coeff_l));
        }
    }

    None
}

fn extract_coeff_trig_sq(node: &Node, func_name: &str) -> Option<(ExactNum, Vec<Node>)> {
    if let Some(args) = is_trig_squared(node, func_name) {
        return Some((ExactNum::one(), args));
    }
    if let Node::Multiply(coeff, power) = node {
        if let Node::Num(ref c) = **coeff {
            if let Some(args) = is_trig_squared(power, func_name) {
                return Some((c.clone(), args));
            }
        }
        if let Node::Num(ref c) = **power {
            if let Some(args) = is_trig_squared(coeff, func_name) {
                return Some((c.clone(), args));
            }
        }
    }
    None
}

fn is_neg_one(node: &Node) -> bool {
    match node {
        Node::Num(n) => n == &ExactNum::integer(-1),
        Node::Negate(inner) => {
            matches!(&**inner, Node::Num(n) if n.is_one())
        }
        _ => false,
    }
}

fn is_even_integer_expr(node: &Node, env: &Environment) -> bool {
    // 2n, 2*n, k*n where k is even and n is integer
    if let Node::Multiply(left, right) = node {
        match (&**left, &**right) {
            (Node::Num(k), Node::Variable(v)) | (Node::Variable(v), Node::Num(k)) => {
                return k.is_even() && env.assumptions().is_integer(v);
            }
            _ => {}
        }
    }
    // A numeric even integer by itself
    if let Node::Num(n) = node {
        return n.is_even();
    }
    false
}

/// Extract (coefficient, radical) from a term that's either a bare √X or coeff·√X.
fn extract_radical_parts(node: &Node) -> Option<(ExactNum, Node)> {
    match node {
        Node::Sqrt(_) => Some((ExactNum::one(), node.clone())),
        Node::Multiply(left, right) => {
            if let Node::Num(ref coeff) = **left {
                if matches!(**right, Node::Sqrt(_)) {
                    return Some((coeff.clone(), *right.clone()));
                }
            }
            if let Node::Num(ref coeff) = **right {
                if matches!(**left, Node::Sqrt(_)) {
                    return Some((coeff.clone(), *left.clone()));
                }
            }
            None
        }
        _ => None,
    }
}

/// Combine like-radical terms: a√X ± b√X → (a±b)√X
fn try_combine_like_radicals(
    left: &Node,
    right: &Node,
    is_subtract: bool,
    _env: &Environment,
) -> Option<Node> {
    let (l_coeff, l_radical) = extract_radical_parts(left)?;
    let (r_coeff, r_radical) = extract_radical_parts(right)?;

    // Check structural equality of radicals
    if format!("{}", l_radical) != format!("{}", r_radical) {
        return None;
    }

    let combined = if is_subtract {
        l_coeff - r_coeff
    } else {
        l_coeff + r_coeff
    };

    if combined.is_zero() {
        return Some(Node::Num(ExactNum::zero()));
    }
    if combined.is_one() {
        return Some(l_radical);
    }
    let neg_one = ExactNum::integer(-1);
    if combined == neg_one {
        return Some(Node::Negate(Box::new(l_radical)));
    }

    Some(Node::Multiply(
        Box::new(Node::Num(combined)),
        Box::new(l_radical),
    ))
}

/// Combine fractions with the same denominator: a/d ± b/d → (a±b)/d
fn try_combine_fractions(
    left: &Node,
    right: &Node,
    subtract: bool,
    env: &Environment,
) -> Option<Node> {
    let left_is_frac = matches!(left, Node::Divide(_, _));
    let right_is_frac = matches!(right, Node::Divide(_, _));
    if !left_is_frac && !right_is_frac {
        return None;
    }

    let (ln, ld) = match left {
        Node::Divide(n, d) => (n.as_ref().clone(), d.as_ref().clone()),
        _ => (left.clone(), Node::Num(ExactNum::one())),
    };
    let (rn, rd) = match right {
        Node::Divide(n, d) => (n.as_ref().clone(), d.as_ref().clone()),
        _ => (right.clone(), Node::Num(ExactNum::one())),
    };

    if ld == rd {
        let combined_num = if subtract {
            Node::Subtract(Box::new(ln), Box::new(rn))
        } else {
            Node::Add(Box::new(ln), Box::new(rn))
        };
        let simplified_num = combined_num.simplify(env).ok()?;
        let result = Node::Divide(Box::new(simplified_num), Box::new(ld));
        return result.simplify(env).ok();
    }

    let left_num = Node::Multiply(Box::new(ln), Box::new(rd.clone()));
    let right_num = Node::Multiply(Box::new(rn), Box::new(ld.clone()));
    let combined_num = if subtract {
        Node::Subtract(Box::new(left_num), Box::new(right_num))
    } else {
        Node::Add(Box::new(left_num), Box::new(right_num))
    };
    let combined_den = Node::Multiply(Box::new(ld), Box::new(rd));

    let simplified_num = combined_num.simplify(env).ok()?;
    let simplified_den = combined_den.simplify(env).ok()?;

    if let Node::Num(ref n) = simplified_den {
        if n.is_one() {
            return Some(simplified_num);
        }
    }

    let result = Node::Divide(Box::new(simplified_num), Box::new(simplified_den));
    result.simplify(env).ok()
}

/// Collect all multiplicative factors from a nested Multiply tree.
fn collect_multiply_factors<'a>(node: &'a Node, factors: &mut Vec<&'a Node>) {
    if let Node::Multiply(left, right) = node {
        collect_multiply_factors(left, factors);
        collect_multiply_factors(right, factors);
    } else {
        factors.push(node);
    }
}

/// Try to simplify √(product) by extracting numeric square factors and even powers.
/// Returns Some(simplified) if any extraction succeeded, None otherwise.
fn try_simplify_sqrt_product(
    radicand: &Node,
    env: &Environment,
    use_node_sqrt: bool,
) -> Option<Node> {
    let mut factors = Vec::new();
    collect_multiply_factors(radicand, &mut factors);

    if factors.len() < 2 {
        return None;
    }

    let mut outside: Vec<Node> = Vec::new();
    let mut inside: Vec<Node> = Vec::new();
    let mut changed = false;

    for factor in &factors {
        match factor {
            Node::Num(n) => {
                if let Some(val) = n.to_i64() {
                    if val > 0 {
                        let (out, inn) = extract_square_factors(val as u64);
                        if out > 1 {
                            outside.push(Node::Num(ExactNum::integer(out as i64)));
                            changed = true;
                        }
                        if inn > 1 {
                            inside.push(Node::Num(ExactNum::integer(inn as i64)));
                        }
                        continue;
                    }
                }
                inside.push((*factor).clone());
            }
            Node::Power(base, exp) => {
                if let Node::Num(ref e) = **exp {
                    if let Some(exp_val) = e.to_i64() {
                        if exp_val >= 2 && exp_val % 2 == 0 {
                            let half_exp = exp_val / 2;
                            let base_factor = if half_exp == 1 {
                                if let Node::Variable(ref v) = **base {
                                    if env.assumptions().is_nonneg(v) {
                                        *base.clone()
                                    } else {
                                        Node::Abs(base.clone())
                                    }
                                } else {
                                    Node::Abs(base.clone())
                                }
                            } else {
                                Node::Power(
                                    base.clone(),
                                    Box::new(Node::Num(ExactNum::integer(half_exp))),
                                )
                            };
                            outside.push(base_factor);
                            changed = true;
                            continue;
                        }
                        if exp_val > 2 && exp_val % 2 == 1 {
                            let half = exp_val / 2;
                            if let Node::Variable(ref v) = **base {
                                if env.assumptions().is_nonneg(v) {
                                    if half == 1 {
                                        outside.push(*base.clone());
                                    } else {
                                        outside.push(Node::Power(
                                            base.clone(),
                                            Box::new(Node::Num(ExactNum::integer(half))),
                                        ));
                                    }
                                } else {
                                    let abs_base = Node::Abs(base.clone());
                                    if half == 1 {
                                        outside.push(abs_base);
                                    } else {
                                        outside.push(Node::Power(
                                            Box::new(abs_base),
                                            Box::new(Node::Num(ExactNum::integer(half))),
                                        ));
                                    }
                                }
                                inside.push(*base.clone());
                                changed = true;
                                continue;
                            }
                        }
                    }
                }
                inside.push((*factor).clone());
            }
            _ => {
                inside.push((*factor).clone());
            }
        }
    }

    if !changed {
        return None;
    }

    let outside_node = if outside.is_empty() {
        None
    } else {
        Some(
            outside
                .into_iter()
                .reduce(|a, b| Node::Multiply(Box::new(a), Box::new(b)))
                .unwrap(),
        )
    };

    let inside_node = if inside.is_empty() {
        None
    } else {
        let inner = inside
            .into_iter()
            .reduce(|a, b| Node::Multiply(Box::new(a), Box::new(b)))
            .unwrap();
        let sqrt_node = if use_node_sqrt {
            Node::Sqrt(Box::new(inner))
        } else {
            Node::Function("sqrt".to_string(), vec![inner])
        };
        Some(sqrt_node)
    };

    let result = match (outside_node, inside_node) {
        (Some(out), Some(sqrt)) => Node::Multiply(Box::new(out), Box::new(sqrt)),
        (Some(out), None) => out,
        (None, Some(sqrt)) => sqrt,
        (None, None) => Node::Num(ExactNum::one()),
    };

    let simplified = crate::simplify::Simplifiable::simplify(&result, env).unwrap_or(result);
    Some(simplified)
}

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

    #[test]
    fn test_extract_square_factors() {
        assert_eq!(extract_square_factors(12), (2, 3));
        assert_eq!(extract_square_factors(8), (2, 2));
        assert_eq!(extract_square_factors(18), (3, 2));
        assert_eq!(extract_square_factors(7), (1, 7));
        assert_eq!(extract_square_factors(4), (2, 1));
        assert_eq!(extract_square_factors(1), (1, 1));
        assert_eq!(extract_square_factors(72), (6, 2));
        assert_eq!(extract_square_factors(100), (10, 1));
        assert_eq!(extract_square_factors(50), (5, 2));
        assert_eq!(extract_square_factors(0), (0, 0));
    }
}
