uiua_parser 0.19.0-dev.4

use std::{
    cell::RefCell,
    cmp::Ordering,
    collections::HashMap,
    fmt,
    iter::{self, repeat_n},
    mem::{replace, take},
    ops::*,
    slice,
};

use ecow::{EcoVec, eco_vec};
use serde::*;

use crate::{Complex, SUBSCRIPT_DIGITS, ga::*};

#[derive(Debug, Clone, Default, Serialize, Deserialize)]
pub struct Multivector {
    #[serde(default, skip_serializing_if = "EcoVec::is_empty", rename = "c")]
    coefs: EcoVec<f64>,
    /// The geometric algebra flavor
    #[serde(default, skip_serializing_if = "is_vanilla", rename = "f")]
    pub flavor: Flavor,
}

fn is_vanilla(flavor: &Flavor) -> bool {
    matches!(flavor, Flavor::Vanilla)
}

impl Multivector {
    /// Set the maximum number of dimensions
    pub fn set_dims(&mut self, new_dims: u8) {
        let dims = self.dims();
        match dims.cmp(&new_dims) {
            Ordering::Equal => {}
            Ordering::Less => {
                // Expansion
                self.coefs
                    .extend(repeat_n(0.0, (1 << new_dims) - (1 << dims)));
                let slice = self.coefs.make_mut();
                let mut left = 0;
                for d in 0..=dims {
                    let ai = grade_size(dims, d);
                    let bi = grade_size(new_dims, d);
                    slice[left + ai..].rotate_right(bi - ai);
                    left += bi;
                }
                if dims == 3 {
                    self.set_blade(0b101, -self.get_blade(0b101))
                }
            }
            Ordering::Greater => {
                // Contraction
                let slice = self.coefs.make_mut();
                let dim_mask = (1 << new_dims) - 1;
                let mut left = 0;
                let (mask_table, _) = mask_tables(dims);
                for i in 0..(1 << dims) {
                    if mask_table[i] & !dim_mask == 0 {
                        slice[left] = slice[i];
                        left += 1;
                    }
                }
                self.coefs.truncate(1 << new_dims);
                if new_dims == 3 {
                    self.set_blade(0b101, -self.get_blade(0b101))
                }
            }
        }
    }
    fn set_dims_right(&mut self, dims: u8) {
        let low_dims = self.dims();
        if low_dims >= dims {
            return;
        }
        self.coefs
            .extend(repeat_n(0.0, (1 << dims) - (1 << low_dims)));
        let slice = self.coefs.make_mut();
        let mut b_left = 0;
        for d in 0..=low_dims {
            let ai = grade_size(low_dims, d);
            let bi = grade_size(dims, d);
            slice[b_left..].rotate_right(bi - ai);
            b_left += bi;
        }
    }
    pub fn set_flavor(&mut self, flavor: Flavor) {
        match (self.flavor, flavor) {
            (Flavor::Vanilla, Flavor::Projective) => self.set_dims_right(self.dims() + 1),
            _ => {}
        }
        self.flavor = flavor;
    }
    fn conform_to(&mut self, other: &Self) {
        match (self.flavor, other.flavor) {
            (Flavor::Vanilla, Flavor::Projective) => self.set_flavor(Flavor::Projective),
            _ => {}
        }
        let (a_dims, b_dims) = (self.dims(), other.dims());
        if a_dims < b_dims {
            self.set_dims(b_dims);
        }
    }
    fn conform(&mut self, other: &mut Self) {
        match (self.flavor, other.flavor) {
            (Flavor::Vanilla, Flavor::Projective) => self.set_flavor(Flavor::Projective),
            (Flavor::Projective, Flavor::Vanilla) => other.set_flavor(Flavor::Projective),
            _ => {}
        }
        let (a_dims, b_dims) = (self.dims(), other.dims());
        match a_dims.cmp(&other.dims()) {
            Ordering::Equal => {}
            Ordering::Greater => other.set_dims(a_dims),
            Ordering::Less => self.set_dims(b_dims),
        }
    }
    /// Get the maximum number of dimensions
    pub fn dims(&self) -> u8 {
        (self.coefs.len().max(1) as f64).log2() as u8
    }
    /// Create a multivector from the real and imaginary parts of a complex number
    pub fn complex(re: impl Into<f64>, im: impl Into<f64>) -> Self {
        Self::all([re.into(), 0.0, 0.0, im.into()])
    }
    /// Create a multivector from all of its coefficients
    pub fn all(coefs: impl Into<EcoVec<f64>>) -> Self {
        let coefs = coefs.into();
        assert!(
            coefs.len().max(1).is_power_of_two(),
            "Must have power of 2 number of coefficients, but there are {}",
            coefs.len()
        );
        Multivector {
            coefs,
            flavor: Flavor::Vanilla,
        }
    }
    /// Create a multivector from its vector coefficients
    pub fn vector(vector: impl Into<EcoVec<f64>>, flavor: Flavor) -> Self {
        let mut coefs = vector.into();
        let dims = coefs.len();
        coefs.insert(0, 0.0);
        coefs.extend(repeat_n(0.0, (1 << dims) - 1 - dims));
        Multivector { coefs, flavor }
    }
    pub fn vga_vector(vector: impl Into<EcoVec<f64>>) -> Self {
        Self::vector(vector.into(), Flavor::Vanilla)
    }
    pub fn pga_vector(vector: impl Into<EcoVec<f64>>) -> Self {
        Self::vector(vector.into(), Flavor::Projective)
    }
    /// Create a multivector from its n-1 blade coefficients
    pub fn n_1_blades(blades: impl Into<EcoVec<f64>>, flavor: Flavor) -> Self {
        let mut coefs = blades.into();
        let dims = coefs.len();
        coefs.insert(0, 0.0);
        coefs.extend(repeat_n(0.0, (1 << dims) - 1 - dims));
        coefs.make_mut().rotate_left(dims + 2);
        Multivector { coefs, flavor }
    }
    pub fn vga_n_1_blades(blades: impl Into<EcoVec<f64>>) -> Self {
        Self::n_1_blades(blades.into(), Flavor::Vanilla)
    }
    pub fn pga_n_1_blades(blades: impl Into<EcoVec<f64>>) -> Self {
        Self::n_1_blades(blades.into(), Flavor::Projective)
    }
    /// Create a multivector from some low-grade, even, or all blades
    pub fn blades_left(
        dims: u8,
        blades: impl Into<EcoVec<f64>>,
        flavor: Flavor,
    ) -> Result<Self, EcoVec<f64>> {
        Self::blades_impl(dims, blades.into(), false, flavor)
    }
    /// Create a multivector from some high-grade, even, or all blades
    pub fn blades_right(
        dims: u8,
        blades: impl Into<EcoVec<f64>>,
        flavor: Flavor,
    ) -> Result<Self, EcoVec<f64>> {
        Self::blades_impl(dims, blades.into(), true, flavor)
    }
    fn blades_impl(
        dims: u8,
        mut blades: EcoVec<f64>,
        odd: bool,
        flavor: Flavor,
    ) -> Result<Self, EcoVec<f64>> {
        let blade_count = 1 << dims;
        Ok(if blades.len() == blade_count {
            Self::all(blades)
        } else if blades.len() * 2 == blade_count {
            blades.extend(repeat_n(0.0, blade_count / 2));
            let slice = blades.make_mut();
            let mut left = 0;
            for g in 0..=dims {
                let i = grade_size(dims, g);
                if g % 2 != odd as u8 {
                    slice[left..].rotate_right(i)
                }
                left += i;
            }
            Self::all(blades).flavor(flavor)
        } else {
            let (start, end) = if odd {
                (((dims + 1) as f32 / 2.0).floor() as u8, dims + 1)
            } else {
                (0, ((dims + 1) as f32 / 2.0).ceil() as u8)
            };
            let mut left: usize = (0..start).map(|d| grade_size(dims, d)).sum();
            for g in start..end {
                let grade_size = grade_size(dims, g);
                if grade_size == blades.len() {
                    blades.extend(repeat_n(0.0, blade_count - grade_size));
                    blades.make_mut().rotate_right(left);
                    return Ok(Self::all(blades).flavor(flavor));
                }
                left += grade_size;
            }
            return Err(blades);
        })
    }
    pub fn scalar(mut dims: u8, f: f64, flavor: Flavor) -> Self {
        dims += flavor.dim_adjustment();
        let mut coefs = eco_vec![0.0; 1 << dims];
        coefs.make_mut()[0] = f;
        Multivector { coefs, flavor }
    }
    /// Create a unit pseudoscalar
    pub fn pseudo_unit(dims: u8, flavor: Flavor) -> Self {
        Self::pseudoscalar(dims, 1.0, flavor)
    }
    pub fn vga_pseudo_unit(dims: u8) -> Self {
        Self::pseudo_unit(dims, Flavor::Vanilla)
    }
    pub fn pga_pseudo_unit(dims: u8) -> Self {
        Self::pseudo_unit(dims, Flavor::Projective)
    }
    pub fn vga_pseudoscalar(dims: u8, n: f64) -> Self {
        Self::pseudoscalar(dims, n, Flavor::Vanilla)
    }
    pub fn pga_pseudoscalar(dims: u8, n: f64) -> Self {
        Self::pseudoscalar(dims, n, Flavor::Projective)
    }
    /// Create a pseudoscalar
    pub fn pseudoscalar(mut dims: u8, n: f64, flavor: Flavor) -> Self {
        dims += flavor.dim_adjustment();
        let mut coefs = eco_vec![0.0; 1 << dims];
        *coefs.make_mut().last_mut().unwrap() = n;
        Multivector { coefs, flavor }
    }
    /// Set the flavor
    pub fn flavor(mut self, flavor: Flavor) -> Self {
        self.flavor = flavor;
        self
    }
    /// Make this multivector use Projective Geometric Algebra
    pub fn pga(self) -> Self {
        self.flavor(Flavor::Projective)
    }
    /// Check if the multivector is a scalar
    pub fn is_scalar(&self) -> bool {
        !self.coefs.is_empty() && self.coefs.iter().skip(1).all(|&c| c == 0.0)
    }
    /// Try to get the multivector as a scalar
    pub fn as_scalar(&self) -> Option<f64> {
        self.is_scalar().then(|| self[0])
    }
    /// Try to get the multivector as a complex number
    pub fn as_complex(&self) -> Option<Complex> {
        (self.flavor == Flavor::Vanilla
            && self.coefs.len() == 4
            && self[1] == 0.0
            && self[2] == 0.0)
            .then(|| Complex::new(self[0], self[3]))
    }
    pub fn iter(&self) -> iter::Copied<slice::Iter<'_, f64>> {
        self.coefs.iter().copied()
    }
    pub fn iter_mut(&mut self) -> slice::IterMut<'_, f64> {
        self.coefs.make_mut().iter_mut()
    }
    /// Geometrically dual the multivector
    pub fn dual(&mut self) {
        let dims = self.dims();
        self.product_impl(
            Self::vga_pseudo_unit(dims),
            Flavor::Vanilla,
            MetricMode::All,
        );
    }
    /// Get the geometric dual of the multivector
    pub fn dualed(mut self) -> Self {
        self.dual();
        self
    }
    /// Geometrically antidual the multivector
    pub fn antidual(&mut self) {
        let flavor = take(&mut self.flavor);
        let dims = self.dims();
        let v = replace(self, Self::vga_pseudo_unit(dims));
        *self *= v;
        self.flavor = flavor;
    }
    /// Get the geometric antidual of the multivector
    pub fn antidualed(mut self) -> Self {
        self.antidual();
        self
    }
    /// Geometrically reverse the multivector
    pub fn reverse(&mut self) {
        for (g, c) in blade_grades(self.dims()).zip(self) {
            if g / 2 % 2 == 1 {
                *c = -*c;
            }
        }
    }
    /// Geometrically inverse reverse the multivector
    pub fn inv_reverse(&mut self) {
        for (g, c) in blade_grades(self.dims()).zip(self) {
            if g / 2 % 2 == 0 {
                *c = -*c;
            }
        }
    }
    /// Get the geometric reverse of the multivector
    pub fn reversed(mut self) -> Self {
        self.reverse();
        self
    }
    /// Get the geometric inverse reverse of the multivector
    pub fn inv_reversed(mut self) -> Self {
        self.inv_reverse();
        self
    }
    /// Get the squared magnitude of the multivector
    pub fn squared_magnitude(&self) -> f64 {
        let dims = self.dims();
        let (mask_table, _) = mask_tables(dims);
        (self.iter().enumerate())
            .map(|(i, c)| blade_metric_scalar_only(mask_table[i], dims, self.flavor) * c * c)
            .sum::<f64>()
    }
    /// Get the magnitude of the multivector
    pub fn magnitude(&self) -> f64 {
        self.squared_magnitude().sqrt()
    }
    /// Normalize the multivector
    pub fn normalize(&mut self) {
        let mag = self.clone().magnitude();
        *self /= mag;
        if let Flavor::Projective = self.flavor {
            for mask in [0b1, 0b110] {
                let blade = self.get_blade(mask);
                if blade != 0.0 && blade != 1.0 {
                    for c in self.coefs.make_mut() {
                        *c /= blade;
                    }
                    break;
                }
            }
        }
    }
    /// Get the normalized multivector
    pub fn normalized(mut self) -> Self {
        self.normalize();
        self
    }
    /// Get the rotor from this multivector to another
    pub fn rotor(self, to: Self) -> Self {
        ((self * to.reversed()).normalized() + 1.0).normalized()
    }
    pub fn inner_product(mut self, mut other: Self) -> Self {
        self.conform(&mut other);
        self.product_impl(other, self.flavor, MetricMode::Inner);
        self
    }
    pub fn left_contraction(mut self, mut other: Self) -> Self {
        self.conform(&mut other);
        self.product_impl(other, self.flavor, MetricMode::Left);
        self
    }
    pub fn right_contraction(mut self, mut other: Self) -> Self {
        self.conform(&mut other);
        self.product_impl(other, self.flavor, MetricMode::Right);
        self
    }
    pub fn outer_product(mut self, mut other: Self) -> Self {
        self.conform(&mut other);
        self.product_impl(other, Flavor::NULL, MetricMode::All);
        self
    }
    pub fn regressive_product(self, other: Self) -> Self {
        self.dualed().outer_product(other.dualed()).antidualed()
    }
    fn product(&mut self, mut rhs: Self) {
        self.conform(&mut rhs);
        self.product_impl(rhs, self.flavor, MetricMode::All)
    }
    fn square(&mut self) {
        self.product_impl(self.clone(), self.flavor, MetricMode::All)
    }
    fn squared(mut self) -> Self {
        self.square();
        self
    }
    fn product_impl(&mut self, rhs: Self, flavor: Flavor, mode: MetricMode) {
        let (a, b) = (self, rhs);
        let dims = a.dims();
        let mut new_coefs = eco_vec![0.0; a.coefs.len()];
        let slice = new_coefs.make_mut();
        let (mask_table, inv_mask_table) = mask_tables(dims);
        let metrics = blade_metrics(dims, flavor, mode);
        for (bi, &b_mask) in mask_table.iter().enumerate() {
            for (ai, &a_mask) in mask_table.iter().enumerate() {
                let metric = metrics[a_mask * mask_table.len() + b_mask];
                if metric == 0.0 {
                    continue;
                }
                let c_mask = a_mask ^ b_mask;
                let ci = inv_mask_table[c_mask];
                slice[ci] += metric * a[ai] * b[bi];
                // println!(
                //     "metric:{}{metric}, a: {}, b: {}",
                //     if a[ai] != 0.0 && b[bi] != 0.0 {
                //         "*"
                //     } else {
                //         " "
                //     },
                //     a[ai],
                //     b[bi]
                // );
            }
        }
        for c in slice {
            let fract = c.abs().fract();
            if fract < 1e-14 || 1.0 - fract < 1e-14 {
                *c = c.round();
            }
        }
        a.coefs = new_coefs;
    }
    /// Get an arbitrary blade given a mask of which bases are present in it
    pub fn get_blade(&self, mask: usize) -> f64 {
        let dims = self.dims();
        let dim_mask = (1 << dims) - 1;
        if (mask & !dim_mask).count_ones() > 0 {
            0.0
        } else {
            let (_, inv_mask_table) = mask_tables(dims);
            self[inv_mask_table[mask]]
        }
    }
    pub fn set_blade(&mut self, mask: usize, coef: f64) {
        let dims = self.dims();
        let dim_mask = (1 << dims) - 1;
        if (mask & !dim_mask).count_ones() > 0 {
            self.set_dims(dims + 1);
            self.set_blade(mask, coef);
        } else {
            let (_, inv_mask_table) = mask_tables(dims);
            self[inv_mask_table[mask]] = coef;
        }
    }
    pub fn grades(&self) -> impl Iterator<Item = &[f64]> {
        let dims = self.dims();
        let mut start = 0;
        (0..=dims).map(move |g| {
            let size = grade_size(dims, g);
            let slice = &self.coefs[start..][..size];
            start += size;
            slice
        })
    }
    /// Get the scalar and set it to 0
    pub fn take_scalar(&mut self) -> f64 {
        take(&mut self[0])
    }
    #[doc(hidden)]
    pub fn make_scalar_only(&mut self) {
        for c in self.coefs.make_mut().iter_mut().skip(1) {
            *c = 0.0;
        }
    }
    fn nanify(&mut self) {
        self.coefs.make_mut().fill(0.0);
        self.set_blade(0b0, f64::NAN);
    }
    pub fn inverted(mut self) -> Self {
        if let Some(f) = self.as_scalar() {
            self.set_blade(0, 1.0 / f);
            self
        } else if let Some(mut c) = self.as_complex() {
            c = c.recip();
            self.set_blade(0b0, c.re);
            self.set_blade(0b11, c.im);
            self
        } else {
            let s = self.get_blade(0b0);
            let mut x = self - s;
            if let Some(sqr) = x.clone().squared().as_scalar() {
                (-x + s) / (s * s - sqr)
            } else {
                x.nanify();
                x
            }
        }
    }
    pub fn exp(mut self) -> Self {
        if let Some(x) = self.as_scalar() {
            // Scalar
            self[0] = x.exp();
            self
        } else if let Some(c) = self.as_complex() {
            // Complex
            let exp_re = c.re.exp();
            self[0] = exp_re * c.im.cos();
            self[3] = exp_re * c.im.sin();
            self
        } else {
            let s = self.get_blade(0b0);
            let x = self - s;
            if let Some(sqr) = x.clone().squared().as_scalar() {
                if sqr < 0.0 {
                    // Scalar + Euclidan bivector
                    let λ = (-sqr).sqrt();
                    (x / λ * λ.sin() + λ.cos()) * s.exp()
                } else if sqr > 0.0 {
                    // Scalar + non-Euclidean bivector
                    let λ = sqr.sqrt();
                    (x / λ * λ.sinh() + λ.cosh()) * s.exp()
                } else {
                    // Scalar + nilpotent
                    (x + 1.0) * s.exp()
                }
            } else {
                // Power series
                self = x + s;
                let mut term = Multivector::scalar(self.dims(), 1.0, self.flavor);
                let mut sum = term.clone();
                for n in 1..100 {
                    term *= self.clone() / n as f64;
                    if term.squared_magnitude() <= f64::EPSILON {
                        break;
                    }
                    sum += term.clone();
                }
                sum
            }
        }
    }
    pub fn ln(mut self) -> Self {
        if let Some(x) = self.as_scalar() {
            // Scalar
            self[0] = x.ln();
            self
        } else if let Some(c) = self.as_complex() {
            // Complex
            let c_ln = c.ln();
            self[0] = c_ln.re;
            self[3] = c_ln.im;
            self
        } else {
            let s = self.get_blade(0b0);
            let mut x = self - s;
            if let Some(sqr) = x.clone().squared().as_scalar() {
                if sqr < 0.0 {
                    // Circular
                    let λ = (-sqr).sqrt();
                    x / λ * λ.atan2(s) + (s * s - sqr).sqrt().ln()
                } else if sqr > 0.0 {
                    // Hyperbolic
                    let λ = sqr.sqrt();
                    x / λ * (λ / s).atanh() + (s * s - sqr).sqrt().ln()
                } else {
                    // Nilpotent
                    x / s + s.ln()
                }
            } else {
                x.nanify();
                x
            }
        }
    }
    pub fn sqrt(self) -> Self {
        self.powf(0.5)
    }
    pub fn powi(mut self, mut power: i32) -> Self {
        match power {
            0 => return Self::scalar(self.dims(), 1.0, self.flavor),
            1 => return self,
            2 => return self.squared(),
            -1 => return self.inverted(),
            -2 => return self.inverted().squared(),
            _ => {}
        }
        if power < 0 {
            self = self.inverted();
            power = -power;
        }
        let mut res = Self::scalar(self.dims(), 1.0, self.flavor);
        while power > 0 {
            if power % 2 == 1 {
                res *= self.clone();
            }
            power /= 2;
            if power > 0 {
                self = self.clone().squared();
            }
        }
        res
    }
    pub fn powf(self, power: f64) -> Self {
        if power.fract() == 0.0 {
            return self.powi(power as i32);
        }
        (self.ln() * power).exp()
    }
    pub fn powmv(mut self, power: Self) -> Self {
        if let Some(p) = power.as_scalar() {
            self.powf(p)
        } else {
            self.nanify();
            self
        }
    }
}

fn blade_metric_scalar_only(a_mask: usize, dims: u8, flavor: Flavor) -> f64 {
    let mut metric = 1.0;
    // Calc sign based on swaps
    let grade = a_mask.count_ones();
    if grade >= 2 {
        let swaps = grade * (grade - 1) / 2;
        if swaps % 2 == 1 {
            metric = -metric;
        }
    }
    // Reverse if necessary
    if grade / 2 % 2 == 1 {
        metric = -metric;
    }
    // Apply flavor metrics
    for i in 0..dims {
        if a_mask & (1 << i) != 0 {
            metric *= flavor.metric(i) as f64;
        }
    }
    metric
}

#[derive(Clone, Copy, PartialEq, Eq, Hash)]
enum MetricMode {
    All,
    Left,
    Right,
    Inner,
}

/// Calculate, cache, and access the blade metrics for a given number of dimensions and GA flavor
fn blade_metrics(dims: u8, flavor: Flavor, mode: MetricMode) -> &'static [f64] {
    type Cache = HashMap<(u8, Flavor, MetricMode), &'static [f64]>;
    thread_local! {
        static CACHE: RefCell<Cache> = Default::default();
    }
    let blade_count = 1 << dims;
    CACHE.with(move |r| {
        *(r.borrow_mut().entry((dims, flavor, mode))).or_insert_with(|| {
            let mut metrics = Vec::with_capacity(blade_count * blade_count);
            for a in 0..blade_count {
                for b in 0..blade_count {
                    let has_metric = match mode {
                        MetricMode::All => true,
                        MetricMode::Left => a == 0 || a & b == a,
                        MetricMode::Right => b == 0 || a & b == b,
                        MetricMode::Inner => a == 0 || b == 0 || a & b == a || a & b == b,
                    };
                    if !has_metric {
                        metrics.push(0.0);
                        continue;
                    }
                    let mut metric = 1.0;
                    if dims == 3 {
                        // Account for reording of e02->e20 to
                        // keep a nice cyclical order of e12,e23,e31.
                        // Such an order does not exist for > 3 dims.
                        let ab = a ^ b;
                        for i in [a, b, ab] {
                            if (i ^ (i >> 1)).count_ones() == dims as u32 {
                                metric = -metric;
                            }
                        }
                    }
                    for i in 0..dims {
                        let bit_i = 1 << i;
                        // Calc sign based on swaps
                        if a & bit_i != 0 {
                            let lower_bits = b & (bit_i - 1);
                            if lower_bits.count_ones() % 2 == 1 {
                                metric = -metric;
                            }
                        }
                        // Apply flavor metrics
                        if (a & bit_i != 0) && (b & bit_i != 0) {
                            metric *= flavor.metric(i) as f64;
                        }
                    }
                    metrics.push(metric);
                }
            }
            metrics.leak()
        })
    })
}

fn eq(a: f64, b: f64) -> bool {
    a == b || a.is_nan() && b.is_nan()
}

impl PartialEq for Multivector {
    fn eq(&self, other: &Self) -> bool {
        self.as_scalar() == other.as_scalar()
            || self.coefs.len() == other.coefs.len()
                && self.coefs.iter().zip(&other.coefs).all(|(a, b)| eq(*a, *b))
    }
}

impl PartialEq<Complex> for Multivector {
    fn eq(&self, other: &Complex) -> bool {
        self.as_complex().is_some_and(|a| a == *other)
    }
}

impl PartialEq<Multivector> for Complex {
    fn eq(&self, other: &Multivector) -> bool {
        other.as_complex().is_some_and(|b| *self == b)
    }
}

impl Eq for Multivector {}

impl PartialOrd for Multivector {
    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
        Some(self.cmp(other))
    }
}

impl Ord for Multivector {
    fn cmp(&self, other: &Self) -> Ordering {
        if let Some(a) = self.as_scalar()
            && let Some(b) = other.as_scalar()
        {
            return (a.partial_cmp(&b)).unwrap_or_else(|| a.is_nan().cmp(&b.is_nan()));
        }
        self.coefs.len().cmp(&other.coefs.len()).then_with(|| {
            self.iter()
                .zip(other)
                .map(|(a, b)| {
                    a.partial_cmp(&b)
                        .unwrap_or_else(|| a.is_nan().cmp(&b.is_nan()))
                })
                .find(|&o| o != Ordering::Equal)
                .unwrap_or(Ordering::Equal)
        })
    }
}

impl PartialEq<[f64]> for Multivector {
    fn eq(&self, other: &[f64]) -> bool {
        self.coefs == other
    }
}
impl<const N: usize> PartialEq<[f64; N]> for Multivector {
    fn eq(&self, other: &[f64; N]) -> bool {
        self.coefs == other
    }
}

impl Neg for Multivector {
    type Output = Self;
    fn neg(mut self) -> Self::Output {
        for c in &mut self {
            *c = -*c;
        }
        self
    }
}

impl AddAssign<f64> for Multivector {
    fn add_assign(&mut self, rhs: f64) {
        if self.coefs.is_empty() {
            self.coefs = eco_vec![rhs]
        } else {
            self[0] += rhs
        }
    }
}
impl Add<f64> for Multivector {
    type Output = Self;
    fn add(mut self, rhs: f64) -> Self::Output {
        self += rhs;
        self
    }
}

impl AddAssign for Multivector {
    fn add_assign(&mut self, mut rhs: Self) {
        if let Some(b) = rhs.as_scalar() {
            self.conform_to(&rhs);
            *self += b;
        } else {
            self.conform(&mut rhs);
            for (a, b) in self.coefs.make_mut().iter_mut().zip(rhs.coefs) {
                *a += b;
            }
        }
    }
}
impl Add for Multivector {
    type Output = Self;
    fn add(mut self, rhs: Self) -> Self::Output {
        self += rhs;
        self
    }
}

impl SubAssign for Multivector {
    fn sub_assign(&mut self, mut rhs: Self) {
        if !self.coefs.is_empty()
            && let Some(b) = rhs.as_scalar()
        {
            self.conform_to(&rhs);
            self[0] -= b;
        } else {
            self.conform(&mut rhs);
            for (a, b) in self.coefs.make_mut().iter_mut().zip(rhs.coefs) {
                *a -= b;
            }
        }
    }
}
impl SubAssign<f64> for Multivector {
    fn sub_assign(&mut self, rhs: f64) {
        *self += -rhs;
    }
}
impl Sub for Multivector {
    type Output = Self;
    fn sub(mut self, rhs: Self) -> Self::Output {
        self -= rhs;
        self
    }
}
impl Sub<f64> for Multivector {
    type Output = Self;
    fn sub(mut self, rhs: f64) -> Self::Output {
        self -= rhs;
        self
    }
}

impl MulAssign for Multivector {
    fn mul_assign(&mut self, rhs: Self) {
        if let Some(b) = rhs.as_scalar() {
            self.conform_to(&rhs);
            *self *= b;
        } else {
            self.product(rhs)
        }
    }
}
impl MulAssign<f64> for Multivector {
    fn mul_assign(&mut self, rhs: f64) {
        if rhs == 1.0 {
            return;
        }
        for c in self {
            *c *= rhs;
        }
    }
}
impl Mul for Multivector {
    type Output = Self;
    fn mul(mut self, rhs: Self) -> Self::Output {
        self *= rhs;
        self
    }
}
impl Mul<f64> for Multivector {
    type Output = Self;
    fn mul(mut self, rhs: f64) -> Self::Output {
        self *= rhs;
        self
    }
}

impl DivAssign for Multivector {
    fn div_assign(&mut self, rhs: Self) {
        if let Some(f) = rhs.as_scalar() {
            *self /= f;
        } else if let Some(c) = rhs.as_complex() {
            *self /= c;
        } else {
            *self *= rhs.inverted()
        }
    }
}
impl DivAssign<Complex> for Multivector {
    #[allow(clippy::suspicious_op_assign_impl)]
    fn div_assign(&mut self, rhs: Complex) {
        *self *= Self::from(rhs.recip());
    }
}
impl DivAssign<f64> for Multivector {
    fn div_assign(&mut self, rhs: f64) {
        *self *= 1.0 / rhs;
    }
}
impl Div<Complex> for Multivector {
    type Output = Self;
    fn div(mut self, rhs: Complex) -> Self::Output {
        self /= rhs;
        self
    }
}
impl Div<f64> for Multivector {
    type Output = Self;
    fn div(mut self, rhs: f64) -> Self::Output {
        self /= rhs;
        self
    }
}
impl Div for Multivector {
    type Output = Self;
    fn div(mut self, rhs: Self) -> Self::Output {
        self /= rhs;
        self
    }
}

#[derive(Debug, Clone)]
pub struct UnableToInvert<M = Multivector>(pub M);

impl Index<usize> for Multivector {
    type Output = f64;
    fn index(&self, index: usize) -> &Self::Output {
        &self.coefs[index]
    }
}
impl IndexMut<usize> for Multivector {
    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
        &mut self.coefs.make_mut()[index]
    }
}
impl<'a> IntoIterator for &'a Multivector {
    type Item = f64;
    type IntoIter = iter::Copied<slice::Iter<'a, f64>>;
    fn into_iter(self) -> Self::IntoIter {
        self.iter()
    }
}
impl<'a> IntoIterator for &'a mut Multivector {
    type Item = &'a mut f64;
    type IntoIter = slice::IterMut<'a, f64>;
    fn into_iter(self) -> Self::IntoIter {
        self.iter_mut()
    }
}

impl From<u8> for Multivector {
    fn from(u: u8) -> Self {
        (u as f64).into()
    }
}
impl From<f64> for Multivector {
    fn from(f: f64) -> Self {
        Multivector::all([f])
    }
}
impl From<Complex> for Multivector {
    fn from(c: Complex) -> Self {
        Multivector::complex(c.re, c.im)
    }
}

impl fmt::Display for Multivector {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        let dims = self.dims();
        let dim_offset = (self.flavor.metric(0) != 0) as usize;
        let (mask_table, _) = mask_tables(dims);
        let mut wrote = false;
        for (i, n) in self.iter().enumerate() {
            if n == 0.0 {
                continue;
            }
            if wrote {
                if n > 0.0 {
                    write!(f, " + ")?;
                } else {
                    write!(f, " - ")?;
                }
                wrote = true;
            } else if n < 0.0 {
                write!(f, "¯")?;
                wrote = true;
            }
            let mask = mask_table[i];
            if n.abs() != 1.0 || mask == 0 {
                write!(f, "{}", n.abs())?;
                wrote = true;
            }
            if mask == 0 {
                continue;
            }
            write!(f, "e")?;
            wrote = true;
            if dims >= 3 && mask == 0b101 {
                for j in (0..dims).rev() {
                    if mask & (1 << j) != 0 {
                        write!(f, "{}", SUBSCRIPT_DIGITS[j as usize + dim_offset])?;
                    }
                }
            } else {
                for j in 0..dims {
                    if mask & (1 << j) != 0 {
                        write!(f, "{}", SUBSCRIPT_DIGITS[j as usize + dim_offset])?;
                    }
                }
            }
        }
        if !wrote {
            match self.flavor {
                Flavor::Vanilla => write!(f, "0e{}", SUBSCRIPT_DIGITS[dims as usize]),
                Flavor::Projective => write!(f, "0e₀"),
                _ => write!(f, "0"),
            }?
        }
        Ok(())
    }
}

#[cfg(test)]
mod test {
    use super::Multivector as Mv;

    #[test]
    fn conform() {
        let a = Mv::all([1.0, 2.0, 3.0, 4.0]);
        let mut b = a.clone();
        b.set_dims(3);
        assert_eq!(b, [1.0, 2.0, 3.0, 0.0, 4.0, 0.0, 0.0, 0.0]);
        b.set_dims(2);
        assert_eq!(a, b);

        let a = Mv::all([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]);
        let mut b = a.clone();
        b.set_dims(4);
        assert_eq!(
            b,
            [
                1.0, 2.0, 3.0, 4.0, 0.0, 5.0, -6.0, 7.0, 0.0, 0.0, 0.0, 8.0, 0.0, 0.0, 0.0, 0.0
            ]
        );
        b.set_dims(3);
        assert_eq!(a, b);

        let mut a = Mv::all([1.0, 2.0, 3.0, 4.0]);
        a.set_dims_right(3);
        assert_eq!(a, [1.0, 0.0, 2.0, 3.0, 0.0, 0.0, 4.0, 0.0]);

        let mut a = Mv::all([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]);
        a.set_dims_right(4);
        assert_eq!(
            a,
            [
                1.0, 0.0, 2.0, 3.0, 4.0, 0.0, 0.0, 0.0, 5.0, 6.0, 7.0, 0.0, 0.0, 0.0, 8.0, 0.0
            ]
        );
    }

    #[test]
    fn product() {
        assert_eq!(
            Mv::complex(1, 2) * Mv::complex(3, 4),
            [-5.0, 0.0, 0.0, 10.0],
        );
        assert_eq!(
            Mv::vga_vector([1.0, 2.0]) * Mv::vga_vector([3.0, 4.0]),
            [11.0, 0.0, 0.0, -2.0],
        );
        assert_eq!(
            Mv::vga_vector([1.0, 2.0, 3.0]) * Mv::vga_vector([4.0, 5.0, 6.0]),
            [32.0, 0.0, 0.0, 0.0, -3.0, 6.0, -3.0, 0.0],
        );
    }

    #[test]
    fn dual() {
        assert_eq!(
            Mv::all([1.0, 2.0, 3.0, 4.0]).dualed(),
            [-4.0, -3.0, 2.0, 1.0],
        );
        assert_eq!(
            Mv::all([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]).dualed(),
            [-8.0, -7.0, -6.0, -5.0, 4.0, 3.0, 2.0, 1.0],
        );
    }

    #[test]
    fn reverse() {
        assert_eq!(
            Mv::all([1.0, 2.0, 3.0, 4.0]).reversed(),
            [1.0, 2.0, 3.0, -4.0],
        );
        assert_eq!(
            Mv::all([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]).reversed(),
            [1.0, 2.0, 3.0, 4.0, -5.0, -6.0, -7.0, -8.0],
        );
    }

    #[test]
    fn get_blade() {
        let mv = Mv::all([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]);
        assert_eq!(mv.get_blade(0b000), 1.0);
        assert_eq!(mv.get_blade(0b001), 2.0);
        assert_eq!(mv.get_blade(0b010), 3.0);
        assert_eq!(mv.get_blade(0b100), 4.0);
        assert_eq!(mv.get_blade(0b011), 5.0);
        assert_eq!(mv.get_blade(0b101), 6.0);
        assert_eq!(mv.get_blade(0b110), 7.0);
        assert_eq!(mv.get_blade(0b111), 8.0);
        assert_eq!(mv.get_blade(0b1000), 0.0);
        assert_eq!(mv.get_blade(0b1010), 0.0);
    }
}