gardesk/garcalc / e8ccf3f

+[package]

+name = "garcalc-cas"

+version.workspace = true

+edition.workspace = true

+authors.workspace = true

+license.workspace = true

+description = "Computer algebra system for garcalc"

+

+[dependencies]

+serde = { workspace = true }

+thiserror = { workspace = true }

+

+[dev-dependencies]

+use thiserror::Error;

+

+/// CAS error types

+#[derive(Debug, Error)]

+pub enum CasError {

+    #[error("Parse error at position {position}: {message}")]

+    Parse { position: usize, message: String },

+

+    #[error("Undefined variable: {0}")]

+    UndefinedVariable(String),

+

+    #[error("Undefined function: {0}")]

+    UndefinedFunction(String),

+

+    #[error("Type error: {0}")]

+    Type(String),

+

+    #[error("Division by zero")]

+    DivisionByZero,

+

+    #[error("Domain error: {0}")]

+    Domain(String),

+

+    #[error("Not implemented: {0}")]

+    NotImplemented(String),

+

+    #[error("Invalid argument: {0}")]

+    InvalidArgument(String),

+}

+

+pub type Result<T> = std::result::Result<T, CasError>;

+//! Expression evaluation

+//!

+//! Supports both numeric and symbolic evaluation modes.

+

+use std::collections::HashMap;

+use std::f64::consts::{E, PI};

+

+use crate::error::{CasError, Result};

+use crate::expr::{Expr, Rational, Sign};

+

+/// Variable bindings for evaluation

+pub type Environment = HashMap<String, Expr>;

+

+/// Expression evaluator

+pub struct Evaluator {

+    env: Environment,

+    /// If true, try to keep results symbolic when possible

+    exact_mode: bool,

+    /// Angle mode for trig functions

+    angle_mode: AngleMode,

+}

+

+#[derive(Debug, Clone, Copy, PartialEq, Eq)]

+pub enum AngleMode {

+    Radians,

+    Degrees,

+}

+

+impl Default for AngleMode {

+    fn default() -> Self {

+        Self::Radians

+    }

+}

+

+impl Default for Evaluator {

+    fn default() -> Self {

+        Self::new()

+    }

+}

+

+impl Evaluator {

+    pub fn new() -> Self {

+        Self {

+            env: Environment::new(),

+            exact_mode: false,

+            angle_mode: AngleMode::Radians,

+        }

+    }

+

+    pub fn with_exact_mode(mut self, exact: bool) -> Self {

+        self.exact_mode = exact;

+        self

+    }

+

+    pub fn with_angle_mode(mut self, mode: AngleMode) -> Self {

+        self.angle_mode = mode;

+        self

+    }

+

+    pub fn set_var(&mut self, name: impl Into<String>, value: Expr) {

+        self.env.insert(name.into(), value);

+    }

+

+    pub fn get_var(&self, name: &str) -> Option<&Expr> {

+        self.env.get(name)

+    }

+

+    pub fn clear_vars(&mut self) {

+        self.env.clear();

+    }

+

+    /// Evaluate an expression to a numeric result

+    pub fn eval(&self, expr: &Expr) -> Result<Expr> {

+        match expr {

+            Expr::Integer(n) => Ok(Expr::Integer(*n)),

+            Expr::Rational(r) => Ok(Expr::Rational(*r)),

+            Expr::Float(x) => Ok(Expr::Float(*x)),

+            Expr::Complex(re, im) => Ok(Expr::Complex(*re, *im)),

+

+            Expr::Symbol(sym) => {

+                // Check for constants

+                match sym.as_str() {

+                    "pi" => Ok(Expr::Float(PI)),

+                    "e" => Ok(Expr::Float(E)),

+                    _ => {

+                        // Check environment

+                        if let Some(value) = self.env.get(sym.as_str()) {

+                            self.eval(value)

+                        } else if self.exact_mode {

+                            // In exact mode, keep undefined symbols

+                            Ok(expr.clone())

+                        } else {

+                            Err(CasError::UndefinedVariable(sym.to_string()))

+                        }

+                    }

+                }

+            }

+

+            Expr::Neg(e) => {

+                let val = self.eval(e)?;

+                self.negate(&val)

+            }

+

+            Expr::Add(terms) => {

+                let mut sum = Expr::integer(0);

+                for term in terms {

+                    let val = self.eval(term)?;

+                    sum = self.add(&sum, &val)?;

+                }

+                Ok(sum)

+            }

+

+            Expr::Mul(factors) => {

+                let mut product = Expr::integer(1);

+                for factor in factors {

+                    let val = self.eval(factor)?;

+                    product = self.multiply(&product, &val)?;

+                }

+                Ok(product)

+            }

+

+            Expr::Pow(base, exp) => {

+                let base_val = self.eval(base)?;

+                let exp_val = self.eval(exp)?;

+                self.power(&base_val, &exp_val)

+            }

+

+            Expr::Func(name, args) => {

+                let evaluated_args: Result<Vec<_>> =

+                    args.iter().map(|a| self.eval(a)).collect();

+                self.call_function(name, &evaluated_args?)

+            }

+

+            Expr::Equation(lhs, rhs) => {

+                let lhs_val = self.eval(lhs)?;

+                let rhs_val = self.eval(rhs)?;

+                Ok(Expr::Equation(Box::new(lhs_val), Box::new(rhs_val)))

+            }

+

+            Expr::Vector(elems) => {

+                let evaluated: Result<Vec<_>> = elems.iter().map(|e| self.eval(e)).collect();

+                Ok(Expr::Vector(evaluated?))

+            }

+

+            Expr::Matrix(rows) => {

+                let evaluated: Result<Vec<Vec<_>>> = rows

+                    .iter()

+                    .map(|row| row.iter().map(|e| self.eval(e)).collect())

+                    .collect();

+                Ok(Expr::Matrix(evaluated?))

+            }

+

+            // Symbolic operations - keep as-is for now (will implement in later sprints)

+            Expr::Derivative { .. }

+            | Expr::Integral { .. }

+            | Expr::Limit { .. }

+            | Expr::Sum { .. }

+            | Expr::Product { .. } => {

+                if self.exact_mode {

+                    Ok(expr.clone())

+                } else {

+                    Err(CasError::NotImplemented(

+                        "symbolic operations in numeric mode".to_string(),

+                    ))

+                }

+            }

+

+            Expr::Inequality { lhs, op, rhs } => {

+                let lhs_val = self.eval(lhs)?;

+                let rhs_val = self.eval(rhs)?;

+                Ok(Expr::Inequality {

+                    lhs: Box::new(lhs_val),

+                    op: *op,

+                    rhs: Box::new(rhs_val),

+                })

+            }

+

+            Expr::Undefined => Ok(Expr::Undefined),

+            Expr::Infinity(sign) => Ok(Expr::Infinity(*sign)),

+        }

+    }

+

+    fn to_f64(&self, expr: &Expr) -> Result<f64> {

+        match expr {

+            Expr::Integer(n) => Ok(*n as f64),

+            Expr::Rational(r) => Ok(r.to_f64()),

+            Expr::Float(x) => Ok(*x),

+            _ => Err(CasError::Type(format!("expected number, got {expr}"))),

+        }

+    }

+

+    fn negate(&self, expr: &Expr) -> Result<Expr> {

+        match expr {

+            Expr::Integer(n) => Ok(Expr::Integer(-n)),

+            Expr::Rational(r) => Ok(Expr::Rational(Rational::new(-r.num, r.den))),

+            Expr::Float(x) => Ok(Expr::Float(-x)),

+            Expr::Complex(re, im) => Ok(Expr::Complex(-re, -im)),

+            Expr::Infinity(Sign::Positive) => Ok(Expr::Infinity(Sign::Negative)),

+            Expr::Infinity(Sign::Negative) => Ok(Expr::Infinity(Sign::Positive)),

+            _ => Ok(Expr::neg(expr.clone())),

+        }

+    }

+

+    fn add(&self, a: &Expr, b: &Expr) -> Result<Expr> {

+        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)),

+            (Expr::Integer(x), Expr::Float(y)) | (Expr::Float(y), Expr::Integer(x)) => {

+                Ok(Expr::Float(*x as f64 + y))

+            }

+            (Expr::Complex(r1, i1), Expr::Complex(r2, i2)) => Ok(Expr::Complex(r1 + r2, i1 + i2)),

+            (Expr::Complex(re, im), Expr::Float(x)) | (Expr::Float(x), Expr::Complex(re, im)) => {

+                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::Rational(r1), Expr::Rational(r2)) => {

+                let num = r1.num * r2.den + r2.num * r1.den;

+                let den = r1.den * r2.den;

+                Ok(Expr::Rational(Rational::new(num, den)))

+            }

+            (Expr::Rational(r), Expr::Integer(n)) | (Expr::Integer(n), Expr::Rational(r)) => {

+                let num = r.num + n * r.den;

+                Ok(Expr::Rational(Rational::new(num, r.den)))

+            }

+            _ => {

+                // Try converting to floats

+                if let (Ok(x), Ok(y)) = (self.to_f64(a), self.to_f64(b)) {

+                    Ok(Expr::Float(x + y))

+                } else if self.exact_mode {

+                    Ok(Expr::add(vec![a.clone(), b.clone()]))

+                } else {

+                    Err(CasError::Type(format!("cannot add {a} and {b}")))

+                }

+            }

+        }

+    }

+

+    fn multiply(&self, a: &Expr, b: &Expr) -> Result<Expr> {

+        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)),

+            (Expr::Integer(x), Expr::Float(y)) | (Expr::Float(y), Expr::Integer(x)) => {

+                Ok(Expr::Float(*x as f64 * y))

+            }

+            (Expr::Complex(r1, i1), Expr::Complex(r2, i2)) => {

+                // (a+bi)(c+di) = (ac-bd) + (ad+bc)i

+                Ok(Expr::Complex(r1 * r2 - i1 * i2, r1 * i2 + i1 * r2))

+            }

+            (Expr::Complex(re, im), Expr::Float(x)) | (Expr::Float(x), Expr::Complex(re, im)) => {

+                Ok(Expr::Complex(re * x, im * x))

+            }

+            (Expr::Complex(re, im), Expr::Integer(n))

+            | (Expr::Integer(n), Expr::Complex(re, im)) => {

+                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(r), Expr::Integer(n)) | (Expr::Integer(n), Expr::Rational(r)) => {

+                Ok(Expr::Rational(Rational::new(r.num * n, r.den)))

+            }

+            _ => {

+                if let (Ok(x), Ok(y)) = (self.to_f64(a), self.to_f64(b)) {

+                    Ok(Expr::Float(x * y))

+                } else if self.exact_mode {

+                    Ok(Expr::mul(vec![a.clone(), b.clone()]))

+                } else {

+                    Err(CasError::Type(format!("cannot multiply {a} and {b}")))

+                }

+            }

+        }

+    }

+

+    fn power(&self, base: &Expr, exp: &Expr) -> Result<Expr> {

+        // Special cases

+        if exp.is_zero() {

+            return Ok(Expr::integer(1));

+        }

+        if exp.is_one() {

+            return Ok(base.clone());

+        }

+        if base.is_zero() {

+            return Ok(Expr::integer(0));

+        }

+        if base.is_one() {

+            return Ok(Expr::integer(1));

+        }

+

+        match (base, exp) {

+            (Expr::Integer(b), Expr::Integer(e)) => {

+                if *e >= 0 {

+                    Ok(Expr::Integer(b.pow(*e as u32)))

+                } else {

+                    // Negative exponent -> rational or float

+                    let denom = b.pow((-e) as u32);

+                    if self.exact_mode {

+                        Ok(Expr::Rational(Rational::new(1, denom)))

+                    } else {

+                        Ok(Expr::Float(1.0 / denom as f64))

+                    }

+                }

+            }

+            (Expr::Float(b), Expr::Integer(e)) => Ok(Expr::Float(b.powi(*e as i32))),

+            (Expr::Float(b), Expr::Float(e)) => Ok(Expr::Float(b.powf(*e))),

+            (Expr::Integer(b), Expr::Float(e)) => Ok(Expr::Float((*b as f64).powf(*e))),

+            _ => {

+                if let (Ok(b), Ok(e)) = (self.to_f64(base), self.to_f64(exp)) {

+                    Ok(Expr::Float(b.powf(e)))

+                } else if self.exact_mode {

+                    Ok(Expr::pow(base.clone(), exp.clone()))

+                } else {

+                    Err(CasError::Type(format!("cannot compute {base}^{exp}")))

+                }

+            }

+        }

+    }

+

+    fn call_function(&self, name: &str, args: &[Expr]) -> Result<Expr> {

+        // Get numeric argument if single-arg function

+        let arg = if args.len() == 1 {

+            self.to_f64(&args[0]).ok()

+        } else {

+            None

+        };

+

+        // Angle conversion for trig functions

+        let angle = |x: f64| match self.angle_mode {

+            AngleMode::Radians => x,

+            AngleMode::Degrees => x.to_radians(),

+        };

+

+        let from_angle = |x: f64| match self.angle_mode {

+            AngleMode::Radians => x,

+            AngleMode::Degrees => x.to_degrees(),

+        };

+

+        match (name, args.len(), arg) {

+            // Trigonometric

+            ("sin", 1, Some(x)) => Ok(Expr::Float(angle(x).sin())),

+            ("cos", 1, Some(x)) => Ok(Expr::Float(angle(x).cos())),

+            ("tan", 1, Some(x)) => Ok(Expr::Float(angle(x).tan())),

+            ("asin", 1, Some(x)) => Ok(Expr::Float(from_angle(x.asin()))),

+            ("acos", 1, Some(x)) => Ok(Expr::Float(from_angle(x.acos()))),

+            ("atan", 1, Some(x)) => Ok(Expr::Float(from_angle(x.atan()))),

+            ("sinh", 1, Some(x)) => Ok(Expr::Float(x.sinh())),

+            ("cosh", 1, Some(x)) => Ok(Expr::Float(x.cosh())),

+            ("tanh", 1, Some(x)) => Ok(Expr::Float(x.tanh())),

+            ("asinh", 1, Some(x)) => Ok(Expr::Float(x.asinh())),

+            ("acosh", 1, Some(x)) => Ok(Expr::Float(x.acosh())),

+            ("atanh", 1, Some(x)) => Ok(Expr::Float(x.atanh())),

+

+            // Exponential/logarithmic

+            ("exp", 1, Some(x)) => Ok(Expr::Float(x.exp())),

+            ("ln", 1, Some(x)) => Ok(Expr::Float(x.ln())),

+            ("log", 1, Some(x)) => Ok(Expr::Float(x.log10())),

+            ("log10", 1, Some(x)) => Ok(Expr::Float(x.log10())),

+            ("log2", 1, Some(x)) => Ok(Expr::Float(x.log2())),

+

+            // Roots

+            ("sqrt", 1, Some(x)) => {

+                if x >= 0.0 {

+                    Ok(Expr::Float(x.sqrt()))

+                } else {

+                    Ok(Expr::Complex(0.0, (-x).sqrt()))

+                }

+            }

+            ("cbrt", 1, Some(x)) => Ok(Expr::Float(x.cbrt())),

+

+            // Other

+            ("abs", 1, Some(x)) => Ok(Expr::Float(x.abs())),

+            ("floor", 1, Some(x)) => Ok(Expr::Integer(x.floor() as i64)),

+            ("ceil", 1, Some(x)) => Ok(Expr::Integer(x.ceil() as i64)),

+            ("round", 1, Some(x)) => Ok(Expr::Integer(x.round() as i64)),

+            ("sign", 1, Some(x)) => Ok(Expr::Integer(if x > 0.0 {

+                1

+            } else if x < 0.0 {

+                -1

+            } else {

+                0

+            })),

+

+            // Factorial

+            ("factorial", 1, _) => {

+                if let Expr::Integer(n) = &args[0] {

+                    if *n < 0 {

+                        Err(CasError::Domain("factorial of negative number".to_string()))

+                    } else if *n > 20 {

+                        // Use Stirling's approximation for large n

+                        Ok(Expr::Float(gamma(*n as f64 + 1.0)))

+                    } else {

+                        Ok(Expr::Integer(factorial(*n as u64) as i64))

+                    }

+                } else if let Some(x) = arg {

+                    Ok(Expr::Float(gamma(x + 1.0)))

+                } else {

+                    Err(CasError::Type("factorial requires a number".to_string()))

+                }

+            }

+

+            // Two-argument functions

+            ("atan2", 2, _) => {

+                let y = self.to_f64(&args[0])?;

+                let x = self.to_f64(&args[1])?;

+                Ok(Expr::Float(from_angle(y.atan2(x))))

+            }

+            ("pow", 2, _) => self.power(&args[0], &args[1]),

+            ("mod", 2, _) => {

+                let a = self.to_f64(&args[0])?;

+                let b = self.to_f64(&args[1])?;

+                Ok(Expr::Float(a % b))

+            }

+            ("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));

+                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));

+                Ok(Expr::Float(max))

+            }

+            ("gcd", 2, _) => {

+                if let (Expr::Integer(a), Expr::Integer(b)) = (&args[0], &args[1]) {

+                    Ok(Expr::Integer(gcd(*a, *b)))

+                } else {

+                    Err(CasError::Type("gcd requires integers".to_string()))

+                }

+            }

+            ("lcm", 2, _) => {

+                if let (Expr::Integer(a), Expr::Integer(b)) = (&args[0], &args[1]) {

+                    Ok(Expr::Integer(lcm(*a, *b)))

+                } else {

+                    Err(CasError::Type("lcm requires integers".to_string()))

+                }

+            }

+

+            _ => {

+                if self.exact_mode {

+                    Ok(Expr::func(name, args.to_vec()))

+                } else {

+                    Err(CasError::UndefinedFunction(name.to_string()))

+                }

+            }

+        }

+    }

+}

+

+// Helper functions

+

+fn factorial(n: u64) -> u64 {

+    (1..=n).product()

+}

+

+fn gcd(a: i64, b: i64) -> i64 {

+    let (a, b) = (a.abs(), b.abs());

+    if b == 0 {

+        a

+    } else {

+        gcd(b, a % b)

+    }

+}

+

+fn lcm(a: i64, b: i64) -> i64 {

+    (a * b).abs() / gcd(a, b)

+}

+

+/// Gamma function approximation (Lanczos)

+fn gamma(x: f64) -> f64 {

+    if x < 0.5 {

+        PI / (PI * x).sin() / gamma(1.0 - x)

+    } else {

+        let x = x - 1.0;

+        let g = 7.0;

+        let c = [

+            0.99999999999980993,

+            676.5203681218851,

+            -1259.1392167224028,

+            771.32342877765313,

+            -176.61502916214059,

+            12.507343278686905,

+            -0.13857109526572012,

+            9.9843695780195716e-6,

+            1.5056327351493116e-7,

+        ];

+

+        let mut sum = c[0];

+        for (i, &ci) in c.iter().enumerate().skip(1) {

+            sum += ci / (x + i as f64);

+        }

+

+        (2.0 * PI).sqrt() * (x + g + 0.5).powf(x + 0.5) * (-(x + g + 0.5)).exp() * sum

+    }

+}

+

+#[cfg(test)]

+mod tests {

+    use super::*;

+    use crate::parser::parse;

+

+    fn eval(input: &str) -> Result<Expr> {

+        let expr = parse(input)?;

+        Evaluator::new().eval(&expr)

+    }

+

+    fn eval_to_f64(input: &str) -> f64 {

+        match eval(input).unwrap() {

+            Expr::Integer(n) => n as f64,

+            Expr::Float(x) => x,

+            other => panic!("expected number, got {other}"),

+        }

+    }

+

+    #[test]

+    fn test_arithmetic() {

+        assert_eq!(eval("2 + 3").unwrap(), Expr::Integer(5));

+        assert_eq!(eval("10 - 4").unwrap(), Expr::Integer(6));

+        assert_eq!(eval("3 * 4").unwrap(), Expr::Integer(12));

+        assert_eq!(eval("2^10").unwrap(), Expr::Integer(1024));

+    }

+

+    #[test]

+    fn test_float() {

+        let result = eval_to_f64("3.14 * 2");

+        assert!((result - 6.28).abs() < 1e-10);

+    }

+

+    #[test]

+    fn test_functions() {

+        let result = eval_to_f64("sin(0)");

+        assert!(result.abs() < 1e-10);

+

+        let result = eval_to_f64("cos(0)");

+        assert!((result - 1.0).abs() < 1e-10);

+

+        let result = eval_to_f64("sqrt(4)");

+        assert!((result - 2.0).abs() < 1e-10);

+

+        let result = eval_to_f64("ln(e)");

+        assert!((result - 1.0).abs() < 1e-10);

+    }

+

+    #[test]

+    fn test_constants() {

+        let result = eval_to_f64("pi");

+        assert!((result - PI).abs() < 1e-10);

+

+        let result = eval_to_f64("e");

+        assert!((result - E).abs() < 1e-10);

+    }

+

+    #[test]

+    fn test_factorial() {

+        assert_eq!(eval("5!").unwrap(), Expr::Integer(120));

+    }

+

+    #[test]

+    fn test_complex_expr() {

+        let result = eval_to_f64("2 * sin(pi/2) + 1");

+        assert!((result - 3.0).abs() < 1e-10);

+    }

+

+    #[test]

+    fn test_variables() {

+        let expr = parse("x + 1").unwrap();

+        let mut evaluator = Evaluator::new();

+        evaluator.set_var("x", Expr::integer(5));

+        let result = evaluator.eval(&expr).unwrap();

+        assert_eq!(result, Expr::Integer(6));

+    }

+}

+use serde::{Deserialize, Serialize};

+use std::fmt;

+

+/// A symbol (variable name)

+#[derive(Debug, Clone, PartialEq, Eq, Hash, Serialize, Deserialize)]

+pub struct Symbol(pub String);

+

+impl Symbol {

+    pub fn new(name: impl Into<String>) -> Self {

+        Self(name.into())

+    }

+

+    pub fn as_str(&self) -> &str {

+        &self.0

+    }

+}

+

+impl fmt::Display for Symbol {

+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {

+        write!(f, "{}", self.0)

+    }

+}

+

+impl From<&str> for Symbol {

+    fn from(s: &str) -> Self {

+        Self(s.to_string())

+    }

+}

+

+/// A rational number (numerator/denominator)

+#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Serialize, Deserialize)]

+pub struct Rational {

+    pub num: i64,

+    pub den: i64,

+}

+

+impl Rational {

+    pub fn new(num: i64, den: i64) -> Self {

+        let g = gcd(num.abs(), den.abs());

+        let sign = if den < 0 { -1 } else { 1 };

+        Self {

+            num: sign * num / g,

+            den: den.abs() / g,

+        }

+    }

+

+    pub fn integer(n: i64) -> Self {

+        Self { num: n, den: 1 }

+    }

+

+    pub fn is_integer(&self) -> bool {

+        self.den == 1

+    }

+

+    pub fn to_f64(self) -> f64 {

+        self.num as f64 / self.den as f64

+    }

+}

+

+impl fmt::Display for Rational {

+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {

+        if self.den == 1 {

+            write!(f, "{}", self.num)

+        } else {

+            write!(f, "{}/{}", self.num, self.den)

+        }

+    }

+}

+

+fn gcd(a: i64, b: i64) -> i64 {

+    if b == 0 {

+        a

+    } else {

+        gcd(b, a % b)

+    }

+}

+

+/// Mathematical expression AST

+#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]

+pub enum Expr {

+    // Atoms

+    /// Integer value

+    Integer(i64),

+    /// Rational number

+    Rational(Rational),

+    /// Floating point

+    Float(f64),

+    /// Complex number (real, imaginary)

+    Complex(f64, f64),

+    /// Symbolic variable

+    Symbol(Symbol),

+

+    // Arithmetic operations

+    /// Negation: -x

+    Neg(Box<Expr>),

+    /// Addition: a + b + c + ...

+    Add(Vec<Expr>),

+    /// Multiplication: a * b * c * ...

+    Mul(Vec<Expr>),

+    /// Power: base^exponent

+    Pow(Box<Expr>, Box<Expr>),

+

+    // Function application

+    /// Function call: f(args...)

+    Func(String, Vec<Expr>),

+

+    // Calculus (symbolic)

+    /// Derivative: d/dx(expr) with order

+    Derivative {

+        expr: Box<Expr>,

+        var: Symbol,

+        order: u32,

+    },

+    /// Integral: ∫ expr dx, optionally with bounds

+    Integral {

+        expr: Box<Expr>,

+        var: Symbol,

+        lower: Option<Box<Expr>>,

+        upper: Option<Box<Expr>>,

+    },

+    /// Limit: lim_{x→point} expr

+    Limit {

+        expr: Box<Expr>,

+        var: Symbol,

+        point: Box<Expr>,

+        direction: Option<LimitDirection>,

+    },

+    /// Summation: Σ_{var=lower}^{upper} expr

+    Sum {

+        expr: Box<Expr>,

+        var: Symbol,

+        lower: Box<Expr>,

+        upper: Box<Expr>,

+    },

+    /// Product: Π_{var=lower}^{upper} expr

+    Product {

+        expr: Box<Expr>,

+        var: Symbol,

+        lower: Box<Expr>,

+        upper: Box<Expr>,

+    },

+

+    // Algebra

+    /// Equation: lhs = rhs

+    Equation(Box<Expr>, Box<Expr>),

+    /// Inequality: lhs < rhs (or >, <=, >=, !=)

+    Inequality {

+        lhs: Box<Expr>,

+        op: InequalityOp,

+        rhs: Box<Expr>,

+    },

+

+    // Linear algebra

+    /// Matrix: rows of expressions

+    Matrix(Vec<Vec<Expr>>),

+    /// Vector: list of expressions

+    Vector(Vec<Expr>),

+

+    // Special values

+    /// Undefined result

+    Undefined,

+    /// Positive or negative infinity

+    Infinity(Sign),

+}

+

+/// Direction for limits

+#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]

+pub enum LimitDirection {

+    Left,

+    Right,

+}

+

+/// Sign for infinity

+#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]

+pub enum Sign {

+    Positive,

+    Negative,

+}

+

+/// Inequality operators

+#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]

+pub enum InequalityOp {

+    Lt,

+    Le,

+    Gt,

+    Ge,

+    Ne,

+}

+

+impl Expr {

+    // Constructors for common expressions

+

+    pub fn integer(n: i64) -> Self {

+        Self::Integer(n)

+    }

+

+    pub fn float(x: f64) -> Self {

+        Self::Float(x)

+    }

+

+    pub fn symbol(name: impl Into<String>) -> Self {

+        Self::Symbol(Symbol::new(name))

+    }

+

+    pub fn neg(expr: Expr) -> Self {

+        Self::Neg(Box::new(expr))

+    }

+

+    pub fn add(terms: Vec<Expr>) -> Self {

+        if terms.len() == 1 {

+            terms.into_iter().next().unwrap()

+        } else {

+            Self::Add(terms)

+        }

+    }

+

+    pub fn sub(a: Expr, b: Expr) -> Self {

+        Self::add(vec![a, Self::neg(b)])

+    }

+

+    pub fn mul(factors: Vec<Expr>) -> Self {

+        if factors.len() == 1 {

+            factors.into_iter().next().unwrap()

+        } else {

+            Self::Mul(factors)

+        }

+    }

+

+    pub fn div(a: Expr, b: Expr) -> Self {

+        Self::mul(vec![a, Self::pow(b, Self::integer(-1))])

+    }

+

+    pub fn pow(base: Expr, exp: Expr) -> Self {

+        Self::Pow(Box::new(base), Box::new(exp))

+    }

+

+    pub fn func(name: impl Into<String>, args: Vec<Expr>) -> Self {

+        Self::Func(name.into(), args)

+    }

+

+    // Common mathematical constants

+    pub fn pi() -> Self {

+        Self::symbol("pi")

+    }

+

+    pub fn e() -> Self {

+        Self::symbol("e")

+    }

+

+    pub fn i() -> Self {

+        Self::Complex(0.0, 1.0)

+    }

+

+    // Predicate methods

+

+    pub fn is_zero(&self) -> bool {

+        match self {

+            Self::Integer(0) => true,

+            Self::Float(x) if *x == 0.0 => true,

+            _ => false,

+        }

+    }

+

+    pub fn is_one(&self) -> bool {

+        match self {

+            Self::Integer(1) => true,

+            Self::Float(x) if *x == 1.0 => true,

+            _ => false,

+        }

+    }

+

+    pub fn is_negative_one(&self) -> bool {

+        match self {

+            Self::Integer(-1) => true,

+            Self::Float(x) if *x == -1.0 => true,

+            _ => false,

+        }

+    }

+

+    pub fn is_number(&self) -> bool {

+        matches!(

+            self,

+            Self::Integer(_) | Self::Rational(_) | Self::Float(_) | Self::Complex(_, _)

+        )

+    }

+

+    pub fn is_symbol(&self) -> bool {

+        matches!(self, Self::Symbol(_))

+    }

+

+    /// Check if expression contains a variable

+    pub fn contains_var(&self, var: &Symbol) -> bool {

+        match self {

+            Self::Symbol(s) => s == var,

+            Self::Neg(e) => e.contains_var(var),

+            Self::Add(terms) => terms.iter().any(|t| t.contains_var(var)),

+            Self::Mul(factors) => factors.iter().any(|f| f.contains_var(var)),

+            Self::Pow(base, exp) => base.contains_var(var) || exp.contains_var(var),

+            Self::Func(_, args) => args.iter().any(|a| a.contains_var(var)),

+            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::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::Vector(elems) => elems.iter().any(|e| e.contains_var(var)),

+            _ => false,

+        }

+    }

+}

+

+impl fmt::Display for Expr {

+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {

+        match self {

+            Self::Integer(n) => write!(f, "{n}"),

+            Self::Rational(r) => write!(f, "{r}"),

+            Self::Float(x) => {

+                if x.fract() == 0.0 && x.abs() < 1e15 {

+                    write!(f, "{}", *x as i64)

+                } else {

+                    write!(f, "{x}")

+                }

+            }

+            Self::Complex(re, im) => {

+                if *re == 0.0 {

+                    if *im == 1.0 {

+                        write!(f, "i")

+                    } else if *im == -1.0 {

+                        write!(f, "-i")

+                    } else {

+                        write!(f, "{im}i")

+                    }

+                } else if *im >= 0.0 {

+                    write!(f, "{re}+{im}i")

+                } else {

+                    write!(f, "{re}{im}i")

+                }

+            }

+            Self::Symbol(s) => write!(f, "{s}"),

+            Self::Neg(e) => write!(f, "-{e}"),

+            Self::Add(terms) => {

+                if terms.is_empty() {

+                    write!(f, "0")

+                } else {

+                    write!(f, "(")?;

+                    for (i, term) in terms.iter().enumerate() {

+                        if i > 0 {

+                            write!(f, "+")?;

+                        }

+                        write!(f, "{term}")?;

+                    }

+                    write!(f, ")")

+                }

+            }

+            Self::Mul(factors) => {

+                if factors.is_empty() {

+                    write!(f, "1")

+                } else {

+                    write!(f, "(")?;

+                    for (i, factor) in factors.iter().enumerate() {

+                        if i > 0 {

+                            write!(f, "*")?;

+                        }

+                        write!(f, "{factor}")?;

+                    }

+                    write!(f, ")")

+                }

+            }

+            Self::Pow(base, exp) => write!(f, "{base}^{exp}"),

+            Self::Func(name, args) => {

+                write!(f, "{name}(")?;

+                for (i, arg) in args.iter().enumerate() {

+                    if i > 0 {

+                        write!(f, ", ")?;

+                    }

+                    write!(f, "{arg}")?;

+                }

+                write!(f, ")")

+            }

+            Self::Derivative { expr, var, order } => {

+                if *order == 1 {

+                    write!(f, "d/d{var}({expr})")

+                } else {

+                    write!(f, "d^{order}/d{var}^{order}({expr})")

+                }

+            }

+            Self::Integral {

+                expr,

+                var,

+                lower,

+                upper,

+            } => {

+                if let (Some(l), Some(u)) = (lower, upper) {

+                    write!(f, "integrate({expr}, {var}, {l}, {u})")

+                } else {

+                    write!(f, "integrate({expr}, {var})")

+                }

+            }

+            Self::Limit {

+                expr,

+                var,

+                point,

+                direction,

+            } => {

+                let dir = match direction {

+                    Some(LimitDirection::Left) => "-",

+                    Some(LimitDirection::Right) => "+",

+                    None => "",

+                };

+                write!(f, "lim({expr}, {var}, {point}{dir})")

+            }

+            Self::Sum {

+                expr,

+                var,

+                lower,

+                upper,

+            } => write!(f, "sum({expr}, {var}, {lower}, {upper})"),

+            Self::Product {

+                expr,

+                var,

+                lower,

+                upper,

+            } => write!(f, "product({expr}, {var}, {lower}, {upper})"),

+            Self::Equation(lhs, rhs) => write!(f, "{lhs} = {rhs}"),

+            Self::Inequality { lhs, op, rhs } => {

+                let op_str = match op {

+                    InequalityOp::Lt => "<",

+                    InequalityOp::Le => "<=",

+                    InequalityOp::Gt => ">",

+                    InequalityOp::Ge => ">=",

+                    InequalityOp::Ne => "!=",

+                };

+                write!(f, "{lhs} {op_str} {rhs}")

+            }

+            Self::Matrix(rows) => {

+                write!(f, "[")?;

+                for (i, row) in rows.iter().enumerate() {

+                    if i > 0 {

+                        write!(f, "; ")?;

+                    }

+                    for (j, elem) in row.iter().enumerate() {

+                        if j > 0 {

+                            write!(f, ", ")?;

+                        }

+                        write!(f, "{elem}")?;

+                    }

+                }

+                write!(f, "]")

+            }

+            Self::Vector(elems) => {

+                write!(f, "[")?;

+                for (i, elem) in elems.iter().enumerate() {

+                    if i > 0 {

+                        write!(f, ", ")?;

+                    }

+                    write!(f, "{elem}")?;

+                }

+                write!(f, "]")

+            }

+            Self::Undefined => write!(f, "undefined"),

+            Self::Infinity(Sign::Positive) => write!(f, "infinity"),

+            Self::Infinity(Sign::Negative) => write!(f, "-infinity"),

+        }

+    }

+}

+

+#[cfg(test)]

+mod tests {

+    use super::*;

+

+    #[test]

+    fn test_rational() {

+        let r = Rational::new(4, 6);

+        assert_eq!(r.num, 2);

+        assert_eq!(r.den, 3);

+        assert_eq!(r.to_string(), "2/3");

+    }

+

+    #[test]

+    fn test_expr_display() {

+        let expr = Expr::add(vec![

+            Expr::symbol("x"),

+            Expr::mul(vec![Expr::integer(2), Expr::symbol("y")]),

+        ]);

+        assert_eq!(expr.to_string(), "(x+(2*y))");

+    }

+

+    #[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")));

+    }

+}

+//! garcalc-cas: Computer Algebra System for garcalc

+//!

+//! A custom CAS engine providing:

+//! - Expression parsing and representation

+//! - Numeric and symbolic evaluation

+//! - Symbolic differentiation and integration

+//! - Equation solving

+//! - Limits and series expansions

+

+pub mod expr;

+pub mod parser;

+pub mod eval;

+pub mod error;

+

+pub use expr::{Expr, Symbol, Rational};

+pub use parser::Parser;

+pub use eval::Evaluator;

+pub use error::{CasError, Result};

+//! Expression parser using recursive descent

+//!

+//! Grammar:

+//! ```text

+//! expr       = equation

+//! equation   = additive ("=" additive)?

+//! additive   = multiplicative (("+"|"-") multiplicative)*

+//! multiplicative = power (("*"|"/") power)*

+//! power      = unary ("^" power)?

+//! unary      = "-" unary | postfix

+//! postfix    = primary ("!")*

+//! primary    = NUMBER | SYMBOL | "(" expr ")" | function | "[" matrix "]"

+//! function   = SYMBOL "(" args? ")"

+//! args       = expr ("," expr)*

+//! matrix     = row (";" row)*

+//! row        = expr ("," expr)*

+//! ```

+

+use crate::error::{CasError, Result};

+use crate::expr::{Expr, Symbol};

+

+/// Token types

+#[derive(Debug, Clone, PartialEq)]

+pub enum Token {

+    Number(f64),

+    Symbol(String),

+    Plus,

+    Minus,

+    Star,

+    Slash,

+    Caret,

+    Bang,

+    Eq,

+    LParen,

+    RParen,

+    LBracket,

+    RBracket,

+    Comma,

+    Semicolon,

+    Eof,

+}

+

+/// Tokenizer for mathematical expressions

+pub struct Lexer<'a> {

+    input: &'a str,

+    pos: usize,

+}

+

+impl<'a> Lexer<'a> {

+    pub fn new(input: &'a str) -> Self {

+        Self { input, pos: 0 }

+    }

+

+    fn peek_char(&self) -> Option<char> {

+        self.input[self.pos..].chars().next()

+    }

+

+    fn advance(&mut self) -> Option<char> {

+        let c = self.peek_char()?;

+        self.pos += c.len_utf8();

+        Some(c)

+    }

+

+    fn skip_whitespace(&mut self) {

+        while let Some(c) = self.peek_char() {

+            if c.is_whitespace() {

+                self.advance();

+            } else {

+                break;

+            }

+        }

+    }

+

+    fn read_number(&mut self) -> Token {

+        let start = self.pos;

+        let mut has_dot = false;

+        let mut has_exp = false;

+

+        while let Some(c) = self.peek_char() {

+            if c.is_ascii_digit() {

+                self.advance();

+            } else if c == '.' && !has_dot && !has_exp {

+                has_dot = true;

+                self.advance();

+            } else if (c == 'e' || c == 'E') && !has_exp {

+                has_exp = true;

+                self.advance();

+                // Handle optional sign after exponent

+                if let Some(sign) = self.peek_char() {

+                    if sign == '+' || sign == '-' {

+                        self.advance();

+                    }

+                }

+            } else {

+                break;

+            }

+        }

+

+        let num_str = &self.input[start..self.pos];

+        let value: f64 = num_str.parse().unwrap_or(f64::NAN);

+        Token::Number(value)

+    }

+

+    fn read_symbol(&mut self) -> Token {

+        let start = self.pos;

+

+        while let Some(c) = self.peek_char() {

+            if c.is_alphanumeric() || c == '_' {

+                self.advance();

+            } else {

+                break;

+            }

+        }

+

+        let name = self.input[start..self.pos].to_string();

+        Token::Symbol(name)

+    }

+

+    pub fn next_token(&mut self) -> Result<Token> {

+        self.skip_whitespace();

+

+        let Some(c) = self.peek_char() else {

+            return Ok(Token::Eof);

+        };

+

+        let token = match c {

+            '+' => {

+                self.advance();

+                Token::Plus

+            }

+            '-' => {

+                self.advance();

+                Token::Minus

+            }

+            '*' => {

+                self.advance();

+                Token::Star

+            }

+            '/' => {

+                self.advance();

+                Token::Slash

+            }

+            '^' => {

+                self.advance();

+                Token::Caret

+            }

+            '!' => {

+                self.advance();

+                Token::Bang

+            }

+            '=' => {

+                self.advance();

+                Token::Eq

+            }

+            '(' => {

+                self.advance();

+                Token::LParen

+            }

+            ')' => {

+                self.advance();

+                Token::RParen

+            }

+            '[' => {

+                self.advance();

+                Token::LBracket

+            }

+            ']' => {

+                self.advance();

+                Token::RBracket

+            }

+            ',' => {

+                self.advance();

+                Token::Comma

+            }

+            ';' => {

+                self.advance();

+                Token::Semicolon

+            }

+            _ if c.is_ascii_digit() || c == '.' => self.read_number(),

+            _ if c.is_alphabetic() || c == '_' => self.read_symbol(),

+            _ => {

+                return Err(CasError::Parse {

+                    position: self.pos,

+                    message: format!("unexpected character: '{c}'"),

+                });

+            }

+        };

+

+        Ok(token)

+    }

+

+    pub fn position(&self) -> usize {

+        self.pos

+    }

+}

+

+/// Recursive descent parser

+pub struct Parser<'a> {

+    lexer: Lexer<'a>,

+    current: Token,

+}

+

+impl<'a> Parser<'a> {

+    pub fn new(input: &'a str) -> Result<Self> {

+        let mut lexer = Lexer::new(input);

+        let current = lexer.next_token()?;

+        Ok(Self { lexer, current })

+    }

+

+    fn advance(&mut self) -> Result<Token> {

+        let prev = std::mem::replace(&mut self.current, self.lexer.next_token()?);

+        Ok(prev)

+    }

+

+    fn expect(&mut self, expected: Token) -> Result<()> {

+        if self.current == expected {

+            self.advance()?;

+            Ok(())

+        } else {

+            Err(CasError::Parse {

+                position: self.lexer.position(),

+                message: format!("expected {expected:?}, found {:?}", self.current),

+            })

+        }

+    }

+

+    /// Parse a complete expression

+    pub fn parse(&mut self) -> Result<Expr> {

+        let expr = self.parse_equation()?;

+        if self.current != Token::Eof {

+            return Err(CasError::Parse {

+                position: self.lexer.position(),

+                message: format!("unexpected token: {:?}", self.current),

+            });

+        }

+        Ok(expr)

+    }

+

+    fn parse_equation(&mut self) -> Result<Expr> {

+        let lhs = self.parse_additive()?;

+        if self.current == Token::Eq {

+            self.advance()?;

+            let rhs = self.parse_additive()?;

+            Ok(Expr::Equation(Box::new(lhs), Box::new(rhs)))

+        } else {

+            Ok(lhs)

+        }

+    }

+

+    fn parse_additive(&mut self) -> Result<Expr> {

+        let mut terms = vec![self.parse_multiplicative()?];

+

+        loop {

+            match &self.current {

+                Token::Plus => {

+                    self.advance()?;

+                    terms.push(self.parse_multiplicative()?);

+                }

+                Token::Minus => {

+                    self.advance()?;

+                    terms.push(Expr::neg(self.parse_multiplicative()?));

+                }

+                _ => break,

+            }

+        }

+

+        Ok(Expr::add(terms))

+    }

+

+    fn parse_multiplicative(&mut self) -> Result<Expr> {

+        let mut factors = vec![self.parse_power()?];

+

+        loop {

+            match &self.current {

+                Token::Star => {

+                    self.advance()?;

+                    factors.push(self.parse_power()?);

+                }

+                Token::Slash => {

+                    self.advance()?;

+                    let divisor = self.parse_power()?;

+                    factors.push(Expr::pow(divisor, Expr::integer(-1)));

+                }

+                // Implicit multiplication: 2x, xy, x(y+1)

+                Token::Symbol(_) | Token::LParen | Token::Number(_) => {

+                    // Check if this could be implicit multiplication

+                    if let Some(last) = factors.last() {

+                        if matches!(

+                            last,

+                            Expr::Symbol(_)

+                                | Expr::Integer(_)

+                                | Expr::Float(_)

+                                | Expr::Func(_, _)

+                        ) {

+                            factors.push(self.parse_power()?);

+                            continue;

+                        }

+                    }

+                    break;

+                }

+                _ => break,

+            }

+        }

+

+        Ok(Expr::mul(factors))

+    }

+

+    fn parse_power(&mut self) -> Result<Expr> {

+        let base = self.parse_unary()?;

+        if self.current == Token::Caret {

+            self.advance()?;

+            let exp = self.parse_power()?; // Right associative

+            Ok(Expr::pow(base, exp))

+        } else {

+            Ok(base)

+        }

+    }

+

+    fn parse_unary(&mut self) -> Result<Expr> {

+        if self.current == Token::Minus {

+            self.advance()?;

+            let expr = self.parse_unary()?;

+            Ok(Expr::neg(expr))

+        } else {

+            self.parse_postfix()

+        }

+    }

+

+    fn parse_postfix(&mut self) -> Result<Expr> {

+        let mut expr = self.parse_primary()?;

+

+        while self.current == Token::Bang {

+            self.advance()?;

+            expr = Expr::func("factorial", vec![expr]);

+        }

+

+        Ok(expr)

+    }

+

+    fn parse_primary(&mut self) -> Result<Expr> {

+        match &self.current {

+            Token::Number(n) => {

+                let n = *n;

+                self.advance()?;

+                // Convert to integer if possible

+                if n.fract() == 0.0 && n.abs() < i64::MAX as f64 {

+                    Ok(Expr::integer(n as i64))

+                } else {

+                    Ok(Expr::float(n))

+                }

+            }

+            Token::Symbol(name) => {

+                let name = name.clone();

+                self.advance()?;

+

+                // Check for function call

+                if self.current == Token::LParen {

+                    self.advance()?;

+                    let args = self.parse_args()?;

+                    self.expect(Token::RParen)?;

+                    Ok(self.handle_special_function(&name, args))

+                } else {

+                    // Constants

+                    match name.as_str() {

+                        "pi" | "PI" => Ok(Expr::symbol("pi")),

+                        "e" | "E" => Ok(Expr::symbol("e")),

+                        "i" => Ok(Expr::Complex(0.0, 1.0)),

+                        "inf" | "infinity" => Ok(Expr::Infinity(crate::expr::Sign::Positive)),

+                        _ => Ok(Expr::symbol(name)),

+                    }

+                }

+            }

+            Token::LParen => {

+                self.advance()?;

+                let expr = self.parse_equation()?;

+                self.expect(Token::RParen)?;

+                Ok(expr)

+            }

+            Token::LBracket => {

+                self.advance()?;

+                let matrix = self.parse_matrix()?;

+                self.expect(Token::RBracket)?;

+                Ok(matrix)

+            }

+            _ => Err(CasError::Parse {

+                position: self.lexer.position(),

+                message: format!("unexpected token: {:?}", self.current),

+            }),

+        }

+    }

+

+    fn parse_args(&mut self) -> Result<Vec<Expr>> {

+        if self.current == Token::RParen {

+            return Ok(vec![]);

+        }

+

+        let mut args = vec![self.parse_equation()?];

+

+        while self.current == Token::Comma {

+            self.advance()?;

+            args.push(self.parse_equation()?);

+        }

+

+        Ok(args)

+    }

+

+    fn parse_matrix(&mut self) -> Result<Expr> {

+        let mut rows = vec![self.parse_row()?];

+

+        while self.current == Token::Semicolon {

+            self.advance()?;

+            rows.push(self.parse_row()?);

+        }

+

+        // Single row = vector, multiple rows = matrix

+        if rows.len() == 1 {

+            Ok(Expr::Vector(rows.into_iter().next().unwrap()))

+        } else {

+            Ok(Expr::Matrix(rows))

+        }

+    }

+

+    fn parse_row(&mut self) -> Result<Vec<Expr>> {

+        let mut elems = vec![self.parse_equation()?];

+

+        while self.current == Token::Comma {

+            self.advance()?;

+            elems.push(self.parse_equation()?);

+        }

+

+        Ok(elems)

+    }

+