vecmat 0.7.9

//! Quaternion pairs with a central complex, dual, or split unit.

use super::Quaternion;
use crate::traits::Dot;
use core::ops::{Add, Div, Mul, Neg, Sub};
use num_traits::{Float, Num, One, Zero};

/// An element `a + e*b`, where `a` and `b` are ordinary quaternions and
/// the central unit satisfies `e*e = K`.
///
/// `K` must be `-1` (complex), `0` (dual), or `1` (split). Different values
/// are different types and cannot accidentally be multiplied together.
/// Both quaternions use scalar-first component order. This is an algebra
/// element, not necessarily a normalized isometry.
///
/// Quaternion conjugation, scalar conjugation (`e -> -e`), and their
/// combination are distinct operations. The quaternion norm is itself a
/// central pair, so ordinary eight-dimensional Euclidean normalization is
/// not appropriate.
#[repr(C)]
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct QuaternionPair<T, const K: i8> {
    real: Quaternion<T>,
    unreal: Quaternion<T>,
}

impl<T, const K: i8> QuaternionPair<T, K> {
    /// Construct an algebra element. Panics if `K` is not `-1`, `0`, or `1`.
    pub fn new(real: Quaternion<T>, unreal: Quaternion<T>) -> Self {
        assert!(
            (-1..=1).contains(&K),
            "quaternion-pair K must be -1, 0, or 1"
        );
        Self { real, unreal }
    }

    pub fn into_parts(self) -> (Quaternion<T>, Quaternion<T>) {
        (self.real, self.unreal)
    }

    /// Construct from `[a.w, a.x, a.y, a.z, b.w, b.x, b.y, b.z]`.
    pub fn from_array(values: [T; 8]) -> Self {
        let [aw, ax, ay, az, bw, bx, by, bz] = values;
        Self::new(
            Quaternion::new(aw, ax, ay, az),
            Quaternion::new(bw, bx, by, bz),
        )
    }

    /// Return `[a.w, a.x, a.y, a.z, b.w, b.x, b.y, b.z]`.
    pub fn into_array(self) -> [T; 8] {
        let [aw, ax, ay, az] = self.real.into_array();
        let [bw, bx, by, bz] = self.unreal.into_array();
        [aw, ax, ay, az, bw, bx, by, bz]
    }
}

impl<T: Copy, const K: i8> QuaternionPair<T, K> {
    pub fn real(&self) -> Quaternion<T> {
        self.real
    }

    pub fn unreal(&self) -> Quaternion<T> {
        self.unreal
    }
}

impl<T: Neg<Output = T>, const K: i8> QuaternionPair<T, K> {
    /// Quaternion conjugation: `(conj(a), conj(b))`. Reverses products.
    pub fn quaternion_conj(self) -> Self {
        Self::new(self.real.conj(), self.unreal.conj())
    }

    /// Conjugation of the central unit: `(a, -b)`. Preserves product order.
    pub fn scalar_conj(self) -> Self {
        Self::new(self.real, -self.unreal)
    }

    /// Combined conjugation: `(conj(a), -conj(b))`. Reverses products.
    pub fn combined_conj(self) -> Self {
        self.quaternion_conj().scalar_conj()
    }
}

impl<T: Neg<Output = T>, const K: i8> Neg for QuaternionPair<T, K> {
    type Output = Self;
    fn neg(self) -> Self {
        Self::new(-self.real, -self.unreal)
    }
}

impl<T: Add<Output = T>, const K: i8> Add for QuaternionPair<T, K> {
    type Output = Self;
    fn add(self, other: Self) -> Self {
        Self::new(self.real + other.real, self.unreal + other.unreal)
    }
}

impl<T: Sub<Output = T>, const K: i8> Sub for QuaternionPair<T, K> {
    type Output = Self;
    fn sub(self, other: Self) -> Self {
        Self::new(self.real - other.real, self.unreal - other.unreal)
    }
}

impl<T: Num + Copy, const K: i8> Mul for QuaternionPair<T, K> {
    type Output = Self;
    fn mul(self, other: Self) -> Self {
        let ac = self.real * other.real;
        let real = match K {
            -1 => ac - self.unreal * other.unreal,
            0 => ac,
            1 => ac + self.unreal * other.unreal,
            _ => unreachable!("validated quaternion-pair K"),
        };
        Self::new(real, self.real * other.unreal + self.unreal * other.real)
    }
}

impl<T: Mul<Output = T> + Copy, const K: i8> Mul<T> for QuaternionPair<T, K> {
    type Output = Self;
    fn mul(self, scalar: T) -> Self {
        Self::new(self.real * scalar, self.unreal * scalar)
    }
}

impl<T: Div<Output = T> + Copy, const K: i8> Div<T> for QuaternionPair<T, K> {
    type Output = Self;
    fn div(self, scalar: T) -> Self {
        Self::new(self.real / scalar, self.unreal / scalar)
    }
}

impl<T: Zero, const K: i8> Zero for QuaternionPair<T, K> {
    fn zero() -> Self {
        Self::new(Quaternion::zero(), Quaternion::zero())
    }
    fn is_zero(&self) -> bool {
        self.real.is_zero() && self.unreal.is_zero()
    }
}

impl<T: Zero, const K: i8> Default for QuaternionPair<T, K> {
    fn default() -> Self {
        Self::zero()
    }
}

impl<T: Num + Copy, const K: i8> One for QuaternionPair<T, K> {
    fn one() -> Self {
        Self::new(Quaternion::one(), Quaternion::zero())
    }
}

impl<T: Num + Copy, const K: i8> QuaternionPair<T, K> {
    /// The central coefficients `(n, m)` of `self * self.quaternion_conj()`:
    /// `n = |a|² + K |b|²`, `m = 2 dot(a,b)`.
    pub fn norm_sqr(self) -> (T, T) {
        let aa = self.real.norm_sqr();
        let n = match K {
            -1 => aa - self.unreal.norm_sqr(),
            0 => aa,
            1 => aa + self.unreal.norm_sqr(),
            _ => unreachable!("validated quaternion-pair K"),
        };
        let ab = self.real.dot(self.unreal);
        (n, ab + ab)
    }
}

impl<T: Float, const K: i8> QuaternionPair<T, K> {
    pub fn is_finite(self) -> bool {
        self.into_array().into_iter().all(Float::is_finite)
    }

    /// Invert a nonsingular element. Returns `None` for singular/nonfinite
    /// inputs or when intermediate arithmetic overflows.
    pub fn try_inv(self) -> Option<Self> {
        if !self.is_finite() {
            return None;
        }
        let (n, m) = self.norm_sqr();
        let det = match K {
            -1 => n * n + m * m,
            0 => n * n,
            1 => n * n - m * m,
            _ => unreachable!("validated quaternion-pair K"),
        };
        if det == T::zero() || !det.is_finite() {
            return None;
        }
        let zero = T::zero();
        let inv_norm = Self::new(
            Quaternion::new(n / det, zero, zero, zero),
            Quaternion::new(-m / det, zero, zero, zero),
        );
        let inverse = self.quaternion_conj() * inv_norm;
        inverse.is_finite().then_some(inverse)
    }

    /// Divide by a central square root of the quaternion norm.
    ///
    /// The result has `norm_sqr() = (1, 0)` up to roundoff. Complex pairs
    /// use the principal complex square root; dual and split pairs use
    /// the positive real branches. Returns `None` for singular/nonfinite
    /// inputs or when intermediate arithmetic overflows. This selects a
    /// normalization branch, not a canonical sign for the double cover.
    pub fn try_normalize(self) -> Option<Self> {
        if !self.is_finite() {
            return None;
        }
        let (n, m) = self.norm_sqr();
        let one = T::one();
        let two = one + one;
        let (r, s) = match K {
            -1 => {
                let root = num_complex::Complex::new(n, m).sqrt();
                let det = root.norm_sqr();
                (root.re / det, -root.im / det)
            }
            0 => {
                let r = one / n.sqrt();
                (r, -m * r / (two * n))
            }
            1 => {
                let plus = one / (n + m).sqrt();
                let minus = one / (n - m).sqrt();
                ((plus + minus) / two, (plus - minus) / two)
            }
            _ => unreachable!("validated quaternion-pair K"),
        };
        if !r.is_finite() || !s.is_finite() {
            return None;
        }
        let zero = T::zero();
        let inv_root = Self::new(
            Quaternion::new(r, zero, zero, zero),
            Quaternion::new(s, zero, zero, zero),
        );
        let normalized = self * inv_root;
        normalized.is_finite().then_some(normalized)
    }
}

#[cfg(test)]
#[path = "tests/quaternion_pair.rs"]
mod tests;