got more stuff to work?

src/entity.rs

@@ -44,7 +44,7 @@ impl<T: Clone, TRegion: Region<T> + Clone> Constrainable<T, TRegion> {
     }
 
     pub fn resolve_with(&mut self, eqns: &Eqns) -> bool {
-        let resolved_constraints = self.constraints.clone().evaluate_with(eqns).simplify();
+        let resolved_constraints = self.constraints.clone().substitute(eqns).simplify();
         if let Some(n) = resolved_constraints.nearest(&self.value) {
             self.value = n;
             true
@@ -54,7 +54,7 @@ impl<T: Clone, TRegion: Region<T> + Clone> Constrainable<T, TRegion> {
     }
 
     pub fn evaluate_with(&mut self, eqns: &Eqns) -> bool {
-        self.constraints = self.constraints.clone().evaluate_with(eqns).simplify();
+        self.constraints = self.constraints.clone().substitute(eqns).simplify();
         if let Some(n) = self.constraints.nearest(&self.value) {
             self.value = n;
             true

src/main.rs

@@ -103,7 +103,7 @@ fn main() {
     }
     let e1 = Eqn::new(Expr::Unkn(u1), Expr::Const(1.));
     let e2 = Eqn::new(Expr::Unkn(u2), Expr::Const(1.));
-    let eqns = Eqns(vec![e1, e2]);
+    let eqns = Eqns(vec![e1, e2].into());
     for p in &mut points {
         p.borrow_mut().resolve_with(&eqns);
     }

src/math/eqn.rs

@@ -1,8 +1,10 @@
+use std::borrow::{Borrow, Cow};
+use std::collections::btree_map::Entry;
+use std::collections::BTreeMap;
 use std::fmt;
 
 use super::expr::Expr;
 use super::unknown::*;
-use super::Scalar;
 use crate::math::expr::Expr::*;
 
 #[derive(Clone, Debug, PartialEq)]
@@ -43,62 +45,71 @@ impl Eqn {
         Eqn(self.0.simplify(), self.1.simplify())
     }
 
-    pub fn solve(&self, for_u: Unknown) -> Option<Expr> {
+    pub fn substitute(self, eqns: &Eqns) -> Eqn {
-        use Expr::*;
+        Self(self.0.substitute(eqns), self.1.substitute(eqns))
-        if !self.has_unknown(for_u) {
+    }
-            return None;
+}
-        }
+
-        trace!("solve: {}", self);
+pub fn solve_eqn(eqn: &Eqn, for_u: Unknown) -> Option<Expr> {
-        let (l, r) = (
+    use Expr::*;
-            self.0
+    if !eqn.has_unknown(for_u) {
-                .clone().distribute()
+        return None;
-                .simplify(),
+    }
-            self.1
+    trace!("solve: {}", eqn);
-                .clone().distribute()
+    let (mut l, mut r) = (eqn.0.clone().distribute().simplify(), eqn.1.clone().distribute().simplify());
-                .simplify(),
+    if l == r {
-        );
+        return None
-        let (mut l, mut r) = ord_by_unkn(l, r, for_u)?;
+    };
-        loop {
+    l = (l - r).distribute().simplify();
-            trace!("solve iter: {} == {}", l, r);
+    r = Expr::Const(0.);
-            let (new_l, new_r): (Expr, Expr) = match l {
+    loop {
-                Unkn(u) => return if u == for_u { Some(r.simplify()) } else { None },
+        trace!("solve iter: {} == {}", l, r);
-                Sum(es) => {
+        let (new_l, new_r): (Expr, Expr) = match l {
-                    let (us, not_us): (Vec<_>, Vec<_>) =
+            Unkn(u) => return if u == for_u { Some(r.simplify()) } else { None },
-                        es.into_iter().partition(|e| e.has_unknown(for_u));
+            Sum(es) => {
-                    if us.len() != 1 {
+                let (us, not_us): (Vec<_>, Vec<_>) =
-                        return None;
+                    es.into_iter().partition(|e| e.has_unknown(for_u));
-                    }
+                if us.len() != 1 {
-                    (
+                    return None;
-                        Sum(us).simplify(),
-                        Expr::new_minus(r, Sum(not_us)).simplify(),
-                    )
                 }
-                Product(es) => {
+                (
-                    let (us, not_us): (Vec<_>, Vec<_>) =
+                    us.into_iter().next().unwrap(),
-                        es.into_iter().partition(|e| e.has_unknown(for_u));
+                    if not_us.len() == 0 {
-                    if us.len() != 1 {
+                        r
-                        return None;
+                    } else {
-                    }
+                        Expr::new_minus(r, Sum(not_us))
-                    (
+                    },
-                        Product(us).simplify(),
+                )
-                        Expr::new_div(r, Product(not_us)).simplify(),
+            }
-                    )
+            Product(es) => {
+                let (us, not_us): (Vec<_>, Vec<_>) =
+                    es.into_iter().partition(|e| e.has_unknown(for_u));
+                if us.len() != 1 {
+                    return None;
                 }
-                Neg(v) => (*v, Expr::new_neg(r)),
+                (
-                Div(num, den) => {
+                    us.into_iter().next().unwrap(),
-                    let (nu, du) = (num.has_unknown(for_u), den.has_unknown(for_u));
+                    if not_us.len() == 0 {
-                    match (nu, du) {
+                        r
-                        (true, false) => (*num, Expr::new_product(r, *den)),
+                    } else {
-                        (false, true) => (Expr::new_product(r, *den), *num),
+                        Expr::new_div(r, Product(not_us))
-                        (true, true) => return None, // TODO: simplify
+                    },
-                        (false, false) => return None,
+                )
-                    }
+            }
+            Neg(v) => (*v, Expr::new_neg(r)),
+            Div(num, den) => {
+                let (nu, du) = (num.has_unknown(for_u), den.has_unknown(for_u));
+                match (nu, du) {
+                    (true, false) => (*num, Expr::new_product(r, *den)),
+                    (false, true) => (Expr::new_product(r, *den), *num),
+                    (true, true) => return None, // TODO: simplify
+                    (false, false) => return None,
                 }
-                Const(_) => return None,
+            }
-            };
+            Const(_) => return None,
-            l = new_l;
+        };
-            r = new_r;
+        l = new_l.distribute().simplify();
-        }
+        r = new_r.distribute().simplify();
     }
 }
 
@@ -109,9 +120,9 @@ impl fmt::Display for Eqn {
 }
 
 #[derive(Clone, Debug, PartialEq)]
-pub struct Eqns(pub Vec<Eqn>);
+pub struct Eqns<'a>(pub Cow<'a, [Eqn]>);
 
-impl Unknowns for Eqns {
+impl<'a> Unknowns for Eqns<'a> {
     fn unknowns(&self) -> UnknownSet {
         self.0.iter().flat_map(|eqn: &Eqn| eqn.unknowns()).collect()
     }
@@ -123,11 +134,11 @@ impl Unknowns for Eqns {
     }
 }
 
-impl fmt::Display for Eqns {
+impl<'a> fmt::Display for Eqns<'a> {
     fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
-        write!(f, "{{ ");
+        write!(f, "{{ ")?;
         let mut first = true;
-        for eq in &self.0 {
+        for eq in (&self.0 as &Borrow<[Eqn]>).borrow() {
             if first {
                 first = false;
             } else {
@@ -139,8 +150,91 @@ impl fmt::Display for Eqns {
     }
 }
 
+pub fn solve_eqns<'a>(eqns: &Eqns<'a>, for_vars: &[Unknown]) -> Option<Eqns<'static>> {
+    let eqn_unknowns = eqns.unknowns();
+    let has_all_vars = for_vars.iter().all(|u| eqn_unknowns.contains(u));
+    if !has_all_vars {
+        // return None;
+    }
+    let mut var_sols = BTreeMap::<Unknown, Vec<Expr>>::new();
+    for u in &eqn_unknowns {
+        for eq in (&eqns.0 as &Borrow<[Eqn]>).borrow() {
+            let sol = solve_eqn(eq, *u);
+            if let Some(sol) = sol {
+                //println!("sol {0} -> {2} from {1}", u, eq, sol);
+                match var_sols.entry(*u) {
+                    Entry::Vacant(vac) => {
+                        vac.insert(vec![sol]);
+                    }
+                    Entry::Occupied(mut occ) => {
+                        if !occ.get().contains(&sol) {
+                            occ.get_mut().push(sol);
+                        };
+                    }
+                };
+            } else {
+                //println!("no sol for {0} in {1}", u, eq);
+            }
+        }
+    }
+    let mut new_eqns: Vec<Eqn> = Vec::new();
+    let mut all_have_sols = true;
+    for e in &var_sols {
+        let mut has_sol = false;
+        let mut last_expr = Unkn(*e.0);
+        println!("{}", e.0);
+        for s in e.1 {
+            println!(" == {}", s);
+            new_eqns.push(Eqn(last_expr, s.clone()));
+            last_expr = s.clone();
+            match s {
+                Const(_) => { has_sol = true; }
+                _ => {}
+            };
+        }
+        if !has_sol {
+            all_have_sols = false
+        }
+    }
+    let mut new_eqns = Eqns(new_eqns.into());
+    println!("new eqns: {}", new_eqns);
+    if all_have_sols {
+        return Some(new_eqns);
+    }
+    return solve_eqns(&new_eqns, for_vars);
+    /*
+    let mut more_eqns: Vec<Eqn> = Vec::new();
+    let mut new_orig_eqns = eqns.clone();
+    for eq in new_orig_eqns.0.to_mut() {
+        let u = match eq.0 {
+            Unkn(u) => u,
+            _ => continue,
+        };
+        let my_new_eqns = Eqns(
+            new_eqns
+                .0
+                .iter()
+                .cloned()
+                .filter(|e| !e.1.has_unknown(u))
+                .collect::<Vec<_>>()
+                .into(),
+        );
+        let r =
+            eq.1.clone()
+                .substitute(&my_new_eqns)
+                .distribute()
+                .simplify();
+        more_eqns.push(Eqn(Unkn(u), r));
+    }
+    let mut more_eqns = Eqns(more_eqns.into());
+    println!("more eqns: {}", more_eqns);
+    */
+    None
+}
+
 #[cfg(test)]
 mod tests {
+    use super::super::Scalar;
     use super::*;
 
     #[test]
@@ -171,7 +265,7 @@ mod tests {
     }
 
     #[test]
-    fn test_eqn_solve() {
+    fn test_solve_eqn() {
         use Expr::*;
         let _ = env_logger::try_init();
         let u1 = Unknown(1);
@@ -179,22 +273,22 @@ mod tests {
         let e2 = Const(1.);
 
         let eqn = Eqn(e1.clone(), e2.clone());
-        assert_eq!(eqn.solve(u1), Some(Const(1.)));
+        assert_eq!(solve_eqn(&eqn, u1), Some(Const(1.)));
         let eqn = Eqn(e2.clone(), e1.clone());
-        assert_eq!(eqn.solve(u1), Some(Const(1.)));
+        assert_eq!(solve_eqn(&eqn, u1), Some(Const(1.)));
 
         let e3: Expr = Expr::from(1.) + 1. + 2.;
         let eqn = Eqn(e1.clone(), e3.clone());
-        assert_eq!(eqn.solve(u1), Some(Const(4.)));
+        assert_eq!(solve_eqn(&eqn, u1), Some(Const(4.)));
         let e3 = Expr::from(1.) - 1.;
         let eqn = Eqn(e1.clone(), e3.clone());
-        assert_eq!(eqn.solve(u1), Some(Const(0.)));
+        assert_eq!(solve_eqn(&eqn, u1), Some(Const(0.)));
 
         let e1 = Expr::from(2.) / (Expr::from(1.) - 4.);
         let e2 = Expr::new_minus(Const(1.), Unkn(u1));
         let eqn = Eqn(e1, e2);
         info!("eqn: {} => {}", eqn, eqn.clone().simplify());
-        let e = eqn.solve(u1).unwrap();
+        let e = solve_eqn(&eqn, u1).unwrap();
         assert!(const_expr(e.clone()).is_some());
         assert!(relative_eq!(const_expr(e.clone()).unwrap(), 5. / 3.));
 
@@ -208,7 +302,7 @@ mod tests {
             e2.clone().simplify()
         );
         let eqn = Eqn(e1, e2);
-        let e = eqn.solve(u1).unwrap();
+        let e = solve_eqn(&eqn, u1).unwrap();
         assert!(const_expr(e.clone()).is_some());
         assert!(relative_eq!(const_expr(e.clone()).unwrap(), -6.));
 
@@ -229,7 +323,7 @@ mod tests {
             e2.clone().distribute().simplify()
         );
         let eqn = Eqn(e1, e2);
-        let e = eqn.solve(u1).unwrap();
+        let e = solve_eqn(&eqn, u1).unwrap();
         assert!(const_expr(e.clone()).is_some());
         assert!(relative_eq!(const_expr(e.clone()).unwrap(), -8.));
 
@@ -240,31 +334,82 @@ mod tests {
             e1,
             e2,
             e1.clone().distribute().simplify(),
-            e2.clone().distribute().simplify().simplify()
+            e2.clone().distribute().simplify()
         );
         let eqn = Eqn(e1, e2);
-        let e = eqn.solve(u1).unwrap();
+        let e = solve_eqn(&eqn, u1).unwrap();
         assert!(const_expr(e.clone()).is_some());
         assert!(relative_eq!(const_expr(e.clone()).unwrap(), -9.));
     }
 
     #[test]
-    fn test_eqns_solve() {
+    fn test_solve_eqn2() {
-       use Expr::*;
+        let _ = env_logger::try_init();
+
+        let u1 = Unknown(1);
+        let e1 = Expr::Unkn(u1);
+        let e2 = (u1 * 2.) - 2.;
+        info!(
+            "e1==e2: {}=={} => {}=={}",
+            e1,
+            e2,
+            e1.clone().distribute().simplify(),
+            e2.clone().distribute().simplify(),
+        );
+        let eqn = Eqn(e1, e2);
+        let e = solve_eqn(&eqn, u1).unwrap();
+        assert!(const_expr(e.clone()).is_some());
+        assert!(relative_eq!(const_expr(e.clone()).unwrap(), 2.));
+    }
+
+    #[test]
+    fn test_solve_equations() {
+        let _ = env_logger::try_init();
+
+        use Expr::*;
 
         let x = Unknown(1);
         let y = Unknown(2);
         let t1 = Unknown(3);
         let t2 = Unknown(4);
 
-        let eqns = Eqns(vec![
+        let eqns = Eqns(
-            Eqn::new(x.into(), t1 / 2.),
+            vec![
-            Eqn::new(y.into(), t1 / 2.),
+                Eqn::new(x.into(), t1 * 1.),
-            Eqn::new(x.into(), Const(1.0) - t2 / 2.),
+                Eqn::new(y.into(), t1 / 2.),
-            Eqn::new(y.into(), t2 / 2.),
+                Eqn::new(x.into(), Const(0.0) + t2 / 2.),
-        ]);
+                Eqn::new(y.into(), t2 / 2.),
+            ]
+            .into(),
+        );
         println!("eqns: {}", eqns);
-        let sol = eqns.solve(&[t1, t2]).unwrap();
+        let sol = solve_eqns(&eqns, &[t1, t2]).unwrap();
+        println!("sols: {}", sol);
+    }
+
+    #[test]
+    fn test_solve_equations2() {
+        let _ = env_logger::try_init();
+
+        use Expr::*;
+
+        let x = Unknown(1);
+        let y = Unknown(2);
+        let t1 = Unknown(3);
+        let t2 = Unknown(4);
+
+        let eqns = Eqns(
+            vec![
+                Eqn::new(x.into(), t1 * 1.),
+                Eqn::new(y.into(), t1 / 2.),
+                Eqn::new(x.into(), Const(1.0) + t2 / 1.),
+                Eqn::new(y.into(), t2 / 2.),
+            ]
+            .into(),
+        );
+        println!("eqns: {}", eqns);
+        let sol = solve_eqns(&eqns, &[t1, t2]).unwrap();
+        println!("sols: {}", sol);
     }
 
 }

src/math/expr.rs

@@ -1,9 +1,11 @@
+use std::borrow::Borrow;
 use std::collections::BTreeMap;
 use std::fmt;
 
-use super::eqn::Eqns;
+use super::eqn::{Eqn, Eqns};
 use super::unknown::{Unknown, UnknownSet, Unknowns};
 use super::Scalar;
+use std::f64::NAN;
 
 #[derive(Clone, Debug, PartialEq)]
 pub enum Expr {
@@ -29,9 +31,9 @@ impl Unknowns for Exprs {
     }
 }
 
-fn write_separated_exprs(es: &Exprs, f: &mut fmt::Formatter, sep: &str) -> fmt::Result {
+fn write_separated_exprs(es: &[Expr], f: &mut fmt::Formatter, sep: &str) -> fmt::Result {
-    let mut is_first = true;
     write!(f, "(")?;
+    let mut is_first = es.len() > 1;
     for e in es {
         if is_first {
             is_first = false;
@@ -66,16 +68,24 @@ fn sum_fold(l: Expr, r: Expr) -> Expr {
             if comm.is_empty() {
                 Expr::new_sum(Product(l), Product(r))
             } else {
-                Expr::new_product(Product(comm), Expr::new_sum(Product(l), Product(r)))
+                Expr::new_product(Product(comm), Expr::new_sum(Product(l), Product(r)).collapse())
             }
         }
         (Product(mut l), r) | (r, Product(mut l)) => {
             let comm = remove_term(&mut l, &r);
             match comm {
-                Some(_) => Expr::new_product(r, Expr::new_sum(Product(l), Const(1.))),
+                Some(_) => Expr::new_product(r, Expr::new_sum(Product(l), Const(1.))).collapse(),
                 None => Expr::new_sum(Product(l), r),
             }
         }
+        (Div(box ln, box ld), Div(box rn, box rd)) => {
+            if ld == rd {
+                Expr::new_div(Expr::new_sum(ln, rn), ld).collapse()
+            } else {
+                Expr::new_div(Expr::new_sum(Expr::new_product(ln, rd.clone()), Expr::new_product(rn, ld.clone())),
+                    Expr::new_product(ld, rd)).collapse()
+            }
+        }
         (l, r) => Expr::new_sum(l, r),
     }
 }
@@ -84,15 +94,16 @@ fn group_sum(es: Exprs) -> Exprs {
     use Expr::*;
     let mut common: BTreeMap<UnknownSet, Expr> = BTreeMap::new();
     for e in es {
+        match e {
+            Const(c) if relative_eq!(c, 0.) => {
+                continue;
+            }
+            _ => (),
+        };
         let unkns = e.unknowns();
         match common.get_mut(&unkns) {
             None => {
-                match e {
+                common.insert(unkns, e);
-                    Const(c) if relative_eq!(c, 0.) => (),
-                    e => {
-                        common.insert(unkns, e);
-                    }
-                };
             }
             Some(existing) => {
                 match existing {
@@ -110,7 +121,14 @@ fn group_sum(es: Exprs) -> Exprs {
     for c in common.values() {
         trace!("group sum value: {}", c);
     }
-    common.into_iter().map(|(_, v)| v).collect()
+    common
+        .into_iter()
+        .map(|(_, v)| v)
+        .filter(|e| match e {
+            Const(c) if relative_eq!(*c, 0.) => false,
+            _ => true,
+        })
+        .collect()
 }
 
 fn product_fold(l: Expr, r: Expr) -> Expr {
@@ -255,25 +273,25 @@ impl Expr {
         Expr::new_div(Expr::Const(1.), den)
     }
 
-    pub fn is_zero(self) -> bool {
+    pub fn is_zero(&self) -> bool {
         use Expr::*;
-        match self.simplify() {
+        match self.clone().simplify() {
             Const(c) => relative_eq!(c, 0.),
             _ => false,
         }
     }
 
-    pub fn is_one(self) -> bool {
+    pub fn is_one(&self) -> bool {
         use Expr::*;
-        match self.simplify() {
+        match self.clone().simplify() {
             Const(c) => relative_eq!(c, 1.),
             _ => false,
         }
     }
 
-    pub fn evaluate_with(self, eqns: &Eqns) -> Expr {
+    pub fn substitute<'a>(self, eqns: &Eqns<'a>) -> Expr {
         use Expr::*;
-        for eqn in &eqns.0 {
+        for eqn in (&eqns.0 as &Borrow<[Eqn]>).borrow() {
             if self == eqn.0 {
                 return eqn.1.clone();
             }
@@ -281,29 +299,129 @@ impl Expr {
         match self {
             Sum(mut es) => {
                 for e in &mut es {
-                    *e = e.clone().evaluate_with(eqns);
+                    *e = e.clone().substitute(eqns);
                 }
                 Sum(es)
             }
             Product(mut es) => {
                 for e in &mut es {
-                    *e = e.clone().evaluate_with(eqns);
+                    *e = e.clone().substitute(eqns);
                 }
                 Product(es)
             }
             Neg(mut e) => {
-                *e = e.evaluate_with(eqns);
+                *e = e.substitute(eqns);
                 Neg(e)
             }
             Div(mut num, mut den) => {
-                *num = num.evaluate_with(eqns);
+                *num = num.substitute(eqns);
-                *den = den.evaluate_with(eqns);
+                *den = den.substitute(eqns);
                 Div(num, den)
             }
             other => other,
         }
     }
 
+    pub fn collapse(self) -> Expr {
+        use Expr::*;
+        match self {
+            Sum(es) => {
+                let mut new_es: Vec<_> = es
+                    .into_iter()
+                    .map(|e| e.collapse())
+                    .flat_map(|e| match e {
+                        Sum(more_es) => more_es,
+                        other => vec![other],
+                    })
+                    .collect();
+
+                let consts = new_es.drain_filter(|e| match e {
+                    Const(_) => true,
+                    _ => false
+                }).fold(Const(0.), |l, r| match (l, r) {
+                    (Const(lc), Const(rc)) => Const(lc + rc),
+                    _ => unreachable!(),
+                });
+                new_es.push(consts);
+
+                match new_es.len() {
+                    0 => Const(0.),                          // none
+                    1 => new_es.into_iter().next().unwrap(), // one
+                    _ => Sum(new_es),                        // many
+                }
+            }
+            Product(es) => {
+                let new_es: Vec<_> = es
+                    .into_iter()
+                    .map(|e| e.collapse())
+                    .flat_map(|e| match e {
+                        Product(more_es) => more_es,
+                        other => vec![other],
+                    })
+                    .collect();
+                match new_es.len() {
+                    0 => Const(1.),                          // none
+                    1 => new_es.into_iter().next().unwrap(), // one
+                    _ => {
+                        if new_es.iter().any(|e| e.is_zero()) {
+                            Const(0.)
+                        } else {
+                            Product(new_es)
+                        }
+                    },                    // many
+                }
+            }
+            Neg(mut v) => {
+                *v = v.collapse();
+                match v {
+                    box Const(c) => Const(-c),
+                    box Neg(v) => *v,
+                    box Product(mut es) => {
+                        es.push(Const(-1.));
+                        Product(es).collapse()
+                    }
+                    e => Product(vec![Const(-1.), *e]).collapse(),
+                }
+            }
+            Div(mut num, mut den) => {
+                *num = num.collapse();
+                *den = den.collapse();
+                match (num, den) {
+                    (box Const(num), box Const(den)) => Const(num / den),
+                    (box Const(num), den) => {
+                        if relative_eq!(num, 0.) {
+                            Const(0.)
+                        } else {
+                            Div(Box::new(Const(num)), den)
+                        }
+                    }
+                    (num, box Const(den)) => {
+                        if relative_eq!(den, 1.) {
+                            *num
+                        } else {
+                            Expr::new_product(*num, Const(1. / den)).collapse()
+                        }
+                    }
+                    (num, box Div(dennum, denden)) => {
+                        Div(Box::new(Product(vec![*num, *denden])), dennum).collapse()
+                    }
+                    (box Product(mut es), box den) => match es.remove_item(&den) {
+                        Some(_) => Product(es),
+                        None => Expr::new_div(Product(es), den),
+                    },
+                    (num, den) => {
+                        if num == den {
+                            Expr::Const(1.)
+                        } else {
+                            Div(num, den)
+                        }
+                    }
+                }
+            }
+            e => e,
+        }
+    }
+
     pub fn simplify(self) -> Expr {
         use Expr::*;
         match self {
@@ -362,7 +480,7 @@ impl Expr {
                         es.push(Const(-1.));
                         Product(es).simplify()
                     }
-                    e => Product(vec![Const(-1.), *e]),
+                    e => Product(vec![Const(-1.), *e]).simplify(),
                 }
             }
             Div(mut num, mut den) => {
@@ -375,14 +493,14 @@ impl Expr {
                         if relative_eq!(den, 1.) {
                             *num
                         } else {
-                            Expr::new_product(*num, Const(1. / den))
+                            Expr::new_product(*num, Const(1. / den)).simplify()
                         }
                     }
                     (num, box Div(dennum, denden)) => {
                         Div(Box::new(Product(vec![*num, *denden])), dennum).simplify()
                     }
                     (box Product(mut es), box den) => match es.remove_item(&den) {
-                        Some(_) => Product(es),
+                        Some(_) => Product(es).simplify(),
                         None => Expr::new_div(Product(es), den),
                     },
                     (num, den) => {
@@ -394,6 +512,10 @@ impl Expr {
                     }
                 }
             }
+            Const(c) => if c.is_nan() {
+                println!("NAN!");
+                Const(c)
+            } else { Const(c) },
             e => e,
         }
     }
@@ -411,34 +533,34 @@ impl Expr {
                 res
             }
             Product(es) => distribute_product_sums(es),
-            Div(mut num, mut den) => {
+            Div(num, den) => match (num, den) {
-                *num = num.distribute();
+                (box Sum(es), box den) => Sum(es
-                *den = den.distribute();
+                    .into_iter()
-                match (num, den) {
+                    .map(|e| Expr::new_div(e, den.clone()).distribute())
-                    (box Sum(es), box den) => Sum(es
+                    .collect()),
-                        .into_iter()
+                (num, den) => Div(num, den),
-                        .map(|e| Expr::new_div(e, den.clone()))
+            },
-                        .collect()),
-                    (mut num, mut den) => Div(num, den),
-                }
-            }
             Neg(v) => match v {
-                // box Sum(mut l, mut r) => {
+                box Const(c) => Const(-c),
-                //     *l = Neg(l.clone()).distribute();
+                box Sum(mut es) => {
-                //     *r = Neg(r.clone()).distribute();
+                    for e in &mut es {
-                //     Sum(l, r)
+                        *e = Expr::new_neg(e.clone()).distribute();
-                // }
+                    }
-                // box Product(mut l, r) => {
+                    Sum(es)
-                //     *l = Neg(l.clone()).distribute();
+                }
-                //     Product(l, r)
+                box Product(mut es) => {
-                // }
+                    for e in &mut es {
+                        *e = e.clone().distribute();
+                    }
+                    es.push(Const(-1.));
+                    Product(es)
+                }
                 box Neg(v) => v.distribute(),
                 box Div(mut num, mut den) => {
-                    *num = Neg(num.clone()).distribute();
+                    *num = Neg(num.clone());
-                    *den = Neg(den.clone()).distribute();
+                    Div(num, den).distribute()
-                    Div(num, den)
                 }
-                e => Neg(e),
+                e => Neg(Box::new(e.distribute())),
             },
             e => e,
         }

src/math/region.rs

@@ -11,7 +11,7 @@ pub trait GenericRegion {
     fn full() -> Self;
     fn intersection(self, other: Self) -> Self;
     fn simplify(self) -> Self;
-    fn evaluate_with(self, eqns: &eqn::Eqns) -> Self;
+    fn substitute(self, eqns: &eqn::Eqns) -> Self;
 }
 
 pub trait Region<T>: GenericRegion {
@@ -63,12 +63,12 @@ impl GenericRegion for Region1 {
         }
     }
 
-    fn evaluate_with(self, eqns: &eqn::Eqns) -> Self {
+    fn substitute(self, eqns: &eqn::Eqns) -> Self {
         use Region1::*;
         match self {
-            Singleton(n) => Singleton(n.evaluate_with(eqns)),
+            Singleton(n) => Singleton(n.substitute(eqns)),
-            Range(l, u) => Range(l.evaluate_with(eqns), u.evaluate_with(eqns)),
+            Range(l, u) => Range(l.substitute(eqns), u.substitute(eqns)),
-            Intersection(r1, r2) => r1.evaluate_with(eqns).intersection(r2.evaluate_with(eqns)),
+            Intersection(r1, r2) => r1.substitute(eqns).intersection(r2.substitute(eqns)),
             other => other,
         }
     }
@@ -162,7 +162,7 @@ impl Line2 {
     }
 
     pub fn evaluate(&self, t: Value) -> Point2<Value> {
-        self.start.clone() + self.dir.clone() * t
+        self.start.clone() + self.dir * t
     }
 
     pub fn evaluate_extent(&self) -> Option<Point2<Value>> {
@@ -182,7 +182,7 @@ impl Line2 {
 
     pub fn nearest(&self, p: &Point2<Value>) -> Point2<Value> {
         // rotate angle 90 degrees
-        let perp_dir = self.dir.clone() + Rot2::cardinal(1);
+        let perp_dir = self.dir + Rot2::cardinal(1);
         let perp = Line2::new(p.clone(), perp_dir, Region1::Full);
         match self.intersect(&perp) {
             Region2::Singleton(np) => np,
@@ -193,7 +193,7 @@ impl Line2 {
 
     pub fn intersect(&self, other: &Line2) -> Region2 {
         // if the two lines are parallel...
-        let dirs = self.dir.clone() - other.dir.clone();
+        let dirs = self.dir - other.dir;
         if relative_eq!(dirs.sin(), 0.) {
             let starts = self.dir.conj() * (other.start.clone() - self.start.clone());
             return if starts.y.simplify().is_zero() {
@@ -208,16 +208,16 @@ impl Line2 {
         let (a, b) = (self, other);
         let (a_0, a_v, b_0, b_v) = (
             a.start.clone(),
-            a.dir.clone(),
+            a.dir,
             b.start.clone(),
-            b.dir.clone(),
+            b.dir,
         );
         let (a_c, a_s, b_c, b_s) = (a_v.cos(), a_v.sin(), b_v.cos(), b_v.sin());
-        let t_b = (a_0.x.clone() * a_s.clone()
+        let t_b = (a_0.x.clone() * a_s
-            - a_0.y.clone() * a_c.clone()
+            - a_0.y.clone() * a_c
-            - b_0.x.clone() * a_s.clone()
+            - b_0.x.clone() * a_s
-            + b_0.y.clone() * a_c.clone())
+            + b_0.y.clone() * a_c)
-            / (a_s.clone() * b_c.clone() - a_c.clone() * b_s.clone());
+            / (a_s * b_c - a_c * b_s);
         // Region2::Singleton(b.evaluate(t_b))
         trace!("intersect a: {}, b: {}, t_b = {}", a, b, t_b);
         let res = Region2::Line(b.clone().with_extent(Region1::Singleton(t_b.simplify())));
@@ -242,11 +242,11 @@ impl Line2 {
         Region2::Line(new_l)
     }
 
-    pub fn evaluate_with(self, eqns: &eqn::Eqns) -> Self {
+    pub fn substitute(self, eqns: &eqn::Eqns) -> Self {
         Line2 {
-            start: self.start.evaluate_with(eqns),
+            start: self.start.substitute(eqns),
             dir: self.dir,
-            extent: self.extent.evaluate_with(eqns),
+            extent: self.extent.substitute(eqns),
         }
     }
 }
@@ -300,12 +300,12 @@ impl GenericRegion for Region2 {
         }
     }
 
-    fn evaluate_with(self, eqns: &eqn::Eqns) -> Self {
+    fn substitute(self, eqns: &eqn::Eqns) -> Self {
         use Region2::*;
         match self {
-            Singleton(n) => Singleton(n.evaluate_with(eqns)),
+            Singleton(n) => Singleton(n.substitute(eqns)),
-            Line(l) => Line(l.evaluate_with(eqns)),
+            Line(l) => Line(l.substitute(eqns)),
-            Intersection(r1, r2) => r1.evaluate_with(eqns).intersection(r2.evaluate_with(eqns)),
+            Intersection(r1, r2) => r1.substitute(eqns).intersection(r2.substitute(eqns)),
             other => other,
         }
     }

src/math/vec.rs

@@ -98,10 +98,10 @@ impl Point2<Value> {
         }
     }
 
-    pub fn evaluate_with(self, eqns: &Eqns) -> Self {
+    pub fn substitute(self, eqns: &Eqns) -> Self {
         Self {
-            x: self.x.evaluate_with(eqns),
+            x: self.x.substitute(eqns),
-            y: self.y.evaluate_with(eqns),
+            y: self.y.substitute(eqns),
         }
     }
 }
@@ -184,8 +184,8 @@ impl Rot2 {
 
     pub fn from_angle(angle: Scalar) -> Self {
         Self {
-            cos: angle.cos().into(),
+            cos: angle.cos(),
-            sin: angle.sin().into(),
+            sin: angle.sin(),
         }
     }
 
@@ -195,22 +195,10 @@ impl Rot2 {
 
     pub fn cardinal(index: i64) -> Self {
         match index % 4 {
-            0 => Rot2 {
+            0 => Rot2 { cos: 1., sin: 0. },
-                cos: (1.).into(),
+            1 => Rot2 { cos: 0., sin: 1. },
-                sin: (0.).into(),
+            2 => Rot2 { cos: -1., sin: 0. },
-            },
+            3 => Rot2 { cos: 0., sin: -1. },
-            1 => Rot2 {
-                cos: (0.).into(),
-                sin: (1.).into(),
-            },
-            2 => Rot2 {
-                cos: (-1.).into(),
-                sin: (0.).into(),
-            },
-            3 => Rot2 {
-                cos: (0.).into(),
-                sin: (-1.).into(),
-            },
             _ => unreachable!(),
         }
     }
@@ -266,7 +254,7 @@ impl ops::Add<Rot2> for Rot2 {
     type Output = Rot2;
     fn add(self, rhs: Rot2) -> Rot2 {
         Rot2 {
-            cos: self.cos.clone() * rhs.cos.clone() - self.sin.clone() * rhs.sin.clone(),
+            cos: self.cos * rhs.cos - self.sin * rhs.sin,
             sin: self.cos * rhs.sin + self.sin * rhs.cos,
         }
     }
@@ -276,7 +264,7 @@ impl ops::Sub<Rot2> for Rot2 {
     type Output = Rot2;
     fn sub(self, rhs: Rot2) -> Rot2 {
         Rot2 {
-            cos: self.cos.clone() * rhs.cos.clone() + self.sin.clone() * rhs.sin.clone(),
+            cos: self.cos * rhs.cos + self.sin * rhs.sin,
             sin: self.sin * rhs.cos - self.cos * rhs.sin,
         }
     }

src/relation.rs

@@ -58,7 +58,7 @@ impl PointAngle {
         Self::new(
             p1,
             p2,
-            Rot2::from_cos_sin_unchecked((1.).into(), (0.).into()),
+            Rot2::from_cos_sin_unchecked(1., 0.),
         )
     }
 
@@ -66,7 +66,7 @@ impl PointAngle {
         Self::new(
             p1,
             p2,
-            Rot2::from_cos_sin_unchecked((0.).into(), (1.).into()),
+            Rot2::from_cos_sin_unchecked(0., 1.),
         )
     }
 }
@@ -76,7 +76,7 @@ impl Relation for PointAngle {
         use Region2::*;
         let (mut p1, mut p2) = (self.p1.borrow_mut(), self.p2.borrow_mut());
         let constrain_line = |p1: &Point2<Value>, p2: &mut PointEntity| {
-            let line = Region2::Line(Line2::new(p1.clone(), self.angle.clone(), Region1::Full));
+            let line = Region2::Line(Line2::new(p1.clone(), self.angle, Region1::Full));
             trace!(
                 "PointAngle line: {}, p2 constraint: {}",
                 line,
@@ -92,7 +92,7 @@ impl Relation for PointAngle {
             (Singleton(p1), Singleton(p2)) => {
                 // if the angle p1 and p2 form is parallel to self.angle, the result
                 // will have a y component of 0
-                let r = self.angle.clone().conj() * (p2.clone() - p1.clone());
+                let r = self.angle.conj() * (p2.clone() - p1.clone());
                 trace!("angle.cos: {}", r.x);
                 // if relative_eq!(r.y, 0.) {
                 ResolveResult::Constrained