`ca03027`

Extend CAS symbolic sums/products and simplification

Authored by mfwolffe <wolffemf@dukes.jmu.edu> 2 months ago

SHA: ca03027b84896b0fb0d1be860a26cb8d244f25b3
Parents: 61a2065
Tree: 65de631

5 changed files

Status	File	+	-
M	`garcalc-cas/src/eval.rs`	843	54
M	`garcalc-cas/src/expr.rs`	51	3
M	`garcalc-cas/src/lib.rs`	5	5
M	`garcalc-cas/src/parser.rs`	5	5
M	`garcalc-cas/src/symbolic.rs`	177	146

garcalc-cas/src/eval.rsmodified

1085 lines changed — click to load

  //!
  //! Supports both numeric and symbolic evaluation modes.
 +use std::collections::BTreeSet;
  use std::collections::HashMap;
  use std::f64::consts::{E, PI};
              Expr::Func(name, args) => {
                  // Symbolic functions operate on unevaluated expressions
                  match name.as_str() {
 -                    "diff" | "derivative" | "integrate" | "integral" |
 -                    "limit" | "lim" |
 -                    "solve" | "simplify" | "expand" | "factor" |
 -                    "substitute" | "subs" => {
 +                    "diff" | "derivative" | "integrate" | "integral" | "limit" | "lim"
 +                    | "solve" | "simplify" | "expand" | "factor" | "substitute" | "subs" => {
                          self.call_function(name, args)
+                     }
                      _ => {
+             }
              // Symbolic operations - perform symbolic computation then try to evaluate
 -            Expr::Derivative { expr: inner, var, order } => {
 +            Expr::Derivative {
 +                expr: inner,
 +                var,
 +                order,
 +            } => {
                  let result = Differentiator::diff_n(inner, var, *order)?;
                  let simplified = Simplifier::simplify(&result);
                  // Try to evaluate the result
+                 }
+             }
 -            Expr::Integral { expr: inner, var, lower, upper } => {
 +            Expr::Integral {
 +                expr: inner,
 +                var,
 +                lower,
 +                upper,
 +            } => {
                  if let (Some(l), Some(u)) = (lower, upper) {
                      // Definite integral - evaluate to number
                      let result = Integrator::integrate_definite(inner, var, l, u)?;
+                 }
+             }
 -            Expr::Limit { expr: inner, var, point, direction } => {
 +            Expr::Limit {
 +                expr: inner,
 +                var,
 +                point,
 +                direction,
 +            } => {
                  let result = Limits::limit(inner, var, point, *direction)?;
                  let simplified = Simplifier::simplify(&result);
                  self.eval(&simplified)
+             }
 -            Expr::Sum { .. }
 -            | Expr::Product { .. } => {
 -                if self.exact_mode {
 -                    Ok(expr.clone())
 -                } else {
 -                    Err(CasError::NotImplemented(
 -                        "sums/products not yet implemented".to_string(),
 -                    ))
 -                }
 -            }
 +            Expr::Sum {
 +                expr: inner,
 +                var,
 +                lower,
 +                upper,
 +            } => self.eval_sum(inner, var, lower, upper),
++
 +            Expr::Product {
 +                expr: inner,
 +                var,
 +                lower,
 +                upper,
 +            } => self.eval_product(inner, var, lower, upper),
              Expr::Inequality { lhs, op, rhs } => {
                  let lhs_val = self.eval(lhs)?;
                  Ok(Expr::Complex(re + x, *im))
+             }
              (Expr::Complex(re, im), Expr::Integer(n))
 -            | (Expr::Integer(n), Expr::Complex(re, im)) => {
 -                Ok(Expr::Complex(re + *n as f64, *im))
 -            }
 +            | (Expr::Integer(n), Expr::Complex(re, im)) => Ok(Expr::Complex(re + *n as f64, *im)),
              (Expr::Rational(r1), Expr::Rational(r2)) => {
                  let num = r1.num * r2.den + r2.num * r1.den;
                  let den = r1.den * r2.den;
+     }
      fn multiply(&self, a: &Expr, b: &Expr) -> Result<Expr> {
 +        if a.is_one() {
 +            return Ok(b.clone());
 +        }
 +        if b.is_one() {
 +            return Ok(a.clone());
 +        }
 +        if a.is_negative_one() {
 +            return self.negate(b);
 +        }
 +        if b.is_negative_one() {
 +            return self.negate(a);
 +        }
++
          match (a, b) {
              (Expr::Integer(x), Expr::Integer(y)) => Ok(Expr::Integer(x * y)),
              (Expr::Float(x), Expr::Float(y)) => Ok(Expr::Float(x * y)),
                  let n = *n as f64;
                  Ok(Expr::Complex(re * n, im * n))
+             }
 -            (Expr::Rational(r1), Expr::Rational(r2)) => {
 -                Ok(Expr::Rational(Rational::new(r1.num * r2.num, r1.den * r2.den)))
 -            }
 +            (Expr::Rational(r1), Expr::Rational(r2)) => Ok(Expr::Rational(Rational::new(
 +                r1.num * r2.num,
 +                r1.den * r2.den,
 +            ))),
              (Expr::Rational(r), Expr::Integer(n)) | (Expr::Integer(n), Expr::Rational(r)) => {
                  Ok(Expr::Rational(Rational::new(r.num * n, r.den)))
+             }
+             }
              ("min", _, _) if !args.is_empty() => {
                  let values: Result<Vec<f64>> = args.iter().map(|a| self.to_f64(a)).collect();
 -                let min = values?
 -                    .into_iter()
 -                    .fold(f64::INFINITY, |a, b| a.min(b));
 +                let min = values?.into_iter().fold(f64::INFINITY, |a, b| a.min(b));
                  Ok(Expr::Float(min))
+             }
              ("max", _, _) if !args.is_empty() => {
                  let values: Result<Vec<f64>> = args.iter().map(|a| self.to_f64(a)).collect();
 -                let max = values?
 -                    .into_iter()
 -                    .fold(f64::NEG_INFINITY, |a, b| a.max(b));
 +                let max = values?.into_iter().fold(f64::NEG_INFINITY, |a, b| a.max(b));
                  Ok(Expr::Float(max))
+             }
              ("gcd", 2, _) => {
                      let result = Differentiator::diff(&args[0], var)?;
                      Ok(Simplifier::simplify(&result))
                  } else {
 -                    Err(CasError::Type("diff requires variable as second argument".to_string()))
 +                    Err(CasError::Type(
 +                        "diff requires variable as second argument".to_string(),
 +                    ))
+                 }
+             }
              ("diff", 3, _) | ("derivative", 3, _) => {
                      let result = Differentiator::diff_n(&args[0], var, *n as u32)?;
                      Ok(Simplifier::simplify(&result))
                  } else {
 -                    Err(CasError::Type("diff requires variable and integer order".to_string()))
 +                    Err(CasError::Type(
 +                        "diff requires variable and integer order".to_string(),
 +                    ))
+                 }
+             }
                      let result = Integrator::integrate(&args[0], var)?;
                      Ok(Simplifier::simplify(&result))
                  } else {
 -                    Err(CasError::Type("integrate requires variable as second argument".to_string()))
 +                    Err(CasError::Type(
 +                        "integrate requires variable as second argument".to_string(),
 +                    ))
+                 }
+             }
              ("integrate", 4, _) | ("integral", 4, _) => {
                      // Try to evaluate the result numerically
                      self.eval(&Simplifier::simplify(&result))
                  } else {
 -                    Err(CasError::Type("integrate requires variable as second argument".to_string()))
 +                    Err(CasError::Type(
 +                        "integrate requires variable as second argument".to_string(),
 +                    ))
+                 }
+             }
                          Ok(Expr::Vector(solutions))
+                     }
                  } else {
 -                    Err(CasError::Type("solve requires variable as second argument".to_string()))
 +                    Err(CasError::Type(
 +                        "solve requires variable as second argument".to_string(),
 +                    ))
+                 }
+             }
 -            ("simplify", 1, _) => {
 -                Ok(Simplifier::simplify(&args[0]))
 +            ("sum", 4, _) => {
 +                if let Expr::Symbol(var) = &args[1] {
 +                    self.eval_sum(&args[0], var, &args[2], &args[3])
 +                } else {
 +                    Err(CasError::Type(
 +                        "sum requires variable as second argument".to_string(),
 +                    ))
 +                }
+             }
 -            ("expand", 1, _) => {
 -                Ok(Simplifier::simplify(&Simplifier::expand(&args[0])))
 +            ("product", 4, _) | ("prod", 4, _) => {
 +                if let Expr::Symbol(var) = &args[1] {
 +                    self.eval_product(&args[0], var, &args[2], &args[3])
 +                } else {
 +                    Err(CasError::Type(
 +                        "product requires variable as second argument".to_string(),
 +                    ))
 +                }
+             }
 +            ("simplify", 1, _) => Ok(Simplifier::simplify(&args[0])),
++
 +            ("expand", 1, _) => Ok(Simplifier::simplify(&Simplifier::expand(&args[0]))),
++
              ("factor", 1, _) => {
                  // Basic factoring - just return simplified for now
                  // Full factoring is complex, can add later
              ("substitute", 3, _) | ("subs", 3, _) => {
                  // substitute(expr, var, replacement)
                  if let Expr::Symbol(var) = &args[1] {
 -                    Ok(Simplifier::simplify(&Simplifier::substitute(&args[0], var, &args[2])))
 +                    Ok(Simplifier::simplify(&Simplifier::substitute(
 +                        &args[0], var, &args[2],
 +                    )))
                  } else {
 -                    Err(CasError::Type("substitute requires variable as second argument".to_string()))
 +                    Err(CasError::Type(
 +                        "substitute requires variable as second argument".to_string(),
 +                    ))
+                 }
+             }
                      let result = Limits::limit(&args[0], var, &args[2], None)?;
                      self.eval(&Simplifier::simplify(&result))
                  } else {
 -                    Err(CasError::Type("limit requires variable as second argument".to_string()))
 +                    Err(CasError::Type(
 +                        "limit requires variable as second argument".to_string(),
 +                    ))
+                 }
+             }
              ("limit", 4, _) | ("lim", 4, _) => {
                      let result = Limits::limit(&args[0], var, &args[2], direction)?;
                      self.eval(&Simplifier::simplify(&result))
                  } else {
 -                    Err(CasError::Type("limit requires variable as second argument".to_string()))
 +                    Err(CasError::Type(
 +                        "limit requires variable as second argument".to_string(),
 +                    ))
+                 }
+             }
                  if let Expr::Matrix(rows) = &args[0] {
                      self.matrix_transpose(rows)
                  } else {
 -                    Err(CasError::Type("transpose requires a matrix argument".to_string()))
 +                    Err(CasError::Type(
 +                        "transpose requires a matrix argument".to_string(),
 +                    ))
+                 }
+             }
                  if let Expr::Matrix(rows) = &args[0] {
                      self.matrix_trace(rows)
                  } else {
 -                    Err(CasError::Type("trace requires a matrix argument".to_string()))
 +                    Err(CasError::Type(
 +                        "trace requires a matrix argument".to_string(),
 +                    ))
+                 }
+             }
                  if let (Expr::Matrix(a), Expr::Matrix(b)) = (&args[0], &args[1]) {
                      self.matrix_mul(a, b)
                  } else {
 -                    Err(CasError::Type("matmul requires two matrix arguments".to_string()))
 +                    Err(CasError::Type(
 +                        "matmul requires two matrix arguments".to_string(),
 +                    ))
+                 }
+             }
              ("identity", 1, Some(n)) => {
                  let n = n as usize;
                  if n == 0 || n > 100 {
 -                    return Err(CasError::EvaluationError("identity matrix size must be 1-100".to_string()));
 +                    return Err(CasError::EvaluationError(
 +                        "identity matrix size must be 1-100".to_string(),
 +                    ));
+                 }
                  let mut rows = Vec::with_capacity(n);
                  for i in 0..n {
              return Err(CasError::EvaluationError("empty matrix".to_string()));
+         }
          if rows.iter().any(|r| r.len() != n) {
 -            return Err(CasError::EvaluationError("det requires square matrix".to_string()));
 +            return Err(CasError::EvaluationError(
 +                "det requires square matrix".to_string(),
 +            ));
+         }
          // Convert to f64 for numerical computation
              return Err(CasError::EvaluationError("empty matrix".to_string()));
+         }
          if rows.iter().any(|r| r.len() != n) {
 -            return Err(CasError::EvaluationError("inv requires square matrix".to_string()));
 +            return Err(CasError::EvaluationError(
 +                "inv requires square matrix".to_string(),
 +            ));
+         }
          // Convert to f64
              return Err(CasError::EvaluationError("empty matrix".to_string()));
+         }
          if rows.iter().any(|r| r.len() != n) {
 -            return Err(CasError::EvaluationError("trace requires square matrix".to_string()));
 +            return Err(CasError::EvaluationError(
 +                "trace requires square matrix".to_string(),
 +            ));
+         }
          let mut sum = 0.0;
          if b.len() != n {
              return Err(CasError::EvaluationError(format!(
                  "matrix dimensions don't match for multiplication: {}x{} * {}x{}",
 -                m, n, b.len(), p
 +                m,
 +                n,
 +                b.len(),
 +                p
              )));
+         }
          // Convert to f64
 -        let a_num: Vec<Vec<f64>> = a.iter()
 -            .map(|row| row.iter().map(|e| self.to_f64(e)).collect::<Result<Vec<_>>>())
 +        let a_num: Vec<Vec<f64>> = a
 +            .iter()
 +            .map(|row| {
 +                row.iter()
 +                    .map(|e| self.to_f64(e))
 +                    .collect::<Result<Vec<_>>>()
 +            })
              .collect::<Result<Vec<_>>>()?;
 -        let b_num: Vec<Vec<f64>> = b.iter()
 -            .map(|row| row.iter().map(|e| self.to_f64(e)).collect::<Result<Vec<_>>>())
 +        let b_num: Vec<Vec<f64>> = b
 +            .iter()
 +            .map(|row| {
 +                row.iter()
 +                    .map(|e| self.to_f64(e))
 +                    .collect::<Result<Vec<_>>>()
 +            })
              .collect::<Result<Vec<_>>>()?;
          let mut result = Vec::with_capacity(m);
          Ok(Expr::Matrix(result))
+     }
++
 +    fn eval_integer_bound(&self, bound: &Expr) -> Result<i64> {
 +        let value = self.eval(bound)?;
 +        match value {
 +            Expr::Integer(n) => Ok(n),
 +            Expr::Rational(r) if r.den == 1 => Ok(r.num),
 +            Expr::Float(x) => {
 +                if !x.is_finite() {
 +                    return Err(CasError::Type("bound must be a finite number".to_string()));
 +                }
 +                let rounded = x.round();
 +                if (x - rounded).abs() < 1e-10
 +                    && rounded >= i64::MIN as f64
 +                    && rounded <= i64::MAX as f64
 +                {
 +                    Ok(rounded as i64)
 +                } else {
 +                    Err(CasError::Type(format!(
 +                        "bound must be an integer, got {}",
 +                        Expr::Float(x)
 +                    )))
 +                }
 +            }
 +            _ => Err(CasError::Type(format!(
 +                "bound must be an integer, got {value}"
 +            ))),
 +        }
 +    }
++
 +    fn eval_sum(
 +        &self,
 +        body: &Expr,
 +        var: &crate::expr::Symbol,
 +        lower: &Expr,
 +        upper: &Expr,
 +    ) -> Result<Expr> {
 +        let inferred_var;
 +        let active_var = if body.contains_var(var) {
 +            var
 +        } else {
 +            inferred_var = Self::infer_iteration_var(body);
 +            inferred_var.as_ref().unwrap_or(var)
 +        };
++
 +        match self.eval_discrete_series(body, active_var, lower, upper, false) {
 +            Ok(value) => Ok(value),
 +            Err(_err) => {
 +                if let Some(symbolic) = Self::symbolic_sum(body, active_var, lower, upper) {
 +                    return Ok(Simplifier::simplify(&symbolic));
 +                }
++
 +                Ok(Expr::Sum {
 +                    expr: Box::new(body.clone()),
 +                    var: active_var.clone(),
 +                    lower: Box::new(lower.clone()),
 +                    upper: Box::new(upper.clone()),
 +                })
 +            }
 +        }
 +    }
++
 +    fn eval_product(
 +        &self,
 +        body: &Expr,
 +        var: &crate::expr::Symbol,
 +        lower: &Expr,
 +        upper: &Expr,
 +    ) -> Result<Expr> {
 +        let inferred_var;
 +        let active_var = if body.contains_var(var) {
 +            var
 +        } else {
 +            inferred_var = Self::infer_iteration_var(body);
 +            inferred_var.as_ref().unwrap_or(var)
 +        };
++
 +        match self.eval_discrete_series(body, active_var, lower, upper, true) {
 +            Ok(value) => Ok(value),
 +            Err(_err) => {
 +                if let Some(symbolic) = Self::symbolic_product(body, active_var, lower, upper) {
 +                    return Ok(Simplifier::simplify(&symbolic));
 +                }
++
 +                Ok(Expr::Product {
 +                    expr: Box::new(body.clone()),
 +                    var: active_var.clone(),
 +                    lower: Box::new(lower.clone()),
 +                    upper: Box::new(upper.clone()),
 +                })
 +            }
 +        }
 +    }
++
 +    fn infer_iteration_var(body: &Expr) -> Option<crate::expr::Symbol> {
 +        let mut vars = BTreeSet::new();
 +        Self::collect_symbols(body, &mut vars);
 +        if vars.len() == 1 {
 +            vars.into_iter().next().map(crate::expr::Symbol::new)
 +        } else {
 +            None
 +        }
 +    }
++
 +    fn collect_symbols(expr: &Expr, out: &mut BTreeSet<String>) {
 +        match expr {
 +            Expr::Symbol(s) => {
 +                out.insert(s.as_str().to_string());
 +            }
 +            Expr::Neg(inner) => Self::collect_symbols(inner, out),
 +            Expr::Add(terms) | Expr::Mul(terms) | Expr::Vector(terms) => {
 +                for term in terms {
 +                    Self::collect_symbols(term, out);
 +                }
 +            }
 +            Expr::Pow(base, exp) => {
 +                Self::collect_symbols(base, out);
 +                Self::collect_symbols(exp, out);
 +            }
 +            Expr::Func(_, args) => {
 +                for arg in args {
 +                    Self::collect_symbols(arg, out);
 +                }
 +            }
 +            Expr::Derivative { expr, .. } => Self::collect_symbols(expr, out),
 +            Expr::Integral {
 +                expr, lower, upper, ..
 +            } => {
 +                Self::collect_symbols(expr, out);
 +                if let Some(lo) = lower {
 +                    Self::collect_symbols(lo, out);
 +                }
 +                if let Some(hi) = upper {
 +                    Self::collect_symbols(hi, out);
 +                }
 +            }
 +            Expr::Limit { expr, point, .. } => {
 +                Self::collect_symbols(expr, out);
 +                Self::collect_symbols(point, out);
 +            }
 +            Expr::Sum {
 +                expr, lower, upper, ..
 +            }
 +            | Expr::Product {
 +                expr, lower, upper, ..
 +            } => {
 +                Self::collect_symbols(expr, out);
 +                Self::collect_symbols(lower, out);
 +                Self::collect_symbols(upper, out);
 +            }
 +            Expr::Equation(lhs, rhs) => {
 +                Self::collect_symbols(lhs, out);
 +                Self::collect_symbols(rhs, out);
 +            }
 +            Expr::Inequality { lhs, rhs, .. } => {
 +                Self::collect_symbols(lhs, out);
 +                Self::collect_symbols(rhs, out);
 +            }
 +            Expr::Matrix(rows) => {
 +                for row in rows {
 +                    for elem in row {
 +                        Self::collect_symbols(elem, out);
 +                    }
 +                }
 +            }
 +            Expr::Integer(_)
 +            | Expr::Rational(_)
 +            | Expr::Float(_)
 +            | Expr::Complex(_, _)
 +            | Expr::Undefined
 +            | Expr::Infinity(_) => {}
 +        }
 +    }
++
 +    fn symbolic_sum(
 +        body: &Expr,
 +        var: &crate::expr::Symbol,
 +        lower: &Expr,
 +        upper: &Expr,
 +    ) -> Option<Expr> {
 +        // Finite count for symbolic bounds: upper - lower + 1.
 +        let count = Expr::add(vec![
 +            Expr::sub(upper.clone(), lower.clone()),
 +            Expr::Integer(1),
 +        ]);
++
 +        if !body.contains_var(var) {
 +            return Some(Expr::mul(vec![body.clone(), count]));
 +        }
++
 +        match body {
 +            Expr::Symbol(s) if s == var => Some(Self::sum_linear(lower, upper)),
++
 +            Expr::Pow(base, exp) if matches!(base.as_ref(), Expr::Symbol(s) if s == var) => {
 +                match exp.as_ref() {
 +                    Expr::Integer(0) => Some(count),
 +                    Expr::Integer(1) => Some(Self::sum_linear(lower, upper)),
 +                    Expr::Integer(2) => Some(Self::sum_square(lower, upper)),
 +                    Expr::Integer(3) => Some(Self::sum_cube(lower, upper)),
 +                    _ => None,
 +                }
 +            }
++
 +            Expr::Neg(inner) => Self::symbolic_sum(inner, var, lower, upper).map(Expr::neg),
++
 +            Expr::Add(terms) => {
 +                let mut summed_terms = Vec::with_capacity(terms.len());
 +                for term in terms {
 +                    summed_terms.push(Self::symbolic_sum(term, var, lower, upper)?);
 +                }
 +                Some(Expr::add(summed_terms))
 +            }
++
 +            Expr::Mul(factors) => {
 +                let (mut independent, dependent): (Vec<Expr>, Vec<Expr>) =
 +                    factors.iter().cloned().partition(|f| !f.contains_var(var));
++
 +                if dependent.is_empty() {
 +                    independent.push(count);
 +                    return Some(Expr::mul(independent));
 +                }
++
 +                if dependent.len() == 1 {
 +                    let dep_sum = Self::symbolic_sum(&dependent[0], var, lower, upper)?;
 +                    independent.push(dep_sum);
 +                    return Some(Expr::mul(independent));
 +                }
++
 +                None
 +            }
++
 +            _ => None,
 +        }
 +    }
++
 +    fn sum_linear(lower: &Expr, upper: &Expr) -> Expr {
 +        fn triangular(x: Expr) -> Expr {
 +            Expr::mul(vec![
 +                x.clone(),
 +                Expr::add(vec![x, Expr::Integer(1)]),
 +                Expr::Rational(Rational::new(1, 2)),
 +            ])
 +        }
++
 +        let lo_minus_one = Expr::sub(lower.clone(), Expr::Integer(1));
 +        Expr::sub(triangular(upper.clone()), triangular(lo_minus_one))
 +    }
++
 +    fn sum_square(lower: &Expr, upper: &Expr) -> Expr {
 +        fn square_sum_prefix(x: Expr) -> Expr {
 +            let two_x_plus_one = Expr::add(vec![
 +                Expr::mul(vec![Expr::Integer(2), x.clone()]),
 +                Expr::Integer(1),
 +            ]);
 +            Expr::mul(vec![
 +                x.clone(),
 +                Expr::add(vec![x, Expr::Integer(1)]),
 +                two_x_plus_one,
 +                Expr::Rational(Rational::new(1, 6)),
 +            ])
 +        }
++
 +        let lo_minus_one = Expr::sub(lower.clone(), Expr::Integer(1));
 +        Expr::sub(
 +            square_sum_prefix(upper.clone()),
 +            square_sum_prefix(lo_minus_one),
 +        )
 +    }
++
 +    fn sum_cube(lower: &Expr, upper: &Expr) -> Expr {
 +        fn cube_sum_prefix(x: Expr) -> Expr {
 +            let tri = Expr::mul(vec![
 +                x.clone(),
 +                Expr::add(vec![x, Expr::Integer(1)]),
 +                Expr::Rational(Rational::new(1, 2)),
 +            ]);
 +            Expr::pow(tri, Expr::Integer(2))
 +        }
++
 +        let lo_minus_one = Expr::sub(lower.clone(), Expr::Integer(1));
 +        Expr::sub(
 +            cube_sum_prefix(upper.clone()),
 +            cube_sum_prefix(lo_minus_one),
 +        )
 +    }
++
 +    fn symbolic_product(
 +        body: &Expr,
 +        var: &crate::expr::Symbol,
 +        lower: &Expr,
 +        upper: &Expr,
 +    ) -> Option<Expr> {
 +        let count = Expr::add(vec![
 +            Expr::sub(upper.clone(), lower.clone()),
 +            Expr::Integer(1),
 +        ]);
++
 +        if !body.contains_var(var) {
 +            return Some(Expr::pow(body.clone(), count));
 +        }
++
 +        if let Some(factored) = Self::factor_simple_product_body(body, var) {
 +            return Self::symbolic_product(&factored, var, lower, upper);
 +        }
++
 +        match body {
 +            Expr::Symbol(s) if s == var => Some(Self::factorial_range(lower, upper)),
 +            Expr::Neg(inner) => {
 +                let inner_product = Self::symbolic_product(inner, var, lower, upper)?;
 +                Some(Expr::mul(vec![
 +                    Expr::pow(Expr::Integer(-1), count),
 +                    inner_product,
 +                ]))
 +            }
 +            Expr::Mul(factors) => {
 +                let (independent, dependent): (Vec<Expr>, Vec<Expr>) =
 +                    factors.iter().cloned().partition(|f| !f.contains_var(var));
++
 +                let mut parts = Vec::new();
 +                if !independent.is_empty() {
 +                    parts.push(Expr::pow(Expr::mul(independent), count.clone()));
 +                }
++
 +                for dep in dependent {
 +                    parts.push(Self::symbolic_product(&dep, var, lower, upper)?);
 +                }
++
 +                Some(Expr::mul(parts))
 +            }
 +            Expr::Add(_) => Self::product_linear_term(body, var, lower, upper),
 +            Expr::Pow(base, exp) => {
 +                let Expr::Integer(power) = exp.as_ref() else {
 +                    return None;
 +                };
++
 +                if matches!(base.as_ref(), Expr::Symbol(s) if s == var) {
 +                    return Some(Expr::pow(
 +                        Self::factorial_range(lower, upper),
 +                        Expr::Integer(*power),
 +                    ));
 +                }
++
 +                if let Some(linear_base_product) =
 +                    Self::product_linear_term(base, var, lower, upper)
 +                {
 +                    return Some(Expr::pow(linear_base_product, Expr::Integer(*power)));
 +                }
++
 +                None
 +            }
 +            _ => None,
 +        }
 +    }
++
 +    fn factorial_range(lower: &Expr, upper: &Expr) -> Expr {
 +        Self::factorial_range_shifted(lower, upper, 0).unwrap_or_else(|| {
 +            let upper_fact = Expr::func("factorial", vec![upper.clone()]);
 +            if matches!(lower, Expr::Integer(1)) {
 +                upper_fact
 +            } else {
 +                let lower_minus_one = Expr::sub(lower.clone(), Expr::Integer(1));
 +                let lower_fact = Expr::func("factorial", vec![lower_minus_one]);
 +                Expr::mul(vec![upper_fact, Expr::pow(lower_fact, Expr::Integer(-1))])
 +            }
 +        })
 +    }
++
 +    fn product_linear_term(
 +        expr: &Expr,
 +        var: &crate::expr::Symbol,
 +        lower: &Expr,
 +        upper: &Expr,
 +    ) -> Option<Expr> {
 +        let shift = Self::extract_linear_shift(expr, var)?;
 +        Self::factorial_range_shifted(lower, upper, shift)
 +    }
++
 +    fn extract_linear_shift(expr: &Expr, var: &crate::expr::Symbol) -> Option<i64> {
 +        match expr {
 +            Expr::Symbol(s) if s == var => Some(0),
 +            Expr::Add(terms) => {
 +                let mut saw_var = false;
 +                let mut shift = 0_i64;
 +                for term in terms {
 +                    match term {
 +                        Expr::Symbol(s) if s == var => {
 +                            if saw_var {
 +                                return None;
 +                            }
 +                            saw_var = true;
 +                        }
 +                        Expr::Integer(n) => {
 +                            shift = shift.checked_add(*n)?;
 +                        }
 +                        Expr::Neg(inner) => {
 +                            if let Expr::Integer(n) = inner.as_ref() {
 +                                shift = shift.checked_sub(*n)?;
 +                            } else {
 +                                return None;
 +                            }
 +                        }
 +                        _ => return None,
 +                    }
 +                }
 +                if saw_var {
 +                    Some(shift)
 +                } else {
 +                    None
 +                }
 +            }
 +            _ => None,
 +        }
 +    }
++
 +    fn factorial_range_shifted(lower: &Expr, upper: &Expr, shift: i64) -> Option<Expr> {
 +        let shift_minus_one = shift.checked_sub(1)?;
 +        if let Expr::Integer(lo) = lower {
 +            if lo.checked_add(shift_minus_one)? < 0 {
 +                return None;
 +            }
 +        }
 +        if let Expr::Integer(hi) = upper {
 +            if hi.checked_add(shift)? < 0 {
 +                return None;
 +            }
 +        }
++
 +        let shifted_upper =
 +            Simplifier::simplify(&Expr::add(vec![upper.clone(), Expr::Integer(shift)]));
 +        let shifted_lower_minus_one = Simplifier::simplify(&Expr::add(vec![
 +            lower.clone(),
 +            Expr::Integer(shift_minus_one),
 +        ]));
++
 +        let upper_fact = Expr::func("factorial", vec![shifted_upper]);
 +        if matches!(shifted_lower_minus_one, Expr::Integer(0)) {
 +            Some(upper_fact)
 +        } else {
 +            let lower_fact = Expr::func("factorial", vec![shifted_lower_minus_one]);
 +            Some(Expr::mul(vec![
 +                upper_fact,
 +                Expr::pow(lower_fact, Expr::Integer(-1)),
 +            ]))
 +        }
 +    }
++
 +    fn factor_simple_product_body(body: &Expr, var: &crate::expr::Symbol) -> Option<Expr> {
 +        let Expr::Add(terms) = body else {
 +            return None;
 +        };
 +        if terms.len() != 2 {
 +            return None;
 +        }
++
 +        let mut has_square = false;
 +        let mut linear_sign = 0_i64;
 +        for term in terms {
 +            match term {
 +                Expr::Pow(base, exp)
 +                    if matches!(base.as_ref(), Expr::Symbol(s) if s == var)
 +                        && matches!(exp.as_ref(), Expr::Integer(2)) =>
 +                {
 +                    has_square = true;
 +                }
 +                Expr::Symbol(s) if s == var => linear_sign += 1,
 +                Expr::Neg(inner) if matches!(inner.as_ref(), Expr::Symbol(s) if s == var) => {
 +                    linear_sign -= 1;
 +                }
 +                _ => return None,
 +            }
 +        }
++
 +        if !has_square {
 +            return None;
 +        }
++
 +        let var_expr = Expr::Symbol(var.clone());
 +        match linear_sign {
 +            1 => Some(Expr::mul(vec![
 +                var_expr.clone(),
 +                Expr::add(vec![var_expr, Expr::Integer(1)]),
 +            ])),
 +            -1 => Some(Expr::mul(vec![
 +                var_expr.clone(),
 +                Expr::add(vec![var_expr, Expr::Integer(-1)]),
 +            ])),
 +            _ => None,
 +        }
 +    }
++
 +    fn eval_discrete_series(
 +        &self,
 +        body: &Expr,
 +        var: &crate::expr::Symbol,
 +        lower: &Expr,
 +        upper: &Expr,
 +        is_product: bool,
 +    ) -> Result<Expr> {
 +        const MAX_TERMS: i64 = 100_000;
++
 +        let lo = self.eval_integer_bound(lower)?;
 +        let hi = self.eval_integer_bound(upper)?;
++
 +        if lo > hi {
 +            return Ok(if is_product {
 +                Expr::Integer(1)
 +            } else {
 +                Expr::Integer(0)
 +            });
 +        }
++
 +        let count = hi.saturating_sub(lo).saturating_add(1);
 +        if count > MAX_TERMS {
 +            let kind = if is_product { "product" } else { "sum" };
 +            return Err(CasError::EvaluationError(format!(
 +                "{kind} has too many terms ({count}); limit is {MAX_TERMS}"
 +            )));
 +        }
++
 +        let mut acc = if is_product {
 +            Expr::Integer(1)
 +        } else {
 +            Expr::Integer(0)
 +        };
++
 +        for n in lo..=hi {
 +            let substituted = Simplifier::substitute(body, var, &Expr::Integer(n));
 +            let term = self.eval(&substituted)?;
 +            acc = if is_product {
 +                self.multiply(&acc, &term)?
 +            } else {
 +                self.add(&acc, &term)?
 +            };
 +        }
++
 +        Ok(Simplifier::simplify(&acc))
 +    }
+ }
  // Helper functions
+         }
+     }
 +    fn expr_to_f64(expr: &Expr) -> f64 {
 +        match expr {
 +            Expr::Integer(n) => *n as f64,
 +            Expr::Rational(r) => r.to_f64(),
 +            Expr::Float(x) => *x,
 +            other => panic!("expected numeric expression, got {other}"),
 +        }
 +    }
++
      #[test]
      fn test_arithmetic() {
          assert_eq!(eval("2 + 3").unwrap(), Expr::Integer(5));
          let result = evaluator.eval(&expr).unwrap();
          assert_eq!(result, Expr::Integer(6));
+     }
++
 +    #[test]
 +    fn test_sum_evaluation() {
 +        assert_eq!(eval("sum(n, n, 1, 5)").unwrap(), Expr::Integer(15));
 +    }
++
 +    #[test]
 +    fn test_product_evaluation() {
 +        assert_eq!(eval("product(n, n, 1, 4)").unwrap(), Expr::Integer(24));
 +    }
++
 +    #[test]
 +    fn test_empty_discrete_range() {
 +        assert_eq!(eval("sum(n, n, 5, 1)").unwrap(), Expr::Integer(0));
 +        assert_eq!(eval("product(n, n, 5, 1)").unwrap(), Expr::Integer(1));
 +    }
++
 +    #[test]
 +    fn test_symbolic_solve_function() {
 +        let result = eval("solve(x^2 - 4, x)").unwrap();
 +        if let Expr::Vector(solutions) = result {
 +            assert_eq!(solutions.len(), 2);
 +            let values: Vec<f64> = solutions
 +                .iter()
 +                .map(|s| match s {
 +                    Expr::Integer(n) => *n as f64,
 +                    Expr::Rational(r) => r.to_f64(),
 +                    Expr::Float(x) => *x,
 +                    other => panic!("expected numeric solution, got {other}"),
 +                })
 +                .collect();
 +            assert!(values.iter().any(|v| (*v - 2.0).abs() < 1e-10));
 +            assert!(values.iter().any(|v| (*v + 2.0).abs() < 1e-10));
 +        } else {
 +            panic!("expected vector of solutions");
 +        }
 +    }
++
 +    #[test]
 +    fn test_symbolic_sum_linear_closed_form() {
 +        let symbolic = eval("sum(k, k, 1, n)").unwrap();
 +        let mut evaluator = Evaluator::new();
 +        evaluator.set_var("n", Expr::Integer(10));
 +        let value = evaluator.eval(&symbolic).unwrap();
 +        assert!((expr_to_f64(&value) - 55.0).abs() < 1e-10);
 +    }
++
 +    #[test]
 +    fn test_symbolic_sum_quadratic_closed_form() {
 +        let symbolic = eval("sum(k^2 + k, k, 1, n)").unwrap();
 +        let mut evaluator = Evaluator::new();
 +        evaluator.set_var("n", Expr::Integer(5));
 +        let value = evaluator.eval(&symbolic).unwrap();
 +        assert!((expr_to_f64(&value) - 70.0).abs() < 1e-10);
 +    }
++
 +    #[test]
 +    fn test_symbolic_sum_fallback_for_unreduced_form() {
 +        let symbolic = eval("sum(1/(k-1), k, 1, n)").unwrap();
 +        assert!(matches!(symbolic, Expr::Sum { .. }));
 +    }
++
 +    #[test]
 +    fn test_symbolic_product_constant_closed_form() {
 +        let symbolic = eval("product(2, k, 1, n)").unwrap();
 +        let mut evaluator = Evaluator::new();
 +        evaluator.set_var("n", Expr::Integer(5));
 +        let value = evaluator.eval(&symbolic).unwrap();
 +        assert!((expr_to_f64(&value) - 32.0).abs() < 1e-10);
 +    }
++
 +    #[test]
 +    fn test_symbolic_product_factorial_closed_form() {
 +        let symbolic = eval("product(k, k, 1, n)").unwrap();
 +        let mut evaluator = Evaluator::new();
 +        evaluator.set_var("n", Expr::Integer(5));
 +        let value = evaluator.eval(&symbolic).unwrap();
 +        assert!((expr_to_f64(&value) - 120.0).abs() < 1e-10);
 +    }
++
 +    #[test]
 +    fn test_symbolic_product_linear_shift_closed_form() {
 +        let symbolic = eval("product(k+1, k, 1, n)").unwrap();
 +        let mut evaluator = Evaluator::new();
 +        evaluator.set_var("n", Expr::Integer(5));
 +        let value = evaluator.eval(&symbolic).unwrap();
 +        assert!((expr_to_f64(&value) - 720.0).abs() < 1e-10);
 +    }
++
 +    #[test]
 +    fn test_symbolic_product_mul_decomposition_closed_form() {
 +        let symbolic = eval("product(2*k, k, 1, n)").unwrap();
 +        let mut evaluator = Evaluator::new();
 +        evaluator.set_var("n", Expr::Integer(5));
 +        let value = evaluator.eval(&symbolic).unwrap();
 +        assert!((expr_to_f64(&value) - 3840.0).abs() < 1e-10);
 +    }
++
 +    #[test]
 +    fn test_symbolic_product_simple_quadratic_factoring() {
 +        let symbolic = eval("product(k^2 + k, k, 1, n)").unwrap();
 +        let mut evaluator = Evaluator::new();
 +        evaluator.set_var("n", Expr::Integer(4));
 +        let value = evaluator.eval(&symbolic).unwrap();
 +        assert!((expr_to_f64(&value) - 2880.0).abs() < 1e-10);
 +    }
++
 +    #[test]
 +    fn test_symbolic_product_telescoping_ratio() {
 +        let symbolic = eval("product(k/(k-1), k, 2, n)").unwrap();
 +        let mut evaluator = Evaluator::new();
 +        evaluator.set_var("n", Expr::Integer(6));
 +        let value = evaluator.eval(&symbolic).unwrap();
 +        assert!((expr_to_f64(&value) - 6.0).abs() < 1e-10);
 +    }
++
 +    #[test]
 +    fn test_symbolic_product_fallback_for_unreduced_form() {
 +        let symbolic = eval("product(k^2 + 2, k, 1, n)").unwrap();
 +        assert!(matches!(symbolic, Expr::Product { .. }));
 +    }
++
 +    #[test]
 +    fn test_sum_infers_iteration_var_when_template_var_unused() {
 +        let expr = Expr::Sum {
 +            expr: Box::new(Expr::Symbol(crate::expr::Symbol::new("n"))),
 +            var: crate::expr::Symbol::new("i"),
 +            lower: Box::new(Expr::Integer(1)),
 +            upper: Box::new(Expr::Integer(5)),
 +        };
 +        let result = Evaluator::new().eval(&expr).unwrap();
 +        assert_eq!(result, Expr::Integer(15));
 +    }
++
 +    #[test]
 +    fn test_product_infers_iteration_var_when_template_var_unused() {
 +        let expr = Expr::Product {
 +            expr: Box::new(Expr::Symbol(crate::expr::Symbol::new("n"))),
 +            var: crate::expr::Symbol::new("i"),
 +            lower: Box::new(Expr::Integer(1)),
 +            upper: Box::new(Expr::Integer(5)),
 +        };
 +        let result = Evaluator::new().eval(&expr).unwrap();
 +        assert_eq!(result, Expr::Integer(120));
 +    }
++
 +    #[test]
 +    fn test_exact_multiply_rational_one_identity() {
 +        let mut evaluator = Evaluator::new();
 +        evaluator.exact_mode = true;
++
 +        let expr = Expr::mul(vec![
 +            Expr::Rational(crate::expr::Rational::new(1, 1)),
 +            Expr::symbol("n"),
 +        ]);
 +        let result = evaluator.eval(&expr).unwrap();
++
 +        assert_eq!(result, Expr::symbol("n"));
 +    }
+ }

garcalc-cas/src/expr.rsmodified

      pub fn is_zero(&self) -> bool {
          match self {
              Self::Integer(0) => true,
 +            Self::Rational(r) if r.den != 0 && r.num == 0 => true,
              Self::Float(x) if *x == 0.0 => true,
              _ => false,
+         }
      pub fn is_one(&self) -> bool {
          match self {
              Self::Integer(1) => true,
 +            Self::Rational(r) if r.den != 0 && r.num == r.den => true,
              Self::Float(x) if *x == 1.0 => true,
              _ => false,
+         }
      pub fn is_negative_one(&self) -> bool {
          match self {
              Self::Integer(-1) => true,
 +            Self::Rational(r) if r.den != 0 && r.num == -r.den => true,
              Self::Float(x) if *x == -1.0 => true,
              _ => false,
+         }
              Self::Derivative { expr, .. } => expr.contains_var(var),
              Self::Integral { expr, .. } => expr.contains_var(var),
              Self::Limit { expr, point, .. } => expr.contains_var(var) || point.contains_var(var),
 -            Self::Sum { expr, lower, upper, .. } | Self::Product { expr, lower, upper, .. } => {
 -                expr.contains_var(var) || lower.contains_var(var) || upper.contains_var(var)
 +            Self::Sum {
 +                expr, lower, upper, ..
+             }
 +            | Self::Product {
 +                expr, lower, upper, ..
 +            } => expr.contains_var(var) || lower.contains_var(var) || upper.contains_var(var),
              Self::Equation(lhs, rhs) => lhs.contains_var(var) || rhs.contains_var(var),
              Self::Inequality { lhs, rhs, .. } => lhs.contains_var(var) || rhs.contains_var(var),
 -            Self::Matrix(rows) => rows.iter().any(|row| row.iter().any(|e| e.contains_var(var))),
 +            Self::Matrix(rows) => rows
 +                .iter()
 +                .any(|row| row.iter().any(|e| e.contains_var(var))),
              Self::Vector(elems) => elems.iter().any(|e| e.contains_var(var)),
              _ => false,
+         }
+             }
              Self::Pow(base, exp) => write!(f, "{base}^{exp}"),
              Self::Func(name, args) => {
 +                if name == "factorial" && args.len() == 1 {
 +                    let arg = &args[0];
 +                    if factorial_arg_needs_parens(arg) {
 +                        write!(f, "({arg})!")?;
 +                    } else {
 +                        write!(f, "{arg}!")?;
 +                    }
 +                    return Ok(());
 +                }
++
                  write!(f, "{name}(")?;
                  for (i, arg) in args.iter().enumerate() {
                      if i > 0 {
+     }
+ }
 +fn factorial_arg_needs_parens(arg: &Expr) -> bool {
 +    matches!(arg, Expr::Neg(_) | Expr::Pow(_, _))
 +}
++
  #[cfg(test)]
  mod tests {
      use super::*;
          assert_eq!(expr.to_string(), "(x+(2*y))");
+     }
 +    #[test]
 +    fn test_factorial_display() {
 +        let simple = Expr::func("factorial", vec![Expr::symbol("n")]);
 +        assert_eq!(simple.to_string(), "n!");
++
 +        let grouped = Expr::func(
 +            "factorial",
 +            vec![Expr::add(vec![Expr::symbol("n"), Expr::integer(1)])],
 +        );
 +        assert_eq!(grouped.to_string(), "(n+1)!");
++
 +        let negated = Expr::func("factorial", vec![Expr::neg(Expr::symbol("n"))]);
 +        assert_eq!(negated.to_string(), "(-n)!");
 +    }
++
      #[test]
      fn test_contains_var() {
          let expr = Expr::add(vec![Expr::symbol("x"), Expr::integer(1)]);
          assert!(expr.contains_var(&Symbol::new("x")));
          assert!(!expr.contains_var(&Symbol::new("y")));
+     }
++
 +    #[test]
 +    fn test_numeric_identity_predicates_include_rationals() {
 +        let zero = Expr::Rational(Rational::new(0, 5));
 +        let one = Expr::Rational(Rational::new(2, 2));
 +        let negative_one = Expr::Rational(Rational::new(-3, 3));
++
 +        assert!(zero.is_zero());
 +        assert!(one.is_one());
 +        assert!(negative_one.is_negative_one());
 +    }
+ }

garcalc-cas/src/lib.rsmodified

  //! - Equation solving
  //! - Limits and series expansions
 +pub mod error;
 +pub mod eval;
  pub mod expr;
  pub mod parser;
 -pub mod eval;
 -pub mod error;
  pub mod symbolic;
 -pub use expr::{Expr, Symbol, Rational, LimitDirection, Sign};
 -pub use parser::Parser;
 -pub use eval::Evaluator;
  pub use error::{CasError, Result};
 +pub use eval::Evaluator;
 +pub use expr::{Expr, LimitDirection, Rational, Sign, Symbol};
 +pub use parser::Parser;
  pub use symbolic::{Differentiator, Integrator, Limits, Simplifier, Solver};

garcalc-cas/src/parser.rsmodified

                      if let Some(last) = factors.last() {
                          if matches!(
                              last,
 -                            Expr::Symbol(_)
 -                                | Expr::Integer(_)
 -                                | Expr::Float(_)
 -                                | Expr::Func(_, _)
 +                            Expr::Symbol(_) | Expr::Integer(_) | Expr::Float(_) | Expr::Func(_, _)
                          ) {
                              factors.push(self.parse_power()?);
                              continue;
                          _ => Symbol::new("x"),
                      };
                      let (lower, upper) = if args.len() >= 2 {
 -                        (Some(Box::new(args.remove(0))), Some(Box::new(args.remove(0))))
 +                        (
 +                            Some(Box::new(args.remove(0))),
 +                            Some(Box::new(args.remove(0))),
 +                        )
                      } else {
                          (None, None)
                      };

garcalc-cas/src/symbolic.rsmodified

  //!
  //! Provides symbolic differentiation, integration, simplification, and solving.
 -use crate::expr::{Expr, Symbol};
  use crate::error::{CasError, Result};
 +use crate::expr::{Expr, Symbol};
  /// Symbolic differentiator
  pub struct Differentiator;
              Expr::Func(name, args) => Self::diff_func(name, args, var),
              // Derivative of a derivative: just nest it
 -            Expr::Derivative { expr, var: v, order } => {
 +            Expr::Derivative {
 +                expr,
 +                var: v,
 +                order,
 +            } => {
                  if v == var {
                      Ok(Expr::Derivative {
                          expr: expr.clone(),
                  // 1 / (2 * sqrt(x))
                  Expr::div(
                      Expr::Integer(1),
 -                    Expr::mul(vec![Expr::Integer(2), Expr::func("sqrt", vec![arg.clone()])]),
 +                    Expr::mul(vec![
 +                        Expr::Integer(2),
 +                        Expr::func("sqrt", vec![arg.clone()]),
 +                    ]),
+                 )
+             }
              "cbrt" => {
      fn integrate_impl(expr: &Expr, var: &Symbol) -> Result<Expr> {
          // If expression doesn't contain variable, it's a constant
          if !expr.contains_var(var) {
 -            return Ok(Expr::mul(vec![
 -                expr.clone(),
 -                Expr::Symbol(var.clone()),
 -            ]));
 +            return Ok(Expr::mul(vec![expr.clone(), Expr::Symbol(var.clone())]));
+         }
          match expr {
                  if **base == Expr::Symbol(var.clone()) && !exp.contains_var(var) {
                      // Check for n = -1 case: ∫ x^(-1) dx = ln|x|
                      if exp.is_negative_one() {
 -                        return Ok(Expr::func("ln", vec![Expr::func(
 -                            "abs",
 -                            vec![Expr::Symbol(var.clone())],
 -                        )]));
 +                        return Ok(Expr::func(
 +                            "ln",
 +                            vec![Expr::func("abs", vec![Expr::Symbol(var.clone())])],
 +                        ));
+                     }
                      let n_plus_1 = Expr::add(vec![(**exp).clone(), Expr::Integer(1)]);
              // Trigonometric
              "sin" => Expr::neg(Expr::func("cos", vec![arg.clone()])),
              "cos" => Expr::func("sin", vec![arg.clone()]),
 -            "tan" => Expr::neg(Expr::func("ln", vec![Expr::func(
 -                "abs",
 -                vec![Expr::func("cos", vec![arg.clone()])],
 -            )])),
 -            "cot" => Expr::func("ln", vec![Expr::func(
 -                "abs",
 -                vec![Expr::func("sin", vec![arg.clone()])],
 -            )]),
 -            "sec" => Expr::func("ln", vec![Expr::func(
 -                "abs",
 -                vec![Expr::add(vec![
 -                    Expr::func("sec", vec![arg.clone()]),
 -                    Expr::func("tan", vec![arg.clone()]),
 -                ])],
 -            )]),
 -            "csc" => Expr::neg(Expr::func("ln", vec![Expr::func(
 -                "abs",
 -                vec![Expr::add(vec![
 -                    Expr::func("csc", vec![arg.clone()]),
 -                    Expr::func("cot", vec![arg.clone()]),
 -                ])],
 -            )])),
 +            "tan" => Expr::neg(Expr::func(
 +                "ln",
 +                vec![Expr::func(
 +                    "abs",
 +                    vec![Expr::func("cos", vec![arg.clone()])],
 +                )],
 +            )),
 +            "cot" => Expr::func(
 +                "ln",
 +                vec![Expr::func(
 +                    "abs",
 +                    vec![Expr::func("sin", vec![arg.clone()])],
 +                )],
 +            ),
 +            "sec" => Expr::func(
 +                "ln",
 +                vec![Expr::func(
 +                    "abs",
 +                    vec![Expr::add(vec![
 +                        Expr::func("sec", vec![arg.clone()]),
 +                        Expr::func("tan", vec![arg.clone()]),
 +                    ])],
 +                )],
 +            ),
 +            "csc" => Expr::neg(Expr::func(
 +                "ln",
 +                vec![Expr::func(
 +                    "abs",
 +                    vec![Expr::add(vec![
 +                        Expr::func("csc", vec![arg.clone()]),
 +                        Expr::func("cot", vec![arg.clone()]),
 +                    ])],
 +                )],
 +            )),
              // Hyperbolic
              "sinh" => Expr::func("cosh", vec![arg.clone()]),
  impl Limits {
      /// Compute the limit of an expression as var approaches point
 -    pub fn limit(expr: &Expr, var: &Symbol, point: &Expr, direction: Option<crate::expr::LimitDirection>) -> Result<Expr> {
 +    pub fn limit(
 +        expr: &Expr,
 +        var: &Symbol,
 +        point: &Expr,
 +        direction: Option<crate::expr::LimitDirection>,
 +    ) -> Result<Expr> {
          Self::limit_impl(expr, var, point, direction, 0)
+     }
                              // Handle x^(-n) where n > 1
                              if let Expr::Integer(n) = inner.as_ref() {
                                  if *n > 0 {
 -                                    denominator_parts.push(Expr::pow((**base).clone(), Expr::Integer(*n)));
 +                                    denominator_parts
 +                                        .push(Expr::pow((**base).clone(), Expr::Integer(*n)));
                                      continue;
+                                 }
+                             }
                      };
                      // Check if it's an indeterminate form (0/0 or ∞/∞)
 -                    let num_at_point = Simplifier::simplify(&Simplifier::substitute(&numerator, var, point));
 -                    let denom_at_point = Simplifier::simplify(&Simplifier::substitute(&denominator, var, point));
 +                    let num_at_point =
 +                        Simplifier::simplify(&Simplifier::substitute(&numerator, var, point));
 +                    let denom_at_point =
 +                        Simplifier::simplify(&Simplifier::substitute(&denominator, var, point));
                      let num_zero = Self::is_zero(&num_at_point);
                      let denom_zero = Self::is_zero(&denom_at_point);
                      let denom_inf = matches!(denom_at_point, Expr::Infinity(_));
                      if (num_zero && denom_zero) || (num_inf && denom_inf) {
 -                        return Self::try_lhopital(&numerator, &denominator, var, point, direction, depth);
 +                        return Self::try_lhopital(
 +                            &numerator,
 +                            &denominator,
 +                            var,
 +                            point,
 +                            direction,
 +                            depth,
 +                        );
+                     }
+                 }
                      // 1/0 -> infinity (sign depends on direction)
                      Ok(Expr::Infinity(crate::expr::Sign::Positive))
                  } else {
 -                    Ok(Simplifier::simplify(&Expr::pow(base_limit, Expr::Integer(-1))))
 +                    Ok(Simplifier::simplify(&Expr::pow(
 +                        base_limit,
 +                        Expr::Integer(-1),
 +                    )))
+                 }
+             }
              Expr::Func(name, args) if name == "exp" && args.len() == 1 => {
                  if args[0] == Expr::Symbol(var.clone()) {
                      match sign {
 -                        crate::expr::Sign::Positive => Ok(Expr::Infinity(crate::expr::Sign::Positive)),
 +                        crate::expr::Sign::Positive => {
 +                            Ok(Expr::Infinity(crate::expr::Sign::Positive))
 +                        }
                          crate::expr::Sign::Negative => Ok(Expr::Integer(0)),
+                     }
                  } else {
              let denom_deriv = Differentiator::diff(denom, var)?;
              // Recursive limit of f'/g'
 -            let quotient = Expr::mul(vec![
 -                num_deriv,
 -                Expr::pow(denom_deriv, Expr::Integer(-1)),
 -            ]);
 +            let quotient = Expr::mul(vec![num_deriv, Expr::pow(denom_deriv, Expr::Integer(-1))]);
              Self::limit_impl(&quotient, var, point, direction, depth + 1)
          } else if denom_is_zero && !num_is_zero {
              // Limit is ±∞ or doesn't exist
              // Check for 0 * ∞ or similar patterns in products
              Expr::Mul(factors) => {
                  let has_infinity = factors.iter().any(|f| {
 -                    matches!(f, Expr::Infinity(_))
 -                        || matches!(f, Expr::Rational(r) if r.den == 0)
 +                    matches!(f, Expr::Infinity(_)) || matches!(f, Expr::Rational(r) if r.den == 0)
                  });
                  let has_zero = factors.iter().any(|f| Self::is_zero(f));
                  // 0 * ∞ is indeterminate
+                 }
+             }
              Expr::Mul(factors) => factors.iter().map(|f| Self::degree_in(f, var)).sum(),
 -            Expr::Add(terms) => terms.iter().map(|t| Self::degree_in(t, var)).max().unwrap_or(0),
 +            Expr::Add(terms) => terms
 +                .iter()
 +                .map(|t| Self::degree_in(t, var))
 +                .max()
 +                .unwrap_or(0),
              Expr::Neg(e) => Self::degree_in(e, var),
              _ => 0,
+         }
              Expr::Func(name, args) => {
                  let sargs: Vec<_> = args.iter().map(Self::simplify_impl).collect();
 +                if name == "factorial" && sargs.len() == 1 {
 +                    if let Expr::Integer(n) = sargs[0] {
 +                        if n >= 0 {
 +                            if let Some(value) = Self::factorial_i64_checked(n) {
 +                                return Expr::Integer(value);
 +                            }
 +                        }
 +                    }
 +                }
++
                  // Try to evaluate sqrt of perfect squares
                  if name == "sqrt" && sargs.len() == 1 {
                      if let Expr::Integer(n) = &sargs[0] {
                      final_factors.insert(0, Expr::Integer(numerator));
+                 }
              } else {
 -                final_factors.insert(0, Expr::Rational(crate::expr::Rational::new(numerator, denominator)));
 +                final_factors.insert(
 +                    0,
 +                    Expr::Rational(crate::expr::Rational::new(numerator, denominator)),
 +                );
+             }
+         }
              Expr::Mul(final_factors)
+         }
+     }
++
 +    fn factorial_i64_checked(n: i64) -> Option<i64> {
 +        let n = u64::try_from(n).ok()?;
 +        let mut acc: i64 = 1;
 +        for k in 2..=n {
 +            let step = i64::try_from(k).ok()?;
 +            acc = acc.checked_mul(step)?;
 +        }
 +        Some(acc)
 +    }
+ }
  /// Equation solver
                          "Variable coefficient is zero".to_string(),
                      ));
+                 }
 -                return Ok(vec![Simplifier::simplify(&Expr::div(
 -                    Expr::neg(c),
 -                    b,
 -                ))]);
 +                return Ok(vec![Simplifier::simplify(&Expr::div(Expr::neg(c), b))]);
+             }
              // Quadratic
              return Self::solve_quadratic_with_coeffs(&a, &b, &c);
          // Try many starting points to find all roots
          let mut starting_points: Vec<f64> = vec![
 -            0.0, 1.0, -1.0, 2.0, -2.0, 0.5, -0.5, 3.0, -3.0, 4.0, -4.0, 5.0, -5.0,
 -            0.25, -0.25, 0.75, -0.75, 1.5, -1.5, 2.5, -2.5, 10.0, -10.0, 100.0, -100.0,
 +            0.0, 1.0, -1.0, 2.0, -2.0, 0.5, -0.5, 3.0, -3.0, 4.0, -4.0, 5.0, -5.0, 0.25, -0.25,
 +            0.75, -0.75, 1.5, -1.5, 2.5, -2.5, 10.0, -10.0, 100.0, -100.0,
          ];
          // Add more points based on degree
          for i in 0..20 {
+         }
          // Convert to expressions, cleaning up near-integers
 -        let result: Vec<Expr> = roots.into_iter().map(|r| {
 -            if r.abs() < 1e-10 {
 -                Expr::Integer(0)
 -            } else if (r - r.round()).abs() < 1e-10 {
 -                Expr::Integer(r.round() as i64)
 -            } else {
 -                Expr::Float(r)
 -            }
 -        }).collect();
 +        let result: Vec<Expr> = roots
 +            .into_iter()
 +            .map(|r| {
 +                if r.abs() < 1e-10 {
 +                    Expr::Integer(0)
 +                } else if (r - r.round()).abs() < 1e-10 {
 +                    Expr::Integer(r.round() as i64)
 +                } else {
 +                    Expr::Float(r)
 +                }
 +            })
 +            .collect();
          Some(result)
+     }
      /// Newton-Raphson iteration to find a root
 -    fn newton_raphson(expr: &Expr, var: &Symbol, start: f64, evaluator: &crate::eval::Evaluator) -> Option<f64> {
 +    fn newton_raphson(
 +        expr: &Expr,
 +        var: &Symbol,
 +        start: f64,
 +        evaluator: &crate::eval::Evaluator,
 +    ) -> Option<f64> {
          let deriv = Differentiator::diff(expr, var).ok()?;
          let deriv = Simplifier::simplify(&deriv);
          for _ in 0..max_iter {
              let f_x = {
                  let subst = Simplifier::substitute(expr, var, &Expr::Float(x));
 -                evaluator.eval(&subst).ok().and_then(|e| Self::expr_to_f64(&e))?
 +                evaluator
 +                    .eval(&subst)
 +                    .ok()
 +                    .and_then(|e| Self::expr_to_f64(&e))?
              };
              if f_x.abs() < tolerance {
              let fp_x = {
                  let subst = Simplifier::substitute(&deriv, var, &Expr::Float(x));
 -                evaluator.eval(&subst).ok().and_then(|e| Self::expr_to_f64(&e))?
 +                evaluator
 +                    .eval(&subst)
 +                    .ok()
 +                    .and_then(|e| Self::expr_to_f64(&e))?
              };
              if fp_x.abs() < 1e-15 {
+         }
+     }
 -    fn classify_polynomial(expr: &Expr, var: &Symbol) -> Option<(u32, Expr, Expr)> {
 -        // Very simplified polynomial detection
 -        // Returns (degree, leading coefficient, rest)
 -        match expr {
 -            // x alone: degree 1, coef 1
 -            Expr::Symbol(s) if s == var => Some((1, Expr::Integer(1), Expr::Integer(0))),
+-
 -            // a*x: degree 1, coef a
 -            Expr::Mul(factors) => {
 -                let mut coef = Expr::Integer(1);
 -                let mut has_var = false;
 -                let mut var_power = 1u32;
+-
 -                for f in factors {
 -                    match f {
 -                        Expr::Symbol(s) if s == var => {
 -                            has_var = true;
 -                        }
 -                        Expr::Pow(base, exp) => {
 -                            if let Expr::Symbol(s) = &**base {
 -                                if s == var {
 -                                    if let Expr::Integer(n) = **exp {
 -                                        has_var = true;
 -                                        var_power = n as u32;
 -                                    }
 -                                }
 -                            }
 -                        }
 -                        other if !other.contains_var(var) => {
 -                            coef = Expr::mul(vec![coef, other.clone()]);
 -                        }
 -                        _ => return None,
 -                    }
 -                }
+-
 -                if has_var {
 -                    Some((var_power, Simplifier::simplify(&coef), Expr::Integer(0)))
 -                } else {
 -                    None
 -                }
 -            }
+-
 -            // ax + b or ax^2 + bx + c
 -            Expr::Add(terms) => {
 -                let mut max_degree = 0u32;
 -                let mut leading_coef = Expr::Integer(0);
+-
 -                for term in terms {
 -                    if let Some((deg, coef, _)) = Self::classify_polynomial(term, var) {
 -                        if deg > max_degree {
 -                            max_degree = deg;
 -                            leading_coef = coef;
 -                        }
 -                    } else if !term.contains_var(var) {
 -                        // Constant term, ignore for degree calculation
 -                    } else {
 -                        return None;
 -                    }
 -                }
+-
 -                if max_degree > 0 {
 -                    Some((max_degree, leading_coef, Expr::Integer(0)))
 -                } else {
 -                    None
 -                }
 -            }
+-
 -            _ => None,
 -        }
 -    }
+-
      fn solve_quadratic_with_coeffs(a: &Expr, b: &Expr, c: &Expr) -> Result<Vec<Expr>> {
          // Discriminant: b^2 - 4ac
          let discriminant = Simplifier::simplify(&Expr::sub(
              Expr::add(vec![neg_b.clone(), sqrt_d.clone()]),
              two_a.clone(),
          ));
 -        let x2 = Simplifier::simplify(&Expr::div(
 -            Expr::sub(neg_b, sqrt_d),
 -            two_a,
 -        ));
 +        let x2 = Simplifier::simplify(&Expr::div(Expr::sub(neg_b, sqrt_d), two_a));
          Ok(vec![x1, x2])
+     }
          Self::solve_quadratic_with_coeffs(&a, &b, &c)
+     }
 -    fn extract_quadratic_coefficients(
 -        expr: &Expr,
 -        var: &Symbol,
 -    ) -> Result<(Expr, Expr, Expr)> {
 +    fn extract_quadratic_coefficients(expr: &Expr, var: &Symbol) -> Result<(Expr, Expr, Expr)> {
          // Initialize coefficients
          let mut a = Expr::Integer(0);
          let mut b = Expr::Integer(0);
                  _ => {
                      return Err(CasError::EvaluationError(
                          "Not a quadratic equation".to_string(),
 -                    ))
 +                    ));
+                 }
+             }
+         }
          assert_eq!(result, Expr::symbol("x"));
+     }
 +    #[test]
 +    fn test_simplify_mul_rational_one_identity() {
 +        let expr = Expr::mul(vec![
 +            Expr::Rational(crate::expr::Rational::new(1, 1)),
 +            Expr::symbol("n"),
 +        ]);
 +        let result = Simplifier::simplify(&expr);
 +        assert_eq!(result, Expr::symbol("n"));
 +    }
++
      #[test]
      fn test_simplify_combine_like() {
          // 2x + 3x = 5x
          let solutions = Solver::solve(&expr, &Symbol::new("x")).unwrap();
          assert!(!solutions.is_empty());
+     }
++
 +    #[test]
 +    fn test_simplify_factorial_constant() {
 +        let expr = Expr::func("factorial", vec![Expr::Integer(1)]);
 +        let simplified = Simplifier::simplify(&expr);
 +        assert_eq!(simplified, Expr::Integer(1));
 +    }
++
 +    #[test]
 +    fn test_simplify_product_residual_factorial_inverse() {
 +        let expr = Expr::mul(vec![
 +            Expr::func("factorial", vec![Expr::symbol("n")]),
 +            Expr::func(
 +                "factorial",
 +                vec![Expr::add(vec![Expr::Integer(1), Expr::symbol("n")])],
 +            ),
 +            Expr::pow(
 +                Expr::func("factorial", vec![Expr::Integer(1)]),
 +                Expr::Integer(-1),
 +            ),
 +        ]);
++
 +        let simplified = Simplifier::simplify(&expr);
 +        let rendered = simplified.to_string();
 +        assert!(!rendered.contains("1!"));
 +        assert!(!rendered.contains("factorial"));
 +        assert!(rendered.contains("n!"));
 +    }
+ }