#[macro_export]
macro_rules! format_eq {
($lhs:expr, [$($rhs:literal),* $(,)?]) => {{
format_eq!("{}", $lhs, [$($rhs),*]);
}};
($fmt:literal, $lhs:expr, [$($rhs:literal),* $(,)?]) => {{
#[cfg(feature = "pretty_assertions")]
use pretty_assertions::assert_eq;
let lhs = format!($fmt, $lhs);
let mut rhs = String::with_capacity(lhs.len());
$(
rhs.push_str($rhs);
rhs.push_str("\n");
)*
assert_eq!(lhs, rhs);
}};
}
use core::{
cmp::{Ordering, min},
fmt::{self, Debug, Display, Octal},
iter::{FromIterator, repeat_n},
mem::{swap, take},
num::{NonZero, NonZeroI32},
ops::{
Add, AddAssign, BitAnd, BitOr, BitXor, Div, DivAssign, Mul, MulAssign, Neg, Not, Rem, Shl,
Shr, Sub, SubAssign,
},
};
use std::{
collections::{BTreeMap, BTreeSet},
fmt::{Alignment, LowerHex},
};
pub trait Choose {
type Output;
#[must_use]
fn choose(self, other: Self) -> Self::Output;
}
impl Choose for u32 {
type Output = Self;
fn choose(self, other: Self) -> Self::Output {
let (n, k) = (self, other);
if n < k {
0
} else {
(0..k).fold(1, |r, i| r * (n - i) / (i + 1))
}
}
}
pub trait Factor
where
Self: Copy + Mul<Output = Self> + Div<Output = Self> + PartialEq + Eq + Default,
{
const ZERO: Self;
#[must_use]
fn gcd(self, other: Self) -> Self;
#[inline]
#[must_use]
fn lcm(self, other: Self) -> Self {
self.gcd_lcm(other).1
}
#[inline]
#[must_use]
fn gcd_lcm(self, other: Self) -> (Self, Self) {
if self == Self::ZERO && other == Self::ZERO {
(Self::ZERO, Self::ZERO)
} else {
let gcd = self.gcd(other);
let lcm = self * (other / gcd);
(gcd, lcm)
}
}
#[inline]
#[must_use]
fn gcd_bulk(r: impl IntoIterator<Item = Self>) -> Self {
r.into_iter().reduce(Self::gcd).unwrap_or_default()
}
#[inline]
#[must_use]
fn lcm_bulk(r: impl IntoIterator<Item = Self>) -> Self {
r.into_iter().reduce(Self::lcm).unwrap_or_default()
}
#[must_use]
fn signed(f: impl Fn(Self, Self) -> Self, r: impl IntoIterator<Item = Self>) -> Self
where
Self: Neg<Output = Self> + PartialOrd + Ord,
{
let (acc, neg, len) = r
.into_iter()
.fold((Self::ZERO, 0, 0), |(acc, neg, len), r| {
(
f(acc, r),
if r < Self::ZERO { neg + 1 } else { neg },
len + 1,
)
});
if neg > len / 2 { -acc } else { acc }
}
}
impl Factor for Rational {
const ZERO: Self = Self::ZERO;
#[inline]
fn gcd(self, other: Self) -> Self {
Self {
p: self.p.gcd(other.p),
q: self.q.lcm(other.q),
}
}
#[inline]
fn lcm(self, other: Self) -> Self {
Self {
p: self.p.lcm(other.p),
q: self.q.gcd(other.q),
}
}
#[inline]
fn gcd_lcm(self, other: Self) -> (Self, Self) {
let (g_p, l_p) = self.p.gcd_lcm(other.p);
let (g_q, l_q) = self.q.gcd_lcm(other.q);
(Self { p: g_p, q: l_q }, Self { p: l_p, q: g_q })
}
}
impl Factor for i32 {
const ZERO: Self = 0;
#[inline]
fn gcd(self, other: Self) -> Self {
self.unsigned_abs()
.gcd(other.unsigned_abs())
.try_into()
.expect("GCD is `i32::MIN.unsigned_abs() == i32::MAX.unsigned_abs() + 1`.")
}
}
impl Factor for u32 {
const ZERO: Self = 0;
#[allow(clippy::many_single_char_names, clippy::debug_assert_with_mut_call)]
fn gcd(self, other: Self) -> Self {
let mut a = self;
let mut b = other;
let g = if a == 0 || b == 0 {
a | b
} else {
let mut u = a.trailing_zeros();
let v = b.trailing_zeros();
let x = min(u, v);
b >>= v;
while a != 0 {
a >>= u;
let d = b.abs_diff(a);
u = d.trailing_zeros();
b = min(a, b);
a = d;
}
b << x
};
debug_assert_eq!(g, {
let mut a = self;
let mut b = other;
while b != 0 {
a %= b;
swap(&mut a, &mut b);
}
a
});
g
}
}
pub trait Algebra
where
Self: Copy
+ Clone
+ Eq
+ PartialEq
+ Ord
+ PartialOrd
+ Default
+ Into<Symbol>
+ TryFrom<Symbol, Error = Symbol>
+ Debug
+ Display
+ Mul<Output = (i8, Self)>
+ Not<Output = (i8, Self)>,
{
const N: u32;
#[must_use]
fn basis() -> impl ExactSizeIterator<Item = Self> + DoubleEndedIterator<Item = Self>;
#[must_use]
fn scalar() -> Self;
#[must_use]
fn pseudoscalar() -> Self;
#[must_use]
fn grade(&self) -> u32;
#[must_use]
fn blade_len(&self) -> usize;
#[must_use]
#[inline]
fn rev(self) -> (i8, Self) {
(1 - (self.grade().choose(2) & 1) as i8 * 2, self)
}
}
#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord, Default)]
pub struct Multivector<B: Algebra> {
pub map: BTreeMap<B, Polynomial>,
pub onc: bool,
}
impl<B: Algebra> TryFrom<Tree> for Multivector<B> {
type Error = Symbol;
fn try_from(tree: Tree) -> Result<Self, Self::Error> {
match tree {
Tree::Add(sib) => sib
.into_iter()
.try_fold(Self::zero(), |add, sib| Ok(add + Self::try_from(sib)?)),
Tree::Mul(sib) => sib
.into_iter()
.try_fold(Self::one(), |mul, sib| Ok(mul * Self::try_from(sib)?)),
Tree::Num(num) => Ok(num.into()),
Tree::Sym(sym) => sym.try_into(),
}
}
}
impl<B: Algebra> From<Rational> for Multivector<B> {
fn from(r: Rational) -> Self {
if r.is_zero() {
Self::zero()
} else {
let mut p = Polynomial::default();
p.map.insert(Monomial::one(), r);
let mut v = Self::default();
v.map.insert(B::scalar(), p);
v
}
}
}
impl<B: Algebra> TryFrom<Symbol> for Multivector<B> {
type Error = Symbol;
fn try_from(s: Symbol) -> Result<Self, Self::Error> {
if s.is_vec() {
Ok(Self::new([(Symbol::one(), s.try_into()?)]))
} else {
Ok(Self::new([(s, B::scalar())]))
}
}
}
impl<B: Algebra> Multivector<B> {
#[must_use]
#[inline]
pub fn new<E, S>(iter: E) -> Self
where
E: IntoIterator<Item = (S, B)>,
S: Into<Symbol>,
{
iter.into_iter().map(|(s, b)| ([[s]], b)).collect()
}
#[must_use]
pub fn alt(mut self) -> Self {
self.map.values_mut().for_each(|p| *p = take(p).alt());
self
}
#[must_use]
#[inline]
pub fn pin(self) -> Self {
self.cdm(Symbol::PIN)
}
#[must_use]
#[inline]
pub fn lhs(self) -> Self {
self.cdm(Symbol::LHS)
}
#[must_use]
#[inline]
pub fn rhs(self) -> Self {
self.cdm(Symbol::RHS)
}
#[must_use]
pub fn cdm(mut self, mark: char) -> Self {
self.map.values_mut().for_each(|p| *p = take(p).cdm(mark));
self
}
#[must_use]
pub fn swp(mut self) -> Self {
self.map.values_mut().for_each(|p| *p = take(p).swp());
self
}
#[must_use]
pub fn grades(&self) -> BTreeSet<u32> {
self.map.keys().map(Algebra::grade).collect()
}
#[must_use]
pub fn grade(&self) -> Option<u32> {
let grades = self.grades();
(grades.len() == 1)
.then(|| grades.first().copied())
.flatten()
}
#[must_use]
pub fn is_entity(&self) -> bool {
let mut set = BTreeSet::new();
for p in self.map.values() {
if let Some(m) = p.map.keys().next().filter(|_| p.map.len() == 1)
&& let Some(s) = m.map.keys().next().filter(|_| m.map.len() == 1)
&& set.insert(s)
{
continue;
}
return false;
}
true
}
#[must_use]
pub fn vectors(mut self) -> BTreeMap<u32, Self> {
let mut vectors = BTreeMap::new();
for grade in self.grades() {
let map = self
.map
.extract_if(.., |b, _p| b.grade() == grade)
.collect();
vectors.insert(grade, Self { map, onc: false });
}
vectors
}
#[must_use]
pub fn vector(mut self, grade: u32) -> Self {
self.map.retain(|b, _p| grade == b.grade());
self
}
#[must_use]
pub fn basis_blades(&self) -> BTreeSet<B> {
self.map.keys().copied().collect()
}
#[must_use]
pub fn rev(mut self) -> Self {
self.map.iter_mut().for_each(|(b, p)| {
let (s, _b) = b.rev();
if s < 0 {
*p = -take(p);
}
});
self
}
#[must_use]
#[inline]
pub fn zero() -> Self {
Self::default()
}
#[must_use]
pub fn one() -> Self {
Self::new([(Symbol::one(), B::scalar())])
}
#[must_use]
#[allow(clippy::missing_panics_doc)]
pub fn eval<M, S, R>(mut self, map: M) -> Self
where
M: IntoIterator<Item = (S, R)>,
S: Into<Symbol>,
R: Into<Rational>,
{
fn eval(vec: Vec<&mut Polynomial>, map: &BTreeMap<Symbol, Rational>) {
for old_p in vec {
let mut new_p = Polynomial::default();
for (old_m, old_r) in &old_p.map {
let mut new_m = old_m.clone();
let mut new_r = *old_r;
for (old_s, old_e) in &old_m.map {
if let Some(map_r) = map.get(old_s) {
new_m.map.remove(old_s);
new_r *= map_r.pow(old_e.get());
}
}
if new_m.map.is_empty() {
new_m = Monomial::one();
}
new_p.map.insert(new_m, new_r);
}
*old_p = new_p;
}
}
eval(
self.map.values_mut().collect(),
&map.into_iter().map(|(s, r)| (s.into(), r.into())).collect(),
);
self.map.retain(|_b, p| !p.map.is_empty());
self
}
#[must_use]
#[inline]
pub fn pol(self) -> Self {
self * !Self::one()
}
#[must_use]
pub fn norm_squared(self) -> Self {
self.clone() * self.rev()
}
#[must_use]
pub const fn unit(mut self) -> Self {
self.onc = true;
self
}
#[must_use]
pub fn cond(mut self, lhs: &Self, rhs: &Self) -> Self {
for (lhs_b, lhs_p) in lhs.map.clone() {
let rhs_p = rhs.map.get(&lhs_b).cloned().unwrap_or_default();
let lhs_g = lhs_p.signed_gcd();
let lhs_p = lhs_p / lhs_g;
for p in self.map.values_mut() {
let mut f = Factorization::from(p.clone());
for (p, _r) in f.map.values_mut() {
if p == &lhs_p {
*p = rhs_p.clone();
}
}
*p = f.into();
}
}
self.map.retain(|_b, p| !p.map.is_empty());
self
}
#[must_use]
pub fn ops(&self) -> (usize, usize) {
self.map.values().fold((0, 0), |(v_muls, v_adds), p| {
let f = Factorization::from(p.clone());
let (f_muls, f_adds) =
f.map
.into_iter()
.fold((0, 0), |(mut f_muls, f_adds), (m, (p, r))| {
if !m.is_one() {
f_muls += 1;
}
if !r.abs().is_one() {
f_muls += 1;
}
let (p_muls, p_adds) = p.ops();
(f_muls + p_muls, f_adds + p_adds)
});
(v_muls + f_muls, v_adds + f_adds)
})
}
}
impl<B: Algebra, P, M, S> FromIterator<(P, B)> for Multivector<B>
where
P: IntoIterator<Item = M>,
M: IntoIterator<Item = S>,
S: Into<Symbol>,
{
fn from_iter<V: IntoIterator<Item = (P, B)>>(iter: V) -> Self {
Self {
map: iter
.into_iter()
.map(|(p, b)| (b, Polynomial::from_iter(p)))
.collect(),
onc: false,
}
}
}
impl<B: Algebra> Add for Multivector<B> {
type Output = Self;
#[inline]
fn add(mut self, other: Self) -> Self::Output {
self += other;
self
}
}
impl<B: Algebra> AddAssign for Multivector<B> {
fn add_assign(&mut self, other: Self) {
for (b, p) in other.map {
*self.map.entry(b).or_default() += p;
}
self.map.retain(|_b, p| !p.map.is_empty());
}
}
impl<B: Algebra> Sub for Multivector<B> {
type Output = Self;
#[inline]
fn sub(mut self, other: Self) -> Self::Output {
self -= other;
self
}
}
impl<B: Algebra> SubAssign for Multivector<B> {
fn sub_assign(&mut self, other: Self) {
for (b, p) in other.map {
*self.map.entry(b).or_default() -= p;
}
self.map.retain(|_b, p| !p.map.is_empty());
}
}
impl<B: Algebra> Neg for Multivector<B> {
type Output = Self;
fn neg(mut self) -> Self::Output {
self.map.values_mut().for_each(|p| *p = -take(p));
self
}
}
impl<B: Algebra> Mul for Multivector<B> {
type Output = Self;
fn mul(self, other: Self) -> Self::Output {
let mut map = BTreeMap::<B, Polynomial>::new();
for (&lhs_b, lhs_p) in &self.map {
for (&rhs_b, rhs_p) in &other.map {
let (s, b) = lhs_b * rhs_b;
let p = lhs_p.clone() * rhs_p.clone();
*map.entry(b).or_default() += p * i32::from(s);
}
}
map.retain(|_b, p| !p.map.is_empty());
Self { map, onc: false }
}
}
impl<B: Algebra> MulAssign for Multivector<B> {
fn mul_assign(&mut self, other: Self) {
*self = take(self) * other;
}
}
impl<B: Algebra> Mul<i32> for Multivector<B> {
type Output = Self;
#[inline]
fn mul(mut self, other: i32) -> Self::Output {
self *= other;
self
}
}
impl<B: Algebra> MulAssign<i32> for Multivector<B> {
fn mul_assign(&mut self, other: i32) {
if other == 0 {
self.map = BTreeMap::default();
} else {
self.map.values_mut().for_each(|p| *p *= other);
}
}
}
impl<B: Algebra> Div<i32> for Multivector<B> {
type Output = Self;
#[inline]
fn div(mut self, other: i32) -> Self::Output {
self /= other;
self
}
}
impl<B: Algebra> DivAssign<i32> for Multivector<B> {
fn div_assign(&mut self, other: i32) {
assert_ne!(other, 0, "division by zero");
self.map.values_mut().for_each(|p| *p /= other);
}
}
impl<B: Algebra> BitOr for Multivector<B> {
type Output = Self;
fn bitor(self, other: Self) -> Self::Output {
let mut mv = Self::default();
for (lhs_grade, lhs_vector) in self.vectors() {
for (rhs_grade, rhs_vector) in other.clone().vectors() {
mv += (lhs_vector.clone() * rhs_vector).vector(rhs_grade.abs_diff(lhs_grade));
}
}
mv
}
}
impl<B: Algebra> BitXor for Multivector<B> {
type Output = Self;
fn bitxor(self, other: Self) -> Self::Output {
let mut mv = Self::default();
for (lhs_grade, lhs_vector) in self.vectors() {
for (rhs_grade, rhs_vector) in other.clone().vectors() {
mv += (lhs_vector.clone() * rhs_vector).vector(lhs_grade + rhs_grade);
}
}
mv
}
}
impl<B: Algebra> Not for Multivector<B> {
type Output = Self;
fn not(self) -> Self::Output {
let map = BTreeMap::new();
let map = self.map.into_iter().fold(map, |mut map, (b, p)| {
let (s, b) = !b;
map.insert(b, p * i32::from(s));
map
});
Self { map, onc: false }
}
}
impl<B: Algebra> BitAnd for Multivector<B> {
type Output = Self;
#[inline]
fn bitand(self, other: Self) -> Self::Output {
!!!(!self ^ !other)
}
}
impl<B: Algebra> Rem for Multivector<B> {
type Output = Self;
fn rem(self, other: Self) -> Self::Output {
(self.clone() * other.clone() - other * self) / 2
}
}
impl<B: Algebra> Shl for Multivector<B> {
type Output = Self;
fn shl(mut self, other: Self) -> Self::Output {
let onc = other.onc.then(|| {
let lhs = other.clone().norm_squared();
(self.clone() | (Self::one() - lhs.clone()), lhs, Self::one())
});
if self
.grade()
.zip(other.grade())
.map(|(lhs_grade, rhs_grade)| (lhs_grade * rhs_grade) % 2)
== Some(1)
{
self = -self;
}
let shl = other.clone() * self * other.rev();
if let Some((onc, lhs, rhs)) = onc {
(shl + onc).cond(&lhs, &rhs)
} else {
shl
}
}
}
impl<B: Algebra> Shr for Multivector<B> {
type Output = Self;
fn shr(self, other: Self) -> Self::Output {
let shr = (self | other.clone()) * other.clone().rev();
if other.onc {
shr.cond(&other.norm_squared(), &Self::one())
} else {
shr
}
}
}
impl<B: Algebra> Display for Multivector<B> {
#[allow(clippy::too_many_lines, clippy::cognitive_complexity)]
fn fmt(&self, fmt: &mut fmt::Formatter) -> fmt::Result {
fn traverse<'a>(
fmt: &mut fmt::Formatter,
tree: &Tree,
depth: usize,
grasp: bool,
close: bool,
mut defer: &'a str,
) -> Result<&'a str, fmt::Error> {
let math = fmt.fill() == '$';
let wide = matches!(fmt.align(), Some(Alignment::Center | Alignment::Right)) || math;
let code = fmt.align() == Some(Alignment::Center) || fmt.alternate() || math;
match tree {
Tree::Add(siblings) => {
let grasp = grasp || !(defer.is_empty() || defer == "+" || defer == " + ");
if grasp {
write!(fmt, "{defer}")?;
defer = "";
write!(fmt, "(")?;
}
for (index, sibling) in siblings.iter().enumerate() {
if fmt.align().is_none() || index != 0 {
defer = if wide { " + " } else { "+" };
}
let close = index + 1 != siblings.len();
defer = traverse(fmt, sibling, depth + 1, grasp, close, defer)?;
write!(fmt, "{defer}")?;
}
if grasp {
write!(fmt, ")")?;
}
defer = "";
}
Tree::Mul(siblings) => {
let is_one = siblings.last().and_then(Tree::as_sym).is_some_and(|sym| {
(sym.is_one() || sym.is_scalar()) && depth <= 1 && siblings.len() == 2
});
for (index, sibling) in siblings.iter().enumerate() {
defer = if index == 0 {
defer
} else if code && defer.is_empty() && !is_one {
if math {
" "
} else if wide {
" * "
} else {
"*"
}
} else {
""
};
let grasp = !(index == 0 && is_one);
defer = traverse(fmt, sibling, depth + 1, grasp, close, defer)?;
}
defer = "";
}
Tree::Num(num) => {
let suffix = if fmt.precision().is_some() { ".0" } else { "" };
if num.abs().is_one() {
if num.is_negative() {
if wide && !defer.is_empty() {
write!(fmt, " - ")?;
} else {
write!(fmt, "-")?;
}
} else if fmt.align().is_none() {
write!(fmt, "+")?;
}
defer = if fmt.precision().is_some() {
"1.0"
} else {
"1"
};
} else if num.is_negative() {
if wide && !defer.is_empty() {
write!(fmt, " - ")?;
} else {
write!(fmt, "-")?;
}
let num = num.abs();
write!(fmt, "{num}{suffix}")?;
defer = "";
} else if !num.is_zero() {
write!(fmt, "{defer}{num}{suffix}")?;
defer = "";
}
}
Tree::Sym(sym) => {
write!(fmt, "{defer}")?;
if sym.is_vec() && math && !fmt.sign_aware_zero_pad() {
if defer != " " {
write!(fmt, " ")?;
}
write!(fmt, "& ")?;
}
if !(sym.is_one() || sym.is_scalar()) || depth == 0 {
Display::fmt(sym, fmt)?;
defer = "";
}
if !fmt.sign_aware_zero_pad() && sym.is_vec() {
if close && math {
write!(fmt, " \\\\\n ")?;
if let Some(width) = fmt.width() {
write!(fmt, "{:width$}", "")?;
}
} else {
writeln!(fmt)?;
}
}
}
}
if depth == 0 {
if defer == "1" || defer == "1.0" {
write!(fmt, "{defer}")?;
}
defer = "";
}
Ok(defer)
}
let tree = if fmt.sign_plus() {
Tree::from(self.clone())
} else {
Tree::with_factorization(self.clone(), fmt.sign_minus())
};
let defer = if fmt.align().is_none() { "+" } else { "" };
let math = fmt.fill() == '$';
let wide = matches!(fmt.align(), Some(Alignment::Center | Alignment::Right));
if math && wide && !fmt.sign_aware_zero_pad() {
let align = if fmt.align() == Some(Alignment::Center) {
"[t]"
} else {
""
};
write!(fmt, "\\begin{{aligned}}{align}\n ")?;
if let Some(width) = fmt.width() {
write!(fmt, "{:width$}", "")?;
}
}
traverse(fmt, &tree, 0, false, false, defer)?;
if math && wide && !fmt.sign_aware_zero_pad() {
if let Some(width) = fmt.width() {
write!(fmt, "{:width$}", "")?;
}
writeln!(fmt, "\\end{{aligned}}")?;
}
Ok(())
}
}
impl<B: Algebra> LowerHex for Multivector<B> {
fn fmt(&self, fmt: &mut fmt::Formatter) -> fmt::Result {
for (b, p) in &self.map {
let mut map = BTreeMap::new();
map.insert(B::scalar(), p.clone());
let m = Self { map, onc: self.onc };
if let Some(width) = fmt.width() {
write!(fmt, "{:width$}", "")?;
}
let deref = fmt.align() == Some(Alignment::Center);
let rust = fmt.alternate() || deref;
if fmt.sign_plus() {
if rust {
if deref {
writeln!(fmt, "{b:^} = {m:^+0.1};")?;
} else {
writeln!(fmt, "let {b:#} = {m:>+#0.1};")?;
}
} else {
writeln!(fmt, "{b:#}={m:<+#0}")?;
}
} else if fmt.sign_minus() {
if rust {
if deref {
writeln!(fmt, "{b:^} = {m:^-0.1};")?;
} else {
writeln!(fmt, "let {b:#} = {m:>-#0.1};")?;
}
} else {
writeln!(fmt, "{b:#}={m:<-#0}")?;
}
} else if rust {
if deref {
writeln!(fmt, "{b:^} = {m:^0.1};")?;
} else {
writeln!(fmt, "let {b:#} = {m:>#0.1};")?;
}
} else {
writeln!(fmt, "{b:#}={m:<#0}")?;
}
}
Ok(())
}
}
impl<B: Algebra> Octal for Multivector<B> {
fn fmt(&self, fmt: &mut fmt::Formatter) -> fmt::Result {
fn traverse(
fmt: &mut fmt::Formatter,
tree: &Tree,
inode: usize,
mut index: usize,
) -> Result<usize, fmt::Error> {
let [add, mul] = if fmt.alternate() {
["+", "*"]
} else {
["∑", "∏"]
};
match tree {
Tree::Add(siblings) => {
writeln!(fmt, " n{index} [label=\"{add}\" shape=box]")?;
if inode != index {
writeln!(fmt, " n{inode} -> n{index}")?;
}
let inode = index;
for sibling in siblings {
index = traverse(fmt, sibling, inode, index + 1)?;
}
Ok(index)
}
Tree::Mul(siblings) => {
writeln!(fmt, " n{index} [label=\"{mul}\" shape=box]")?;
if inode != index {
writeln!(fmt, " n{inode} -> n{index}")?;
}
let inode = index;
for sibling in siblings {
index = traverse(fmt, sibling, inode, index + 1)?;
}
Ok(index)
}
Tree::Num(num) => {
writeln!(fmt, " n{index} [label=\"{num}\" shape=circle]")?;
if inode != index {
writeln!(fmt, " n{inode} -> n{index}")?;
}
Ok(index)
}
Tree::Sym(sym) => {
let shape = if sym.is_vec() {
"diamond"
} else {
match sym.cdm {
Symbol::NIL => "ellipse",
Symbol::PIN => "hexagon",
Symbol::LHS => "larrow",
Symbol::RHS => "rarrow",
label => panic!("unknown symbol label `{label}`"),
}
};
if fmt.alternate() {
writeln!(fmt, " n{index} [label=\"{sym:#}\" shape={shape}]")?;
} else {
writeln!(fmt, " n{index} [label=\"{sym}\" shape={shape}]")?;
}
if inode != index {
writeln!(fmt, " n{inode} -> n{index}")?;
}
Ok(index)
}
}
}
writeln!(fmt, "digraph vee {{")?;
if fmt.sign_aware_zero_pad() {
writeln!(fmt, " rankdir=LR")?;
}
let t = if fmt.sign_plus() {
Tree::from(self.clone())
} else {
Tree::with_factorization(self.clone(), fmt.sign_minus())
};
traverse(fmt, &t, 0, 0)?;
writeln!(fmt, "}}")?;
Ok(())
}
}
#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord)]
pub struct Factorization {
pub map: BTreeMap<Monomial, (Polynomial, Rational)>,
pub gcd: Rational,
}
impl Factorization {
#[must_use]
pub fn new(p: Polynomial, signed: bool) -> Self {
let mut f = p
.map
.into_iter()
.fold(Self::default(), |mut f, (mut m, c)| {
let mut g = Monomial {
map: m.map.extract_if(.., |s, _e| s.is_pin()).collect(),
};
if g.map.is_empty() {
g = Monomial::one();
}
if m.map.is_empty() {
m = Monomial::one();
}
let (p, r) = f.map.entry(g).or_default();
*p.map.entry(m).or_default() += c;
*r = Rational::ONE;
f
});
let gcd = if signed {
Polynomial::signed_gcd
} else {
Polynomial::gcd
};
f.map.values_mut().for_each(|(p, r)| {
*r = gcd(p);
*p /= *r;
});
f.gcd = if signed {
Rational::signed(Rational::gcd, f.map.values().map(|(_p, r)| *r))
} else {
Rational::gcd_bulk(f.map.values().map(|(_p, r)| *r))
};
f.map.values_mut().for_each(|(_p, r)| *r /= f.gcd);
f
}
#[must_use]
pub fn is_zero(&self) -> bool {
self.map.is_empty()
|| self.gcd.is_zero()
|| self
.map
.iter()
.all(|(m, (p, c))| m.is_zero() || p.is_zero() || c.is_zero())
}
#[must_use]
pub fn is_one(&self) -> bool {
let mut sum = 0;
self.gcd.is_one()
&& self
.map
.iter()
.filter(|(m, (p, c))| !(m.is_zero() || p.is_zero() || c.is_zero()))
.all(|(m, (p, c))| {
sum += 1;
m.is_one() && p.is_one() && c.is_one()
})
&& sum == 1
}
}
impl Default for Factorization {
fn default() -> Self {
Self {
map: BTreeMap::new(),
gcd: Rational::ONE,
}
}
}
impl From<Polynomial> for Factorization {
#[inline]
fn from(p: Polynomial) -> Self {
Self::new(p, true)
}
}
impl From<Factorization> for Polynomial {
fn from(f: Factorization) -> Self {
let mut p = Self::default();
for (f_m, (f_p, f_r)) in f.map {
for (m, c) in f_p.map {
*p.map.entry(f_m.clone() * m).or_default() += c * f_r * f.gcd;
}
}
p
}
}
#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord, Default)]
pub struct Polynomial {
pub map: BTreeMap<Monomial, Rational>,
}
impl Polynomial {
#[must_use]
pub fn one() -> Self {
let mut map = BTreeMap::new();
map.insert(Monomial::one(), Rational::ONE);
Self { map }
}
#[must_use]
pub fn is_zero(&self) -> bool {
self.map.is_empty() || self.map.iter().all(|(m, c)| m.is_zero() || c.is_zero())
}
#[must_use]
pub fn is_one(&self) -> bool {
self.map
.first_key_value()
.is_some_and(|(m, c)| self.map.len() == 1 && m.is_one() && c.is_one())
}
#[must_use]
pub fn alt(self) -> Self {
let map = BTreeMap::new();
let map = self.map.into_iter().fold(map, |mut map, (s, c)| {
map.insert(s.alt(), c);
map
});
Self { map }
}
#[must_use]
pub fn cdm(self, mark: char) -> Self {
let map = BTreeMap::new();
let map = self.map.into_iter().fold(map, |mut map, (s, c)| {
map.insert(s.cdm(mark), c);
map
});
Self { map }
}
#[must_use]
pub fn swp(self) -> Self {
let map = BTreeMap::new();
let map = self.map.into_iter().fold(map, |mut map, (s, c)| {
map.insert(s.swp(), c);
map
});
Self { map }
}
#[must_use]
pub fn ops(&self) -> (usize, usize) {
(
self.map
.iter()
.map(|(m, r)| {
(m.map
.values()
.filter_map(|e| usize::try_from(e.get()).ok())
.sum::<usize>()
+ usize::from(!r.abs().is_one()))
.saturating_sub(1)
})
.sum::<usize>(),
self.map.len().saturating_sub(1),
)
}
#[must_use]
pub fn gcd(&self) -> Rational {
Rational::gcd_bulk(self.map.values().copied())
}
#[must_use]
pub fn lcm(&self) -> Rational {
Rational::lcm_bulk(self.map.values().copied())
}
pub fn signed_gcd(&self) -> Rational {
Rational::signed(Rational::gcd, self.map.values().copied())
}
pub fn signed_lcd(&self) -> Rational {
Rational::signed(Rational::lcm, self.map.values().copied())
}
}
impl<M, S> FromIterator<M> for Polynomial
where
M: IntoIterator<Item = S>,
S: Into<Symbol>,
{
fn from_iter<P: IntoIterator<Item = M>>(iter: P) -> Self {
Self {
map: iter
.into_iter()
.map(|sym| (Monomial::from_iter(sym), Rational::ONE))
.collect(),
}
}
}
impl Add for Polynomial {
type Output = Self;
#[inline]
fn add(mut self, other: Self) -> Self::Output {
self += other;
self
}
}
impl AddAssign for Polynomial {
fn add_assign(&mut self, other: Self) {
for (m, c) in other.map {
*self.map.entry(m).or_default() += c;
}
self.map.retain(|_m, c| !c.is_zero());
}
}
impl Sub for Polynomial {
type Output = Self;
#[inline]
fn sub(mut self, other: Self) -> Self::Output {
self -= other;
self
}
}
impl SubAssign for Polynomial {
fn sub_assign(&mut self, other: Self) {
for (m, c) in other.map {
*self.map.entry(m).or_default() -= c;
}
self.map.retain(|_m, c| !c.is_zero());
}
}
impl Neg for Polynomial {
type Output = Self;
fn neg(mut self) -> Self::Output {
self.map.values_mut().for_each(|c| *c = -*c);
self
}
}
impl Mul for Polynomial {
type Output = Self;
fn mul(self, other: Self) -> Self::Output {
let mut map = BTreeMap::new();
for (lhs_m, lhs_c) in &self.map {
for (rhs_s, rhs_c) in &other.map {
let m = lhs_m.clone() * rhs_s.clone();
let c = *lhs_c * *rhs_c;
*map.entry(m).or_default() += c;
}
}
map.retain(|_m, c: &mut Rational| !c.is_zero());
Self { map }
}
}
impl MulAssign for Polynomial {
fn mul_assign(&mut self, other: Self) {
*self = take(self) * other;
}
}
impl Mul<Rational> for Polynomial {
type Output = Self;
#[inline]
fn mul(mut self, other: Rational) -> Self::Output {
self *= other;
self
}
}
impl MulAssign<Rational> for Polynomial {
fn mul_assign(&mut self, other: Rational) {
if other.is_zero() {
self.map = BTreeMap::default();
} else {
self.map.values_mut().for_each(|c| *c *= other);
}
}
}
impl Mul<i32> for Polynomial {
type Output = Self;
#[inline]
fn mul(mut self, other: i32) -> Self::Output {
self *= other;
self
}
}
impl MulAssign<i32> for Polynomial {
#[inline]
fn mul_assign(&mut self, other: i32) {
*self *= Rational::from(other);
}
}
impl Div<Rational> for Polynomial {
type Output = Self;
#[inline]
fn div(mut self, other: Rational) -> Self::Output {
self /= other;
self
}
}
impl DivAssign<Rational> for Polynomial {
fn div_assign(&mut self, other: Rational) {
assert!(!other.is_zero(), "division by zero");
self.map.values_mut().for_each(|c| *c /= other);
}
}
impl Div<i32> for Polynomial {
type Output = Self;
#[inline]
fn div(mut self, other: i32) -> Self::Output {
self /= other;
self
}
}
impl DivAssign<i32> for Polynomial {
#[inline]
fn div_assign(&mut self, other: i32) {
*self /= Rational::from(other);
}
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub struct Rational {
p: i32,
q: i32,
}
impl Rational {
pub const ZERO: Self = Self { p: 0, q: 1 };
pub const ONE: Self = Self { p: 1, q: 1 };
#[must_use]
pub fn new(mut p: i32, mut q: i32) -> Self {
if p == 0 {
Self::ZERO
} else if p == q {
Self::ONE
} else {
let mut g = p.gcd(q);
assert_ne!(g, 0, "division by zero");
if q < 0 {
g = -g;
}
p /= g;
q /= g;
Self { p, q }
}
}
#[must_use]
#[inline]
pub const fn p(&self) -> i32 {
self.p
}
#[must_use]
#[inline]
pub const fn q(&self) -> i32 {
self.q
}
#[must_use]
pub fn inv(&self) -> Self {
assert_ne!(self.p, 0, "division by zero");
if self.p < 0 {
Self {
p: -self.p,
q: -self.q,
}
} else {
Self {
p: self.q,
q: self.p,
}
}
}
#[must_use]
#[inline]
pub const fn abs(&self) -> Self {
Self {
p: self.p.abs(),
q: self.q,
}
}
#[must_use]
#[inline]
pub const fn is_negative(&self) -> bool {
self.p.is_negative()
}
#[must_use]
#[inline]
pub const fn is_positive(&self) -> bool {
self.p.is_positive()
}
#[must_use]
#[inline]
pub const fn is_one(&self) -> bool {
self.p == 1 && self.q == 1
}
#[must_use]
#[inline]
pub const fn is_zero(&self) -> bool {
self.p == 0
}
#[must_use]
pub fn pow(&self, exp: i32) -> Self {
let abs = exp.unsigned_abs();
let pow = Self::new(self.p().pow(abs), self.q().pow(abs));
if exp.is_negative() { pow.inv() } else { pow }
}
}
impl Default for Rational {
#[inline]
fn default() -> Self {
Self::ZERO
}
}
impl From<i32> for Rational {
#[inline]
fn from(p: i32) -> Self {
Self { p, q: 1 }
}
}
impl From<(i32, i32)> for Rational {
#[inline]
fn from((p, q): (i32, i32)) -> Self {
Self::new(p, q)
}
}
impl Add for Rational {
type Output = Self;
#[inline]
fn add(self, other: Self) -> Self {
Self::new(self.p * other.q + self.q * other.p, self.q * other.q)
}
}
impl AddAssign for Rational {
#[inline]
fn add_assign(&mut self, other: Self) {
*self = *self + other;
}
}
impl Add<i32> for Rational {
type Output = Self;
#[inline]
fn add(self, other: i32) -> Self {
Self::new(self.p + self.q * other, self.q)
}
}
impl AddAssign<i32> for Rational {
#[inline]
fn add_assign(&mut self, other: i32) {
*self = *self + other;
}
}
impl Add<Rational> for i32 {
type Output = Rational;
#[inline]
fn add(self, other: Rational) -> Rational {
Rational::new(self * other.q + other.p, other.q)
}
}
impl Sub for Rational {
type Output = Self;
#[inline]
fn sub(self, other: Self) -> Self {
Self::new(self.p * other.q - self.q * other.p, self.q * other.q)
}
}
impl SubAssign for Rational {
#[inline]
fn sub_assign(&mut self, other: Self) {
*self = *self - other;
}
}
impl Neg for Rational {
type Output = Self;
#[inline]
fn neg(self) -> Self {
Self {
p: -self.p,
q: self.q,
}
}
}
impl Sub<i32> for Rational {
type Output = Self;
#[inline]
fn sub(self, other: i32) -> Self {
Self::new(self.p - self.q * other, self.q)
}
}
impl SubAssign<i32> for Rational {
#[inline]
fn sub_assign(&mut self, other: i32) {
*self = *self - other;
}
}
impl Sub<Rational> for i32 {
type Output = Rational;
#[inline]
fn sub(self, other: Rational) -> Rational {
Rational::new(self * other.q - other.p, other.q)
}
}
impl Mul for Rational {
type Output = Self;
#[inline]
fn mul(self, other: Self) -> Self {
Self::new(self.p * other.p, self.q * other.q)
}
}
impl MulAssign for Rational {
#[inline]
fn mul_assign(&mut self, other: Self) {
*self = *self * other;
}
}
impl Mul<i32> for Rational {
type Output = Self;
#[inline]
fn mul(self, other: i32) -> Self {
Self::new(self.p * other, self.q)
}
}
impl MulAssign<i32> for Rational {
#[inline]
fn mul_assign(&mut self, other: i32) {
*self = *self * other;
}
}
impl Mul<Rational> for i32 {
type Output = Rational;
#[inline]
fn mul(self, other: Rational) -> Rational {
Rational::new(self * other.p, other.q)
}
}
impl Div for Rational {
type Output = Self;
#[inline]
fn div(self, other: Self) -> Self {
Self::new(self.p * other.q, self.q * other.p)
}
}
impl DivAssign for Rational {
#[inline]
fn div_assign(&mut self, other: Self) {
*self = *self / other;
}
}
impl Div<i32> for Rational {
type Output = Self;
#[inline]
#[allow(clippy::suspicious_arithmetic_impl)]
fn div(self, other: i32) -> Self {
Self::new(self.p, self.q * other)
}
}
impl DivAssign<i32> for Rational {
#[inline]
fn div_assign(&mut self, other: i32) {
*self = *self / other;
}
}
impl Div<Rational> for i32 {
type Output = Rational;
#[inline]
#[allow(clippy::suspicious_arithmetic_impl)]
fn div(self, other: Rational) -> Rational {
Rational::new(self * other.q, other.p)
}
}
impl PartialOrd for Rational {
#[inline]
fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
Some(self.cmp(other))
}
}
impl Ord for Rational {
#[inline]
fn cmp(&self, other: &Self) -> Ordering {
(self.p * other.q).cmp(&(other.p * self.q))
}
}
impl Display for Rational {
fn fmt(&self, fmt: &mut fmt::Formatter) -> fmt::Result {
Display::fmt(&self.p, fmt)?;
if self.q != 1 {
write!(fmt, "/{}", self.q)?;
}
Ok(())
}
}
#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord, Default)]
pub struct Monomial {
pub map: BTreeMap<Symbol, NonZeroI32>,
}
impl Monomial {
#[must_use]
#[allow(clippy::missing_panics_doc)]
pub fn one() -> Self {
let mut map = BTreeMap::new();
map.insert(Symbol::one(), NonZero::new(1).unwrap());
Self { map }
}
#[must_use]
#[inline]
pub fn is_zero(&self) -> bool {
self.map.is_empty()
}
#[must_use]
pub fn is_one(&self) -> bool {
self.map
.iter()
.all(|(s, e)| s.is_one() && e.get().abs() == 1)
}
#[must_use]
pub fn alt(self) -> Self {
let map = BTreeMap::new();
let map = self.map.into_iter().fold(map, |mut map, (s, e)| {
map.insert(s.alt(), e);
map
});
Self { map }
}
#[must_use]
pub fn cdm(self, mark: char) -> Self {
let map = BTreeMap::new();
let map = self.map.into_iter().fold(map, |mut map, (s, e)| {
map.insert(s.cdm(mark), e);
map
});
Self { map }
}
#[must_use]
pub fn swp(self) -> Self {
let map = BTreeMap::new();
let map = self.map.into_iter().fold(map, |mut map, (s, e)| {
map.insert(!s, e);
map
});
Self { map }
}
}
impl<S> FromIterator<S> for Monomial
where
S: Into<Symbol>,
{
fn from_iter<M: IntoIterator<Item = S>>(iter: M) -> Self {
Self {
map: iter
.into_iter()
.map(|s| (s.into(), NonZeroI32::new(1).unwrap()))
.collect(),
}
}
}
impl Mul for Monomial {
type Output = Self;
#[inline]
fn mul(mut self, other: Self) -> Self::Output {
self *= other;
self
}
}
impl MulAssign for Monomial {
fn mul_assign(&mut self, other: Self) {
for (s, rhs_e) in other.map {
if let Some(lhs_e) = self.map.get(&s) {
if let Some(e) = NonZeroI32::new(lhs_e.get() + rhs_e.get()) {
assert!(self.map.insert(s, e).is_some());
} else {
assert!(self.map.remove(&s).is_some());
}
} else {
assert!(self.map.insert(s, rhs_e).is_none());
}
}
self.map.retain(|s, _e| !s.is_one());
}
}
impl Div for Monomial {
type Output = Self;
#[inline]
fn div(mut self, other: Self) -> Self::Output {
self /= other;
self
}
}
impl DivAssign for Monomial {
fn div_assign(&mut self, other: Self) {
for (s, rhs_e) in other.map {
if let Some(lhs_e) = self.map.get(&s) {
if let Some(e) = NonZeroI32::new(lhs_e.get() - rhs_e.get()) {
assert!(self.map.insert(s, e).is_some());
} else {
assert!(self.map.remove(&s).is_some());
}
} else {
assert!(self.map.insert(s, rhs_e).is_none());
}
}
self.map.retain(|s, _e| !s.is_one());
}
}
#[derive(Debug, Clone, Copy, Default)]
pub struct Symbol {
var: char,
alt: char,
cdm: char,
lab: &'static str,
}
impl PartialEq for Symbol {
fn eq(&self, other: &Self) -> bool {
self.var == other.var && self.alt == other.alt && self.cdm == other.cdm
}
}
impl Eq for Symbol {}
impl Ord for Symbol {
fn cmp(&self, other: &Self) -> Ordering {
self.var
.cmp(&other.var)
.then(self.alt.cmp(&other.alt))
.then(self.cdm.cmp(&other.cdm))
}
}
impl PartialOrd for Symbol {
fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
Some(self.cmp(other))
}
}
impl Symbol {
pub const NIL: char = '\0';
pub const VEC: char = '\u{200b}';
pub const ALT: char = '\u{0307}';
pub const PIN: char = '\u{0353}';
pub const LHS: char = '\u{0354}';
pub const RHS: char = '\u{0355}';
#[must_use]
#[inline]
pub const fn one() -> Self {
Self {
var: Self::NIL,
alt: Self::NIL,
cdm: Self::NIL,
lab: "",
}
}
#[must_use]
#[inline]
pub const fn is_vec(&self) -> bool {
self.var == Self::VEC
}
#[must_use]
#[inline]
pub const fn is_scalar(&self) -> bool {
self.is_vec() && self.lab.len() == 1
}
#[must_use]
#[inline]
pub const fn is_pseudoscalar(&self) -> bool {
self.is_vec() && self.alt == Self::ALT
}
#[must_use]
#[inline]
pub const fn is_one(&self) -> bool {
self.var == Self::NIL
}
#[must_use]
#[inline]
pub const fn is_pin(&self) -> bool {
self.cdm == Self::PIN
}
#[must_use]
#[inline]
pub const fn is_alt(&self) -> bool {
self.alt == Self::ALT
}
#[must_use]
#[inline]
pub const fn new(var: char, sym: &'static str) -> Self {
Self {
var,
alt: Self::NIL,
cdm: Self::NIL,
lab: sym,
}
}
#[must_use]
#[inline]
pub const fn alt(mut self) -> Self {
self.alt = Self::ALT;
self
}
#[must_use]
#[inline]
pub(crate) const fn cdm(mut self, cdm: char) -> Self {
self.cdm = cdm;
self
}
#[must_use]
#[inline]
pub const fn pin(self) -> Self {
self.cdm(Self::PIN)
}
#[must_use]
#[inline]
pub const fn lhs(self) -> Self {
self.cdm(Self::LHS)
}
#[must_use]
#[inline]
pub const fn rhs(self) -> Self {
self.cdm(Self::RHS)
}
}
impl From<(&str, &'static str)> for Symbol {
#[inline]
fn from((vars, sym): (&str, &'static str)) -> Self {
let mut vars = vars.chars();
let var = vars.next().unwrap_or_default();
assert_eq!(vars.next(), None, "multi-character symbol");
Self::new(var, sym)
}
}
impl From<(char, &'static str)> for Symbol {
#[inline]
fn from((var, sym): (char, &'static str)) -> Self {
Self::new(var, sym)
}
}
impl Not for Symbol {
type Output = Self;
fn not(self) -> Self {
let var = if self.var.is_lowercase() {
let mut iter = self.var.to_uppercase();
assert_eq!(iter.len(), 1, "no uppercase for {}", self.var);
iter.next().unwrap()
} else {
let mut iter = self.var.to_lowercase();
assert_eq!(iter.len(), 1, "no lowercase for {}", self.var);
iter.next().unwrap()
};
Self {
var,
alt: self.alt,
cdm: self.cdm,
lab: self.lab,
}
}
}
impl Display for Symbol {
fn fmt(&self, fmt: &mut fmt::Formatter) -> fmt::Result {
if self.is_vec() {
if fmt.fill() == '$' {
if self.is_pseudoscalar() {
write!(fmt, "\\boldsymbol{{I}}")?;
} else {
let lab = &self.lab[1..];
if lab.is_empty() {
write!(fmt, "\\boldsymbol{{e}}")?;
} else if lab.len() == 1 {
write!(fmt, "\\boldsymbol{{e}}_{lab}")?;
} else {
write!(fmt, "\\boldsymbol{{e}}_{{{lab}}}")?;
}
}
} else if fmt.align() == Some(Alignment::Center) {
write!(fmt, "o.{}", &self.lab)?;
} else if self.is_scalar() {
write!(fmt, "1")?;
} else if self.is_pseudoscalar() {
write!(fmt, "I")?;
} else {
write!(fmt, "{}", &self.lab)?;
}
} else if !self.is_one() {
if fmt.alternate() || fmt.align() == Some(Alignment::Center) || fmt.fill() == '$' {
let var = match self.cdm {
Self::PIN => 'p',
Self::LHS => 'l',
Self::RHS => 'r',
_ => 'v',
};
let lab = &self.lab[1..];
if fmt.fill() == '$' {
if lab.is_empty() {
write!(fmt, "{var}")?;
} else if lab.len() == 1 {
write!(fmt, "{var}_{lab}")?;
} else {
write!(fmt, "{var}_{{{lab}}}")?;
}
} else if fmt.align() == Some(Alignment::Center) {
write!(fmt, "{var}.e{lab}")?;
} else {
write!(fmt, "{var}{lab}")?;
}
} else {
write!(fmt, "{}", self.var)?;
if self.alt != Self::NIL {
write!(fmt, "{}", self.alt)?;
}
if self.cdm != Self::NIL {
write!(fmt, "{}", self.cdm)?;
}
}
}
Ok(())
}
}
#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord)]
pub enum Tree {
Add(Vec<Self>),
Mul(Vec<Self>),
Num(Rational),
Sym(Symbol),
}
impl Tree {
pub const ZERO: Self = Self::Num(Rational::ZERO);
pub const ONE: Self = Self::Num(Rational::ONE);
#[must_use]
#[allow(clippy::missing_panics_doc)]
pub fn with_factorization<B: Algebra>(v: Multivector<B>, signed: bool) -> Self {
let mut add = Vec::with_capacity(v.map.len());
for (b, p) in v.map {
let f = Factorization::new(p, signed);
if f.is_one() {
add.push(b.into().into());
} else {
add.push(Self::Mul(vec![f.into(), b.into().into()]));
}
}
if add.is_empty() {
Self::ZERO
} else if add.len() == 1 {
add.pop().expect("unreachable")
} else {
Self::Add(add)
}
}
#[must_use]
#[inline]
pub const fn as_num(&self) -> Option<Rational> {
if let Self::Num(num) = self {
Some(*num)
} else {
None
}
}
#[must_use]
#[inline]
pub const fn as_sym(&self) -> Option<Symbol> {
if let Self::Sym(sym) = self {
Some(*sym)
} else {
None
}
}
}
impl Default for Tree {
#[inline]
fn default() -> Self {
Self::ZERO
}
}
impl<B: Algebra> From<Multivector<B>> for Tree {
fn from(v: Multivector<B>) -> Self {
let add = v
.map
.into_iter()
.filter(|(_b, p)| !p.is_zero())
.map(|(b, p)| {
if p.is_one() {
b.into().into()
} else {
Self::Mul(vec![p.into(), b.into().into()])
}
})
.collect::<Vec<Self>>();
if add.is_empty() {
Self::ZERO
} else if add.len() == 1 {
add.into_iter().next().expect("unreachable")
} else {
Self::Add(add)
}
}
}
impl From<Factorization> for Tree {
fn from(f: Factorization) -> Self {
let add = f
.map
.into_iter()
.filter(|(m, (p, c))| !(m.is_zero() || p.is_zero() || c.is_zero()))
.map(|(m, (p, c))| {
let is_mul = [c.is_one(), p.is_one(), m.is_one()].map(bool::not);
let len = is_mul.into_iter().map(usize::from).sum::<usize>();
let mut mul = Vec::with_capacity(len);
if is_mul[0] {
mul.push(c.into());
}
if is_mul[1] {
mul.push(p.into());
}
if is_mul[2] {
mul.push(m.into());
}
if len == 1 {
mul.pop().expect("unreachable")
} else {
Self::Mul(mul)
}
})
.collect::<Vec<Self>>();
if add.is_empty() {
Self::ZERO
} else if add.len() == 1 {
let any = add.into_iter().next().expect("unreachable");
if f.gcd.is_one() {
any
} else if let Self::Num(num) = any {
Self::Num(f.gcd * num)
} else {
Self::Mul(vec![f.gcd.into(), any])
}
} else if f.gcd.is_one() {
Self::Add(add)
} else {
Self::Mul(vec![f.gcd.into(), Self::Add(add)])
}
}
}
impl From<Polynomial> for Tree {
fn from(p: Polynomial) -> Self {
let add = p
.map
.into_iter()
.filter(|(p, c)| !(p.is_zero() || c.is_zero()))
.map(|(m, c)| {
if c.is_one() {
m.into()
} else {
let m = Self::from(m);
if let Self::Mul(mul) = m {
let mut vec = vec![Self::from(c)];
vec.extend(mul);
Self::Mul(vec)
} else if let Self::Num(num) = m {
Self::Num(c * num)
} else {
Self::Mul(vec![c.into(), m])
}
}
})
.collect::<Vec<Self>>();
if add.is_empty() {
Self::ZERO
} else if add.len() == 1 {
add.into_iter().next().expect("unreachable")
} else {
Self::Add(add)
}
}
}
impl From<Monomial> for Tree {
fn from(m: Monomial) -> Self {
let mul = m
.map
.into_iter()
.filter(|(s, _e)| !s.is_one())
.flat_map(|(s, e)| {
repeat_n(s, e.get().try_into().expect("negative exponent")).map(From::from)
})
.collect::<Vec<Self>>();
if mul.is_empty() {
Self::ONE
} else if mul.len() == 1 {
mul.into_iter().next().expect("unreachable")
} else {
Self::Mul(mul)
}
}
}
impl From<Rational> for Tree {
#[inline]
fn from(r: Rational) -> Self {
Self::Num(r)
}
}
impl From<Symbol> for Tree {
#[inline]
fn from(s: Symbol) -> Self {
if s.is_one() { Self::ONE } else { Self::Sym(s) }
}
}
pub mod pga;
pub type PgaE0 = Multivector<pga::PgaE0>;
pub type PgaE1 = Multivector<pga::PgaE1>;
pub type PgaE2 = Multivector<pga::PgaE2>;
pub type PgaE3 = Multivector<pga::PgaE3>;
pub type PgaE4 = Multivector<pga::PgaE4>;
pub type PgaE5 = Multivector<pga::PgaE5>;
pub type PgaE6 = Multivector<pga::PgaE6>;
pub type PgaE7 = Multivector<pga::PgaE7>;
pub type PgaH0 = Multivector<pga::PgaH0>;
pub type PgaH1 = Multivector<pga::PgaH1>;
pub type PgaH2 = Multivector<pga::PgaH2>;
pub type PgaH3 = Multivector<pga::PgaH3>;
pub type PgaH4 = Multivector<pga::PgaH4>;
pub type PgaH5 = Multivector<pga::PgaH5>;
pub type PgaH6 = Multivector<pga::PgaH6>;
pub type PgaH7 = Multivector<pga::PgaH7>;
pub type PgaP0 = Multivector<pga::PgaP0>;
pub type PgaP1 = Multivector<pga::PgaP1>;
pub type PgaP2 = Multivector<pga::PgaP2>;
pub type PgaP3 = Multivector<pga::PgaP3>;
pub type PgaP4 = Multivector<pga::PgaP4>;
pub type PgaP5 = Multivector<pga::PgaP5>;
pub type PgaP6 = Multivector<pga::PgaP6>;
pub type PgaP7 = Multivector<pga::PgaP7>;