diff --git a/src/math.rs b/src/math.rs deleted file mode 100644 index 7e39e5d..0000000 --- a/src/math.rs +++ /dev/null @@ -1,903 +0,0 @@ -pub type Scalar = f64; - -pub type Vec2 = nalgebra::Vector2; -pub type Point2 = nalgebra::Point2; - -pub type Rot2 = nalgebra::UnitComplex; - -pub trait Region { - fn full() -> Self; - fn singleton(value: T) -> Self; - - fn nearest(&self, value: &T) -> Option; - fn contains(&self, value: &T) -> bool; -} - -#[derive(Clone, Debug)] -pub enum Region1 { - Empty, - Singleton(Scalar), - Range(Scalar, Scalar), - Union(Box, Box), - Full, -} - -impl Region for Region1 { - fn full() -> Self { - Region1::Full - } - - fn singleton(value: Scalar) -> Self { - Region1::Singleton(value) - } - - fn contains(&self, n: &Scalar) -> bool { - use Region1::*; - match self { - Empty => false, - Singleton(n1) => relative_eq!(n1, n), - Range(l, u) => *l <= *n && *n <= *u, - Union(r1, r2) => r1.contains(n) || r2.contains(n), - Full => true, - } - } - - fn nearest(&self, s: &Scalar) -> Option { - use Region1::*; - match self { - Empty => None, - Full => Some(*s), - Singleton(n) => Some(*n), - Range(l, u) => match (l < s, s < u) { - (true, true) => Some(*s), - (true, false) => Some(*u), - (false, true) => Some(*l), - _ => None, - }, - Union(r1, r2) => { - let distance = |a: Scalar, b: Scalar| (a - b).abs(); - match (r1.nearest(s), r2.nearest(s)) { - (None, None) => None, - (Some(n), None) | (None, Some(n)) => Some(n), - (Some(n1), Some(n2)) => Some({ - if distance(*s, n1) <= distance(*s, n2) { - n1 - } else { - n2 - } - }), - } - } - } - } -} - -// line starting at start, point at angle dir, with range extent -// ie. start + (cos dir, sin dir) * t for t in extent -#[derive(Clone, Debug)] -pub struct Line2 { - start: Point2, - dir: Rot2, - extent: Region1, -} - -impl Line2 { - pub fn new(start: Point2, dir: Rot2, extent: Region1) -> Self { - Self { start, dir, extent } - } - - pub fn evaluate(&self, t: Scalar) -> Point2 { - self.start + self.dir * Vec2::new(t, 0.) - } - - pub fn nearest(&self, p: &Point2) -> Point2 { - // rotate angle 90 degrees - let perp_dir = self.dir * Rot2::from_cos_sin_unchecked(0., 1.); - let perp = Line2::new(*p, perp_dir, Region1::Full); - if let Region2::Singleton(np) = self.intersect(&perp) { - np - } else { - panic!("Line2::nearest not found!"); - } - } - - pub fn intersect(&self, other: &Line2) -> Region2 { - // if the two lines are parallel... - let dirs = self.dir / other.dir; - if relative_eq!(dirs.sin_angle(), 0.) { - let starts = self.dir.to_rotation_matrix().inverse() * (other.start - self.start); - return if relative_eq!(starts.y, 0.) { - // and they are colinear - Region2::Line(self.clone()) - } else { - // they are parallel and never intersect - Region2::Empty - }; - } - // TODO: respect extent - let (a, b) = (self, other); - let (a_0, a_v, b_0, b_v) = (a.start, a.dir, b.start, b.dir); - let (a_c, a_s, b_c, b_s) = ( - a_v.cos_angle(), - a_v.sin_angle(), - b_v.cos_angle(), - b_v.sin_angle(), - ); - let t_b = (a_0.x * a_s - a_0.y * a_c + a_0.x * a_s + b_0.y * a_c) / (a_s * b_c - a_c * b_s); - Region2::Singleton(b.evaluate(t_b)) - } -} - -#[derive(Clone, Debug)] -pub enum Region2 { - Empty, - // single point at 0 - Singleton(Point2), - Line(Line2), - #[allow(dead_code)] - Union(Box, Box), - Full, -} - -impl Region for Region2 { - fn full() -> Self { - Region2::Full - } - - fn singleton(value: Point2) -> Self { - Region2::Singleton(value) - } - - fn contains(&self, p: &Point2) -> bool { - self.nearest(p).map_or(false, |n| relative_eq!(n, p)) - } - - fn nearest(&self, p: &Point2) -> Option { - use Region2::*; - match self { - Empty => None, - Full => Some(*p), - Singleton(n) => Some(*n), - Line(line) => Some(line.nearest(p)), - Union(r1, r2) => { - use nalgebra::distance; - match (r1.nearest(p), r2.nearest(p)) { - (None, None) => None, - (Some(n), None) | (None, Some(n)) => Some(n), - (Some(n1), Some(n2)) => Some({ - if distance(p, &n1) <= distance(p, &n2) { - n1 - } else { - n2 - } - }), - } - } - } - } -} - -impl Region2 { - pub fn union(r1: Region2, r2: Region2) -> Region2 { - use Region2::*; - match (r1, r2) { - (Empty, r) | (r, Empty) => r, - (Full, _) | (_, Full) => Full, - (r1, r2) => Union(Box::new(r1), Box::new(r2)), - } - } - - pub fn intersect(&self, other: &Region2) -> Region2 { - use Region2::*; - match (self, other) { - (Empty, _) | (_, Empty) => Empty, - (Full, r) | (r, Full) => r.clone(), - (Singleton(n1), Singleton(n2)) => { - if n1 == n2 { - Singleton(*n1) - } else { - Empty - } - } - (Singleton(n), o) | (o, Singleton(n)) => { - if o.contains(n) { - Singleton(*n) - } else { - Empty - } - } - (Line(l1), Line(l2)) => l1.intersect(l2), - (Union(un1, un2), o) | (o, Union(un1, un2)) => { - Self::union(un1.intersect(o), un2.intersect(o)) - } - } - } -} - -pub mod solve { - - use std::collections::{BTreeMap, BTreeSet}; - use std::fmt; - use std::iter::FromIterator; - - use crate::math::Scalar; - - // an unknown variable with an id - #[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord)] - struct Unknown(i64); - - type UnknownSet = BTreeSet; - - trait Unknowns { - fn unknowns(&self) -> UnknownSet; - fn has_unknowns(&self) -> bool; - fn has_unknown(&self, u: Unknown) -> bool; - } - - impl Unknowns for Scalar { - fn unknowns(&self) -> UnknownSet { - UnknownSet::new() - } - fn has_unknowns(&self) -> bool { - false - } - fn has_unknown(&self, _: Unknown) -> bool { - false - } - } - - impl Unknowns for Unknown { - fn unknowns(&self) -> UnknownSet { - FromIterator::from_iter(Some(*self)) - } - fn has_unknowns(&self) -> bool { - true - } - fn has_unknown(&self, u: Unknown) -> bool { - *self == u - } - } - - impl fmt::Display for Unknown { - fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - write!(f, "u{}", self.0) - } - } - - #[derive(Clone, Debug, PartialEq)] - enum Expr { - Unkn(Unknown), - Const(Scalar), - Sum(Exprs), - Neg(Box), - Product(Exprs), - Div(Box, Box), - } - - type Exprs = Vec; - - impl Unknowns for Exprs { - fn unknowns(&self) -> UnknownSet { - self.iter().flat_map(|e: &Expr| e.unknowns()).collect() - } - fn has_unknowns(&self) -> bool { - self.iter().any(|e: &Expr| e.has_unknowns()) - } - fn has_unknown(&self, u: Unknown) -> bool { - self.iter().any(|e: &Expr| e.has_unknown(u)) - } - } - - fn write_separated_exprs(es: &Exprs, f: &mut fmt::Formatter, sep: &str) -> fmt::Result { - let mut is_first = true; - for e in es { - if is_first { - is_first = false; - } else { - write!(f, "{}", sep)? - } - write!(f, "({})", e)? - } - Ok(()) - } - - fn remove_common_terms(l: &mut Vec, r: &mut Vec) -> Vec { - let common: Vec<_> = l.drain_filter(|e| r.contains(e)).collect(); - common.iter().for_each(|e| { - r.remove_item(e); - }); - common - } - - fn remove_term(terms: &mut Vec, term: &Expr) -> Option { - terms.remove_item(term) - } - - fn sum_fold(l: Expr, r: Expr) -> Expr { - use itertools::Itertools; - use Expr::*; - match (l, r) { - (Const(lc), Const(rc)) => Const(lc + rc), - (Const(c), o) | (o, Const(c)) if relative_eq!(c, 0.) => o, - (Product(mut l), Product(mut r)) => { - let comm = remove_common_terms(&mut l, &mut r); - Expr::new_product(Sum(comm), Expr::new_sum(Product(l), Product(r))).simplify() - } - (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.))).simplify() - } - None => Expr::new_sum(Product(l), r), - } - } - (l, r) => Expr::new_sum(l, r), - } - } - - fn group_sum(es: Exprs) -> Exprs { - use Expr::*; - let mut common: BTreeMap = BTreeMap::new(); - for e in es { - let unkns = e.unknowns(); - match common.get_mut(&unkns) { - None => { - match e { - Const(c) if relative_eq!(c, 0.) => (), - e => { - common.insert(unkns, e); - } - }; - } - Some(existing) => { - match existing { - Sum(ref mut es) => { - // already failed at merging, so just add it to the list - es.push(e); - } - other => { - *other = sum_fold(other.clone(), e); - } - }; - } - }; - } - common.into_iter().map(|(_, v)| v).collect() - } - - fn product_fold(l: Expr, r: Expr) -> Expr { - use itertools::Itertools; - use Expr::*; - match (l, r) { - (Const(lc), Const(rc)) => Const(lc * rc), - (Const(c), o) | (o, Const(c)) if relative_eq!(c, 1.) => o, - (Div(num, den), mul) | (mul, Div(num, den)) => { - if mul == *den { - *num - } else { - Expr::Div(Box::new(Expr::Product(vec![*num, mul])), den).simplify() - } - } - (l, r) => Expr::new_product(l, r), - } - } - - fn group_product(es: Exprs) -> Exprs { - use Expr::*; - let mut common: BTreeMap = BTreeMap::new(); - for e in es { - let unkns = e.unknowns(); - match common.get_mut(&unkns) { - None => { - match e { - Const(c) if relative_eq!(c, 1.) => (), - e => { - common.insert(unkns, e); - } - }; - } - Some(existing) => { - match existing { - Sum(ref mut es) => { - // already failed at merging, so just add it to the list - es.push(e); - } - other => *other = product_fold(other.clone(), e), - }; - } - }; - } - common.into_iter().map(|(_, v)| v).collect() - } - - fn distribute_product_sums(mut es: Exprs) -> Expr { - trace!("distribute_product_sums: {}", Product(es.clone())); - use itertools::Itertools; - use Expr::*; - let sums = es - .drain_filter(|e| match e { - Sum(_) => true, - _ => false, - }) - .map(|e| { - trace!("sum in product: {}", e); - match e { - Sum(es) => es, - _ => unreachable!(), - } - }); - let products: Vec<_> = sums.multi_cartesian_product().collect(); - if products.is_empty() { - trace!("no sums to distribute"); - return Product(es); - } - let sums = products - .into_iter() - .map(|mut prod| { - prod.extend(es.clone()); - trace!("prod: {}", Product(prod.clone())); - Product(prod) - }) - .collect(); - Sum(sums) - } - - impl Unknowns for Expr { - fn unknowns(&self) -> UnknownSet { - use Expr::*; - match self { - Unkn(u) => u.unknowns(), - Const(_) => UnknownSet::default(), - Sum(es) | Product(es) => es.unknowns(), - Div(l, r) => l.unknowns().union(&r.unknowns()).cloned().collect(), - Neg(e) => e.unknowns(), - } - } - fn has_unknowns(&self) -> bool { - use Expr::*; - match self { - Unkn(u) => u.has_unknowns(), - Const(_) => false, - Sum(es) | Product(es) => es.has_unknowns(), - Div(l, r) => l.has_unknowns() || r.has_unknowns(), - Neg(e) => e.has_unknowns(), - } - } - fn has_unknown(&self, u: Unknown) -> bool { - use Expr::*; - match self { - Unkn(u1) => u1.has_unknown(u), - Const(_) => false, - Sum(es) | Product(es) => es.has_unknown(u), - Div(l, r) => l.has_unknown(u) || r.has_unknown(u), - Neg(e) => e.has_unknown(u), - } - } - } - - impl Expr { - fn new_sum(e1: Expr, e2: Expr) -> Expr { - Expr::Sum(vec![e1, e2]) - } - fn new_product(e1: Expr, e2: Expr) -> Expr { - Expr::Product(vec![e1, e2]) - } - fn new_neg(e1: Expr) -> Expr { - Expr::Neg(Box::new(e1)) - } - fn new_div(num: Expr, den: Expr) -> Expr { - Expr::Div(Box::new(num), Box::new(den)) - } - fn new_minus(e1: Expr, e2: Expr) -> Expr { - Expr::Sum(vec![e1, Expr::new_neg(e2)]) - } - fn new_inv(den: Expr) -> Expr { - Expr::new_div(Expr::Const(1.), den) - } - - fn is_zero(&self) -> bool { - use Expr::*; - match self { - Const(c) => relative_eq!(*c, 0.), - _ => false, - } - } - - fn is_one(&self) -> bool { - use Expr::*; - match self { - Const(c) => relative_eq!(*c, 1.), - _ => false, - } - } - - fn simplify(self) -> Expr { - use Expr::*; - match self { - Sum(es) => { - let mut new_es: Vec<_> = es - .into_iter() - .map(|e| e.simplify()) - .flat_map(|e| match e { - Sum(more_es) => more_es, - other => vec![other], - }) - .collect(); - let pre_new_es = new_es.clone(); - new_es = group_sum(new_es); - trace!( - "simplify sum {} => {}", - Sum(pre_new_es), - Sum(new_es.clone()) - ); - - 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.simplify()) - .flat_map(|e| match e { - Product(more_es) => more_es, - other => vec![other], - }) - .collect(); - let pre_new_es = new_es.clone(); - let new_es = group_product(new_es); - trace!( - "simplify product {} => {}", - Product(pre_new_es), - Product(new_es.clone()) - ); - match new_es.len() { - 0 => Const(1.), // none - 1 => new_es.into_iter().next().unwrap(), // one - _ => Product(new_es), // many - } - } - Neg(mut v) => { - *v = v.simplify(); - trace!("simplify neg {}", Neg(v.clone())); - match v { - box Const(c) => Const(-c), - box Neg(v) => *v, - e => Product(vec![Const(-1.), *e]), - } - } - Div(mut num, mut den) => { - *num = num.simplify(); - *den = den.simplify(); - trace!("simplify div {}", Div(num.clone(), den.clone())); - match (num, den) { - (box Const(num), box Const(den)) => Const(num / den), - (num, box Const(den)) => { - if relative_eq!(den, 1.) { - *num - } else { - Expr::new_product(*num, Const(1. / den)) - } - } - (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), - None => Expr::new_div(Product(es), den), - }, - (num, den) => { - if num == den { - Expr::Const(1.) - } else { - Div(num, den) - } - } - } - } - e => e, - } - } - - fn distribute(self) -> Expr { - use Expr::*; - trace!("distribute {}", self); - match self { - Sum(mut es) => { - for e in &mut es { - *e = e.clone().distribute(); - } - Sum(es) - } - Product(es) => distribute_product_sums(es), - Div(mut num, mut den) => { - *num = num.distribute(); - *den = den.distribute(); - match (num, den) { - (box Sum(es), box den) => Sum(es - .into_iter() - .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) => { - // *l = Neg(l.clone()).distribute(); - // *r = Neg(r.clone()).distribute(); - // Sum(l, r) - // } - // box Product(mut l, r) => { - // *l = Neg(l.clone()).distribute(); - // Product(l, r) - // } - box Neg(v) => v.distribute(), - box Div(mut num, mut den) => { - *num = Neg(num.clone()).distribute(); - *den = Neg(den.clone()).distribute(); - Div(num, den) - } - e => Neg(e), - }, - e => e, - } - } - } - - impl fmt::Display for Expr { - fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - use Expr::*; - match self { - Unkn(u) => write!(f, "{}", u), - Const(c) => write!(f, "{}", c), - Sum(es) => write_separated_exprs(es, f, " + "), - Product(es) => write_separated_exprs(es, f, " * "), - Div(num, den) => write!(f, "({}) / ({})", num, den), - Neg(e) => write!(f, "-({})", e), - } - } - } - - #[derive(Clone, Debug, PartialEq)] - struct Eqn(Expr, Expr); - - impl Unknowns for Eqn { - fn unknowns(&self) -> UnknownSet { - self.0 - .unknowns() - .union(&self.1.unknowns()) - .cloned() - .collect() - } - fn has_unknowns(&self) -> bool { - self.0.has_unknowns() || self.1.has_unknowns() - } - fn has_unknown(&self, u: Unknown) -> bool { - self.0.has_unknown(u) || self.1.has_unknown(u) - } - } - - fn ord_by_unkn(a: Expr, b: Expr, u: Unknown) -> Option<(Expr, Expr)> { - if a.has_unknown(u) { - Some((a, b)) - } else if b.has_unknown(u) { - Some((b, a)) - } else { - None - } - } - - impl Eqn { - fn simplify(self) -> Eqn { - Eqn(self.0.simplify(), self.1.simplify()) - } - - fn solve(&self, for_u: Unknown) -> Option { - use Expr::*; - if !self.has_unknown(for_u) { - return None; - } - let (l, r) = ( - self.0 - .clone() /*.distribute()*/ - .simplify(), - self.1 - .clone() /*.distribute()*/ - .simplify(), - ); - let (mut l, mut r) = ord_by_unkn(l, r, for_u)?; - loop { - trace!("solve: {} == {}", l, r); - let (new_l, new_r): (Expr, Expr) = match l { - Unkn(u) => return if u == for_u { Some(r.simplify()) } else { None }, - Sum(es) => { - let (us, not_us): (Vec<_>, Vec<_>) = - 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<_>) = - es.into_iter().partition(|e| e.has_unknown(for_u)); - if us.len() != 1 { - return None; - } - ( - Product(us).simplify(), - Expr::new_div(r, Product(not_us)).simplify(), - ) - } - 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, - _ => return None, - }; - l = new_l; - r = new_r; - } - } - } - - impl fmt::Display for Eqn { - fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - write!(f, "{} == {}", self.0, self.1) - } - } - - #[derive(Clone, Debug, PartialEq)] - struct Eqns(Vec); - - impl Unknowns for Eqns { - fn unknowns(&self) -> UnknownSet { - self.0.iter().flat_map(|eqn: &Eqn| eqn.unknowns()).collect() - } - fn has_unknowns(&self) -> bool { - self.0.iter().any(|eqn: &Eqn| eqn.has_unknowns()) - } - fn has_unknown(&self, u: Unknown) -> bool { - self.0.iter().any(|eqn: &Eqn| eqn.has_unknown(u)) - } - } - - #[cfg(test)] - mod tests { - use super::*; - - #[test] - fn test_unknowns() { - let u1 = Unknown(1); - let u2 = Unknown(2); - let u3 = Unknown(3); - assert!(u1.has_unknowns()); - assert!(u2.has_unknowns()); - assert!(u1.has_unknown(u1)); - assert!(!u1.has_unknown(u2)); - assert!(u1.unknowns().contains(&u1)); - assert!(!u2.unknowns().contains(&u1)); - - let e1 = Expr::new_minus(Expr::Unkn(u1), Expr::Unkn(u2)); - assert!(e1.has_unknowns()); - assert!(e1.has_unknown(u1)); - assert!(e1.has_unknown(u2)); - assert!(!e1.has_unknown(u3)); - assert!(e1.unknowns().len() == 2); - } - - fn const_expr(e: Expr) -> Option { - match e { - Expr::Const(c) => Some(c), - _ => None, - } - } - - #[test] - fn test_solve() { - use Expr::*; - let _ = env_logger::try_init(); - let u1 = Unknown(1); - let e1 = Unkn(u1); - let e2 = Const(1.); - - let eqn = Eqn(e1.clone(), e2.clone()); - assert_eq!(eqn.solve(u1), Some(Const(1.))); - let eqn = Eqn(e2.clone(), e1.clone()); - assert_eq!(eqn.solve(u1), Some(Const(1.))); - - let e3 = Expr::new_sum(Const(1.), Expr::new_sum(Const(1.), Const(2.))); - let eqn = Eqn(e1.clone(), e3.clone()); - assert_eq!(eqn.solve(u1), Some(Const(4.))); - let e3 = Expr::new_minus(Const(1.), Const(1.)); - let eqn = Eqn(e1.clone(), e3.clone()); - assert_eq!(eqn.solve(u1), Some(Const(0.))); - - let e1 = Expr::new_div(Const(2.), Expr::new_minus(Const(1.), Const(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(); - assert!(const_expr(e.clone()).is_some()); - assert!(relative_eq!(const_expr(e.clone()).unwrap(), 5. / 3.)); - - let e1 = Expr::new_product(Const(2.), Expr::new_minus(Const(1.), Const(4.))); - let e2 = Expr::new_minus(Expr::new_product(Unkn(u1), Const(2.)), Unkn(u1)); - info!( - "e1==e2: {}=={} => {}=={}", - e1, - e2, - e1.clone().simplify(), - e2.clone().simplify() - ); - let eqn = Eqn(e1, e2); - let e = eqn.solve(u1).unwrap(); - assert!(const_expr(e.clone()).is_some()); - assert!(relative_eq!(const_expr(e.clone()).unwrap(), -6.)); - - let e1 = Expr::new_product(Const(2.), Expr::new_minus(Const(1.), Const(4.))); - let e2 = Expr::new_div( - Expr::new_sum( - Expr::new_product(Unkn(u1), Const(2.)), - Expr::new_product(Unkn(u1), Unkn(u1)), - ), - Unkn(u1), - ); - info!( - "{}=={} distrib=> {}=={}", - e1, - e2, - e1.clone().distribute(), - e2.clone().distribute() - ); - info!( - "{}=={} simplify=> {}=={}", - e1, - e2, - e1.clone().distribute().simplify(), - e2.clone().distribute().simplify() - ); - let eqn = Eqn(e1, e2); - let e = eqn.solve(u1).unwrap(); - assert!(const_expr(e.clone()).is_some()); - assert!(relative_eq!(const_expr(e.clone()).unwrap(), -8.)); - - let e1 = Expr::new_product(Const(2.), Expr::new_minus(Const(1.), Const(4.))); - let e2 = Expr::new_div( - Expr::new_sum( - Expr::new_product(Unkn(u1), Const(2.)), - Expr::new_sum( - Expr::new_sum(Expr::new_product(Unkn(u1), Unkn(u1)), Unkn(u1)), - Expr::new_minus(Const(2.), Const(1. + 1.)), - ), - ), - Unkn(u1), - ); - info!( - "e1==e2: {}=={} => {}=={}", - e1, - e2, - e1.clone().distribute().simplify(), - e2.clone().distribute().simplify().simplify() - ); - let eqn = Eqn(e1, e2); - let e = eqn.solve(u1).unwrap(); - assert!(const_expr(e.clone()).is_some()); - assert!(relative_eq!(const_expr(e.clone()).unwrap(), -9.)); - } - - } -} diff --git a/src/math/mod.rs b/src/math/mod.rs new file mode 100644 index 0000000..88d0fb5 --- /dev/null +++ b/src/math/mod.rs @@ -0,0 +1,217 @@ +pub mod solver; + +pub type Scalar = f64; + +pub type Vec2 = nalgebra::Vector2; +pub type Point2 = nalgebra::Point2; + +pub type Rot2 = nalgebra::UnitComplex; + +pub trait Region { + fn full() -> Self; + fn singleton(value: T) -> Self; + + fn nearest(&self, value: &T) -> Option; + fn contains(&self, value: &T) -> bool; +} + +#[derive(Clone, Debug)] +pub enum Region1 { + Empty, + Singleton(Scalar), + Range(Scalar, Scalar), + Union(Box, Box), + Full, +} + +impl Region for Region1 { + fn full() -> Self { + Region1::Full + } + + fn singleton(value: Scalar) -> Self { + Region1::Singleton(value) + } + + fn contains(&self, n: &Scalar) -> bool { + use Region1::*; + match self { + Empty => false, + Singleton(n1) => relative_eq!(n1, n), + Range(l, u) => *l <= *n && *n <= *u, + Union(r1, r2) => r1.contains(n) || r2.contains(n), + Full => true, + } + } + + fn nearest(&self, s: &Scalar) -> Option { + use Region1::*; + match self { + Empty => None, + Full => Some(*s), + Singleton(n) => Some(*n), + Range(l, u) => match (l < s, s < u) { + (true, true) => Some(*s), + (true, false) => Some(*u), + (false, true) => Some(*l), + _ => None, + }, + Union(r1, r2) => { + let distance = |a: Scalar, b: Scalar| (a - b).abs(); + match (r1.nearest(s), r2.nearest(s)) { + (None, None) => None, + (Some(n), None) | (None, Some(n)) => Some(n), + (Some(n1), Some(n2)) => Some({ + if distance(*s, n1) <= distance(*s, n2) { + n1 + } else { + n2 + } + }), + } + } + } + } +} + +// line starting at start, point at angle dir, with range extent +// ie. start + (cos dir, sin dir) * t for t in extent +#[derive(Clone, Debug)] +pub struct Line2 { + start: Point2, + dir: Rot2, + extent: Region1, +} + +impl Line2 { + pub fn new(start: Point2, dir: Rot2, extent: Region1) -> Self { + Self { start, dir, extent } + } + + pub fn evaluate(&self, t: Scalar) -> Point2 { + self.start + self.dir * Vec2::new(t, 0.) + } + + pub fn nearest(&self, p: &Point2) -> Point2 { + // rotate angle 90 degrees + let perp_dir = self.dir * Rot2::from_cos_sin_unchecked(0., 1.); + let perp = Line2::new(*p, perp_dir, Region1::Full); + if let Region2::Singleton(np) = self.intersect(&perp) { + np + } else { + panic!("Line2::nearest not found!"); + } + } + + pub fn intersect(&self, other: &Line2) -> Region2 { + // if the two lines are parallel... + let dirs = self.dir / other.dir; + if relative_eq!(dirs.sin_angle(), 0.) { + let starts = self.dir.to_rotation_matrix().inverse() * (other.start - self.start); + return if relative_eq!(starts.y, 0.) { + // and they are colinear + Region2::Line(self.clone()) + } else { + // they are parallel and never intersect + Region2::Empty + }; + } + // TODO: respect extent + let (a, b) = (self, other); + let (a_0, a_v, b_0, b_v) = (a.start, a.dir, b.start, b.dir); + let (a_c, a_s, b_c, b_s) = ( + a_v.cos_angle(), + a_v.sin_angle(), + b_v.cos_angle(), + b_v.sin_angle(), + ); + let t_b = (a_0.x * a_s - a_0.y * a_c + a_0.x * a_s + b_0.y * a_c) / (a_s * b_c - a_c * b_s); + Region2::Singleton(b.evaluate(t_b)) + } +} + +#[derive(Clone, Debug)] +pub enum Region2 { + Empty, + // single point at 0 + Singleton(Point2), + Line(Line2), + #[allow(dead_code)] + Union(Box, Box), + Full, +} + +impl Region for Region2 { + fn full() -> Self { + Region2::Full + } + + fn singleton(value: Point2) -> Self { + Region2::Singleton(value) + } + + fn contains(&self, p: &Point2) -> bool { + self.nearest(p).map_or(false, |n| relative_eq!(n, p)) + } + + fn nearest(&self, p: &Point2) -> Option { + use Region2::*; + match self { + Empty => None, + Full => Some(*p), + Singleton(n) => Some(*n), + Line(line) => Some(line.nearest(p)), + Union(r1, r2) => { + use nalgebra::distance; + match (r1.nearest(p), r2.nearest(p)) { + (None, None) => None, + (Some(n), None) | (None, Some(n)) => Some(n), + (Some(n1), Some(n2)) => Some({ + if distance(p, &n1) <= distance(p, &n2) { + n1 + } else { + n2 + } + }), + } + } + } + } +} + +impl Region2 { + pub fn union(r1: Region2, r2: Region2) -> Region2 { + use Region2::*; + match (r1, r2) { + (Empty, r) | (r, Empty) => r, + (Full, _) | (_, Full) => Full, + (r1, r2) => Union(Box::new(r1), Box::new(r2)), + } + } + + pub fn intersect(&self, other: &Region2) -> Region2 { + use Region2::*; + match (self, other) { + (Empty, _) | (_, Empty) => Empty, + (Full, r) | (r, Full) => r.clone(), + (Singleton(n1), Singleton(n2)) => { + if n1 == n2 { + Singleton(*n1) + } else { + Empty + } + } + (Singleton(n), o) | (o, Singleton(n)) => { + if o.contains(n) { + Singleton(*n) + } else { + Empty + } + } + (Line(l1), Line(l2)) => l1.intersect(l2), + (Union(un1, un2), o) | (o, Union(un1, un2)) => { + Self::union(un1.intersect(o), un2.intersect(o)) + } + } + } +} diff --git a/src/math/solver.rs b/src/math/solver.rs new file mode 100644 index 0000000..f114a01 --- /dev/null +++ b/src/math/solver.rs @@ -0,0 +1,684 @@ +use std::collections::{BTreeMap, BTreeSet}; +use std::fmt; +use std::iter::FromIterator; + +use crate::math::Scalar; + +// an unknown variable with an id +#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord)] +struct Unknown(i64); + +type UnknownSet = BTreeSet; + +trait Unknowns { + fn unknowns(&self) -> UnknownSet; + fn has_unknowns(&self) -> bool; + fn has_unknown(&self, u: Unknown) -> bool; +} + +impl Unknowns for Scalar { + fn unknowns(&self) -> UnknownSet { + UnknownSet::new() + } + fn has_unknowns(&self) -> bool { + false + } + fn has_unknown(&self, _: Unknown) -> bool { + false + } +} + +impl Unknowns for Unknown { + fn unknowns(&self) -> UnknownSet { + FromIterator::from_iter(Some(*self)) + } + fn has_unknowns(&self) -> bool { + true + } + fn has_unknown(&self, u: Unknown) -> bool { + *self == u + } +} + +impl fmt::Display for Unknown { + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + write!(f, "u{}", self.0) + } +} + +#[derive(Clone, Debug, PartialEq)] +enum Expr { + Unkn(Unknown), + Const(Scalar), + Sum(Exprs), + Neg(Box), + Product(Exprs), + Div(Box, Box), +} + +type Exprs = Vec; + +impl Unknowns for Exprs { + fn unknowns(&self) -> UnknownSet { + self.iter().flat_map(|e: &Expr| e.unknowns()).collect() + } + fn has_unknowns(&self) -> bool { + self.iter().any(|e: &Expr| e.has_unknowns()) + } + fn has_unknown(&self, u: Unknown) -> bool { + self.iter().any(|e: &Expr| e.has_unknown(u)) + } +} + +fn write_separated_exprs(es: &Exprs, f: &mut fmt::Formatter, sep: &str) -> fmt::Result { + let mut is_first = true; + for e in es { + if is_first { + is_first = false; + } else { + write!(f, "{}", sep)? + } + write!(f, "({})", e)? + } + Ok(()) +} + +fn remove_common_terms(l: &mut Vec, r: &mut Vec) -> Vec { + let common: Vec<_> = l.drain_filter(|e| r.contains(e)).collect(); + common.iter().for_each(|e| { + r.remove_item(e); + }); + common +} + +fn remove_term(terms: &mut Vec, term: &Expr) -> Option { + terms.remove_item(term) +} + +fn sum_fold(l: Expr, r: Expr) -> Expr { + use itertools::Itertools; + use Expr::*; + match (l, r) { + (Const(lc), Const(rc)) => Const(lc + rc), + (Const(c), o) | (o, Const(c)) if relative_eq!(c, 0.) => o, + (Product(mut l), Product(mut r)) => { + let comm = remove_common_terms(&mut l, &mut r); + Expr::new_product(Sum(comm), Expr::new_sum(Product(l), Product(r))).simplify() + } + (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.))).simplify() + } + None => Expr::new_sum(Product(l), r), + } + } + (l, r) => Expr::new_sum(l, r), + } +} + +fn group_sum(es: Exprs) -> Exprs { + use Expr::*; + let mut common: BTreeMap = BTreeMap::new(); + for e in es { + let unkns = e.unknowns(); + match common.get_mut(&unkns) { + None => { + match e { + Const(c) if relative_eq!(c, 0.) => (), + e => { + common.insert(unkns, e); + } + }; + } + Some(existing) => { + match existing { + Sum(ref mut es) => { + // already failed at merging, so just add it to the list + es.push(e); + } + other => { + *other = sum_fold(other.clone(), e); + } + }; + } + }; + } + common.into_iter().map(|(_, v)| v).collect() +} + +fn product_fold(l: Expr, r: Expr) -> Expr { + use itertools::Itertools; + use Expr::*; + match (l, r) { + (Const(lc), Const(rc)) => Const(lc * rc), + (Const(c), o) | (o, Const(c)) if relative_eq!(c, 1.) => o, + (Div(num, den), mul) | (mul, Div(num, den)) => { + if mul == *den { + *num + } else { + Expr::Div(Box::new(Expr::Product(vec![*num, mul])), den).simplify() + } + } + (l, r) => Expr::new_product(l, r), + } +} + +fn group_product(es: Exprs) -> Exprs { + use Expr::*; + let mut common: BTreeMap = BTreeMap::new(); + for e in es { + let unkns = e.unknowns(); + match common.get_mut(&unkns) { + None => { + match e { + Const(c) if relative_eq!(c, 1.) => (), + e => { + common.insert(unkns, e); + } + }; + } + Some(existing) => { + match existing { + Sum(ref mut es) => { + // already failed at merging, so just add it to the list + es.push(e); + } + other => *other = product_fold(other.clone(), e), + }; + } + }; + } + common.into_iter().map(|(_, v)| v).collect() +} + +fn distribute_product_sums(mut es: Exprs) -> Expr { + trace!("distribute_product_sums: {}", Product(es.clone())); + use itertools::Itertools; + use Expr::*; + let sums = es + .drain_filter(|e| match e { + Sum(_) => true, + _ => false, + }) + .map(|e| { + trace!("sum in product: {}", e); + match e { + Sum(es) => es, + _ => unreachable!(), + } + }); + let products: Vec<_> = sums.multi_cartesian_product().collect(); + if products.is_empty() { + trace!("no sums to distribute"); + return Product(es); + } + let sums = products + .into_iter() + .map(|mut prod| { + prod.extend(es.clone()); + trace!("prod: {}", Product(prod.clone())); + Product(prod) + }) + .collect(); + Sum(sums) +} + +impl Unknowns for Expr { + fn unknowns(&self) -> UnknownSet { + use Expr::*; + match self { + Unkn(u) => u.unknowns(), + Const(_) => UnknownSet::default(), + Sum(es) | Product(es) => es.unknowns(), + Div(l, r) => l.unknowns().union(&r.unknowns()).cloned().collect(), + Neg(e) => e.unknowns(), + } + } + fn has_unknowns(&self) -> bool { + use Expr::*; + match self { + Unkn(u) => u.has_unknowns(), + Const(_) => false, + Sum(es) | Product(es) => es.has_unknowns(), + Div(l, r) => l.has_unknowns() || r.has_unknowns(), + Neg(e) => e.has_unknowns(), + } + } + fn has_unknown(&self, u: Unknown) -> bool { + use Expr::*; + match self { + Unkn(u1) => u1.has_unknown(u), + Const(_) => false, + Sum(es) | Product(es) => es.has_unknown(u), + Div(l, r) => l.has_unknown(u) || r.has_unknown(u), + Neg(e) => e.has_unknown(u), + } + } +} + +impl Expr { + fn new_sum(e1: Expr, e2: Expr) -> Expr { + Expr::Sum(vec![e1, e2]) + } + fn new_product(e1: Expr, e2: Expr) -> Expr { + Expr::Product(vec![e1, e2]) + } + fn new_neg(e1: Expr) -> Expr { + Expr::Neg(Box::new(e1)) + } + fn new_div(num: Expr, den: Expr) -> Expr { + Expr::Div(Box::new(num), Box::new(den)) + } + fn new_minus(e1: Expr, e2: Expr) -> Expr { + Expr::Sum(vec![e1, Expr::new_neg(e2)]) + } + fn new_inv(den: Expr) -> Expr { + Expr::new_div(Expr::Const(1.), den) + } + + fn is_zero(&self) -> bool { + use Expr::*; + match self { + Const(c) => relative_eq!(*c, 0.), + _ => false, + } + } + + fn is_one(&self) -> bool { + use Expr::*; + match self { + Const(c) => relative_eq!(*c, 1.), + _ => false, + } + } + + fn simplify(self) -> Expr { + use Expr::*; + match self { + Sum(es) => { + let mut new_es: Vec<_> = es + .into_iter() + .map(|e| e.simplify()) + .flat_map(|e| match e { + Sum(more_es) => more_es, + other => vec![other], + }) + .collect(); + let pre_new_es = new_es.clone(); + new_es = group_sum(new_es); + trace!( + "simplify sum {} => {}", + Sum(pre_new_es), + Sum(new_es.clone()) + ); + + 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.simplify()) + .flat_map(|e| match e { + Product(more_es) => more_es, + other => vec![other], + }) + .collect(); + let pre_new_es = new_es.clone(); + let new_es = group_product(new_es); + trace!( + "simplify product {} => {}", + Product(pre_new_es), + Product(new_es.clone()) + ); + match new_es.len() { + 0 => Const(1.), // none + 1 => new_es.into_iter().next().unwrap(), // one + _ => Product(new_es), // many + } + } + Neg(mut v) => { + *v = v.simplify(); + trace!("simplify neg {}", Neg(v.clone())); + match v { + box Const(c) => Const(-c), + box Neg(v) => *v, + e => Product(vec![Const(-1.), *e]), + } + } + Div(mut num, mut den) => { + *num = num.simplify(); + *den = den.simplify(); + trace!("simplify div {}", Div(num.clone(), den.clone())); + match (num, den) { + (box Const(num), box Const(den)) => Const(num / den), + (num, box Const(den)) => { + if relative_eq!(den, 1.) { + *num + } else { + Expr::new_product(*num, Const(1. / den)) + } + } + (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), + None => Expr::new_div(Product(es), den), + }, + (num, den) => { + if num == den { + Expr::Const(1.) + } else { + Div(num, den) + } + } + } + } + e => e, + } + } + + fn distribute(self) -> Expr { + use Expr::*; + trace!("distribute {}", self); + match self { + Sum(mut es) => { + for e in &mut es { + *e = e.clone().distribute(); + } + Sum(es) + } + Product(es) => distribute_product_sums(es), + Div(mut num, mut den) => { + *num = num.distribute(); + *den = den.distribute(); + match (num, den) { + (box Sum(es), box den) => Sum(es + .into_iter() + .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) => { + // *l = Neg(l.clone()).distribute(); + // *r = Neg(r.clone()).distribute(); + // Sum(l, r) + // } + // box Product(mut l, r) => { + // *l = Neg(l.clone()).distribute(); + // Product(l, r) + // } + box Neg(v) => v.distribute(), + box Div(mut num, mut den) => { + *num = Neg(num.clone()).distribute(); + *den = Neg(den.clone()).distribute(); + Div(num, den) + } + e => Neg(e), + }, + e => e, + } + } +} + +impl fmt::Display for Expr { + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + use Expr::*; + match self { + Unkn(u) => write!(f, "{}", u), + Const(c) => write!(f, "{}", c), + Sum(es) => write_separated_exprs(es, f, " + "), + Product(es) => write_separated_exprs(es, f, " * "), + Div(num, den) => write!(f, "({}) / ({})", num, den), + Neg(e) => write!(f, "-({})", e), + } + } +} + +#[derive(Clone, Debug, PartialEq)] +struct Eqn(Expr, Expr); + +impl Unknowns for Eqn { + fn unknowns(&self) -> UnknownSet { + self.0 + .unknowns() + .union(&self.1.unknowns()) + .cloned() + .collect() + } + fn has_unknowns(&self) -> bool { + self.0.has_unknowns() || self.1.has_unknowns() + } + fn has_unknown(&self, u: Unknown) -> bool { + self.0.has_unknown(u) || self.1.has_unknown(u) + } +} + +fn ord_by_unkn(a: Expr, b: Expr, u: Unknown) -> Option<(Expr, Expr)> { + if a.has_unknown(u) { + Some((a, b)) + } else if b.has_unknown(u) { + Some((b, a)) + } else { + None + } +} + +impl Eqn { + fn simplify(self) -> Eqn { + Eqn(self.0.simplify(), self.1.simplify()) + } + + fn solve(&self, for_u: Unknown) -> Option { + use Expr::*; + if !self.has_unknown(for_u) { + return None; + } + let (l, r) = ( + self.0 + .clone() /*.distribute()*/ + .simplify(), + self.1 + .clone() /*.distribute()*/ + .simplify(), + ); + let (mut l, mut r) = ord_by_unkn(l, r, for_u)?; + loop { + trace!("solve: {} == {}", l, r); + let (new_l, new_r): (Expr, Expr) = match l { + Unkn(u) => return if u == for_u { Some(r.simplify()) } else { None }, + Sum(es) => { + let (us, not_us): (Vec<_>, Vec<_>) = + 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<_>) = + es.into_iter().partition(|e| e.has_unknown(for_u)); + if us.len() != 1 { + return None; + } + ( + Product(us).simplify(), + Expr::new_div(r, Product(not_us)).simplify(), + ) + } + 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, + _ => return None, + }; + l = new_l; + r = new_r; + } + } +} + +impl fmt::Display for Eqn { + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + write!(f, "{} == {}", self.0, self.1) + } +} + +#[derive(Clone, Debug, PartialEq)] +struct Eqns(Vec); + +impl Unknowns for Eqns { + fn unknowns(&self) -> UnknownSet { + self.0.iter().flat_map(|eqn: &Eqn| eqn.unknowns()).collect() + } + fn has_unknowns(&self) -> bool { + self.0.iter().any(|eqn: &Eqn| eqn.has_unknowns()) + } + fn has_unknown(&self, u: Unknown) -> bool { + self.0.iter().any(|eqn: &Eqn| eqn.has_unknown(u)) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_unknowns() { + let u1 = Unknown(1); + let u2 = Unknown(2); + let u3 = Unknown(3); + assert!(u1.has_unknowns()); + assert!(u2.has_unknowns()); + assert!(u1.has_unknown(u1)); + assert!(!u1.has_unknown(u2)); + assert!(u1.unknowns().contains(&u1)); + assert!(!u2.unknowns().contains(&u1)); + + let e1 = Expr::new_minus(Expr::Unkn(u1), Expr::Unkn(u2)); + assert!(e1.has_unknowns()); + assert!(e1.has_unknown(u1)); + assert!(e1.has_unknown(u2)); + assert!(!e1.has_unknown(u3)); + assert!(e1.unknowns().len() == 2); + } + + fn const_expr(e: Expr) -> Option { + match e { + Expr::Const(c) => Some(c), + _ => None, + } + } + + #[test] + fn test_solve() { + use Expr::*; + let _ = env_logger::try_init(); + let u1 = Unknown(1); + let e1 = Unkn(u1); + let e2 = Const(1.); + + let eqn = Eqn(e1.clone(), e2.clone()); + assert_eq!(eqn.solve(u1), Some(Const(1.))); + let eqn = Eqn(e2.clone(), e1.clone()); + assert_eq!(eqn.solve(u1), Some(Const(1.))); + + let e3 = Expr::new_sum(Const(1.), Expr::new_sum(Const(1.), Const(2.))); + let eqn = Eqn(e1.clone(), e3.clone()); + assert_eq!(eqn.solve(u1), Some(Const(4.))); + let e3 = Expr::new_minus(Const(1.), Const(1.)); + let eqn = Eqn(e1.clone(), e3.clone()); + assert_eq!(eqn.solve(u1), Some(Const(0.))); + + let e1 = Expr::new_div(Const(2.), Expr::new_minus(Const(1.), Const(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(); + assert!(const_expr(e.clone()).is_some()); + assert!(relative_eq!(const_expr(e.clone()).unwrap(), 5. / 3.)); + + let e1 = Expr::new_product(Const(2.), Expr::new_minus(Const(1.), Const(4.))); + let e2 = Expr::new_minus(Expr::new_product(Unkn(u1), Const(2.)), Unkn(u1)); + info!( + "e1==e2: {}=={} => {}=={}", + e1, + e2, + e1.clone().simplify(), + e2.clone().simplify() + ); + let eqn = Eqn(e1, e2); + let e = eqn.solve(u1).unwrap(); + assert!(const_expr(e.clone()).is_some()); + assert!(relative_eq!(const_expr(e.clone()).unwrap(), -6.)); + + let e1 = Expr::new_product(Const(2.), Expr::new_minus(Const(1.), Const(4.))); + let e2 = Expr::new_div( + Expr::new_sum( + Expr::new_product(Unkn(u1), Const(2.)), + Expr::new_product(Unkn(u1), Unkn(u1)), + ), + Unkn(u1), + ); + info!( + "{}=={} distrib=> {}=={}", + e1, + e2, + e1.clone().distribute(), + e2.clone().distribute() + ); + info!( + "{}=={} simplify=> {}=={}", + e1, + e2, + e1.clone().distribute().simplify(), + e2.clone().distribute().simplify() + ); + let eqn = Eqn(e1, e2); + let e = eqn.solve(u1).unwrap(); + assert!(const_expr(e.clone()).is_some()); + assert!(relative_eq!(const_expr(e.clone()).unwrap(), -8.)); + + let e1 = Expr::new_product(Const(2.), Expr::new_minus(Const(1.), Const(4.))); + let e2 = Expr::new_div( + Expr::new_sum( + Expr::new_product(Unkn(u1), Const(2.)), + Expr::new_sum( + Expr::new_sum(Expr::new_product(Unkn(u1), Unkn(u1)), Unkn(u1)), + Expr::new_minus(Const(2.), Const(1. + 1.)), + ), + ), + Unkn(u1), + ); + info!( + "e1==e2: {}=={} => {}=={}", + e1, + e2, + e1.clone().distribute().simplify(), + e2.clone().distribute().simplify().simplify() + ); + let eqn = Eqn(e1, e2); + let e = eqn.solve(u1).unwrap(); + assert!(const_expr(e.clone()).is_some()); + assert!(relative_eq!(const_expr(e.clone()).unwrap(), -9.)); + } + +}