nsys-math-utils 0.2.0

//! Abstract traits

use either::Either;
#[cfg(feature = "derive_serdes")]
use serde;

use crate::num_traits as num;
use crate::types::{Affinity, LinearIso, Projectivity, Sign};

/// Projective completion (homogeneous coordinates)
pub trait ProjectiveSpace <S : Field> : AffineSpace <S> {
  /// Affine subspace
  type Patch : AffineSpace <S>;
  /// Construct an augmented matrix for the affinity
  #[expect(clippy::type_complexity)]
  fn homography (affinity :
    Affinity <S, Self::Patch, Self::Patch,
      <<Self::Patch as AffineSpace <S>>::Translation as Module <S>>::LinearEndo>
  ) -> Projectivity <S, Self, Self, <Self::Translation as Module <S>>::LinearEndo> where
    <Self::Translation as Module <S>>::LinearEndo :
      From <<<Self::Patch as AffineSpace <S>>::Translation as Module <S>>::LinearEndo>
  {
    Projectivity::new (LinearIso::new ((*affinity.linear_iso).into()).unwrap())
  }
  /// Return the projective completion (homogeneous coordinates) of the affine point or
  /// vector
  fn homogeneous (point_or_vector :
    Either <Self::Patch, <Self::Patch as AffineSpace <S>>::Translation>
  ) -> Self where Self : From <(Self::Patch, S)> {
    match point_or_vector {
      Either::Left  (point)  => (point, S::one()).into(),
      Either::Right (vector) => (Self::Patch::from_vector (vector), S::zero()).into()
    }
  }
}

/// `AffineSpace` with translations in a (Euclidean) real inner product space
pub trait EuclideanSpace <S : Real> : AffineSpace <S> + MetricSpace <S> { }

/// Space of `Point`s (positions) and `Vector`s (displacements)
pub trait AffineSpace <S : Field> : Point <Self::Translation> {
  type Translation : VectorSpace <S> + GroupAction <Self>;
}

/// Point types convertible to and from a vector type, with difference function that
/// follows from the free and transitive group action
pub trait Point <V> : Sized + std::ops::Sub <Self, Output=V> where
  V : AdditiveGroup + GroupAction <Self>
{
  fn to_vector (self) -> V;
  fn from_vector (vector : V) -> Self;
  fn origin() -> Self {
    Self::from_vector (V::zero())
  }
}

/// (Right) action of a group on a set
pub trait GroupAction <X> : Group {
  fn action (self, x : X) -> X;
}

impl <G> GroupAction <G> for G where G : Group {
  #[expect(clippy::renamed_function_params)]
  fn action (self, g : Self) -> Self {
    Self::operation (g, self)
  }
}

/// Bilinear form on a `VectorSpace`
pub trait InnerProductSpace <S : Field> : VectorSpace <S> + Dot <S> {
  fn inner_product (self, other : Self) -> S {
    self.dot (other)
  }
  fn outer_product (self, other : Self) -> Self::LinearEndo;
}

impl <V, S> NormedVectorSpace <S> for V where
  V : InnerProductSpace <S>,
  S : Field
{
  fn norm_squared (self) -> S {
    self.self_dot()
  }
}

/// Scalar product (bilinear form) on a `Module`
pub trait Dot <S : Ring> : Module <S> {
  fn dot (self, other : Self) -> S;
  /// Self dot product
  fn self_dot (self) -> S {
    self.dot (self)
  }
}

/// `VectorSpace` with vector length/magnitude function
pub trait NormedVectorSpace <S : Field> : InnerProductSpace <S> {
  fn norm_squared (self) -> S;
  fn norm (self) -> S where S : Sqrt {
    self.norm_squared().sqrt()
  }
  #[must_use]
  fn normalize (self) -> Self where S : Sqrt {
    self / self.norm()
  }
  fn unit_sigvec (self) -> Self where S : SignedExt + Sqrt {
    let v = self.sigvec();
    if v.is_zero() {
      v
    } else {
      v.normalize()
    }
  }
}

/// Set of points with distance function
pub trait MetricSpace <S : Field> : NormedVectorSpace <S> {
  fn distance_squared (self, other : Self) -> S {
    (self - other).norm_squared()
  }
  fn distance (self, other : Self) -> S where S : Sqrt {
    self.distance_squared (other).sqrt()
  }
}

/// Module with scalars taken from a `Field`
pub trait VectorSpace <S : Field> : Module <S> + std::ops::Div <S, Output=Self> + Copy {
  // required
  fn map <F> (self, f : F) -> Self where F : FnMut (S) -> S;
  // provided
  /// Map `signum_or_zero` over each element of the given vector
  fn sigvec (self) -> Self where S : SignedExt {
    self.map (SignedExt::signum_or_zero)
  }
}

/// Additive group combined with scalar multiplication
pub trait Module <S : Ring> : AdditiveGroup + std::ops::Mul <S, Output=Self> + Copy {
  /// Linear endomorphism represented by a square matrix type
  type LinearEndo : LinearMap <S, Self, Self> + Ring;
}

/// Module homomorphism
pub trait LinearMap <S, V, W> : std::ops::Mul <V, Output=W> + Copy where
  S : Ring,
  V : Module <S>,
  W : Module <S>
{
  fn determinant (self) -> S;
  fn transpose   (self) -> Self;
}

/// Additional matrix methods
pub trait Matrix <S> : Copy {
  type Rows;
  type Submatrix;
  /// Returns the row-major matrix
  fn rows (self) -> Self::Rows;
  /// Returns the submatrix formed by removing the given (col,row)
  fn submatrix (self, i : usize, j : usize) -> Self::Submatrix;
  /// Returns the matrix formed by filling additional row and column with zeros
  fn fill_zeros (submatrix : Self::Submatrix) -> Self where S : num::Zero;
  /// Computes the determinant of the (i,j)-submatrix
  #[inline]
  fn minor <V, W> (self, i : usize, j : usize) -> S where
    S : Ring,
    V : Module <S>,
    W : Module <S>,
    Self::Submatrix : LinearMap <S, V, W>
  {
    self.submatrix (i, j).determinant()
  }
}

/// Householder transformation
pub trait ElementaryReflector {
  type Vector;
  /// Construct a Householder transformation
  fn elementary_reflector (v : Self::Vector, index : usize) -> Self;
}

/// `OrderedField` with special functions
pub trait Real : OrderedField + Exp + Powf + Sqrt + Trig {
  // TODO: more constants
  fn pi() -> Self;
  fn frac_pi_3() -> Self;
  fn sqrt_3() -> Self;
  fn frac_1_sqrt_3() -> Self;
  fn floor (self) -> Self;
  fn ceil (self) -> Self;
  fn trunc (self) -> Self;
  fn fract (self) -> Self;
}

/// A `Field` with an ordering
pub trait OrderedField : Field + OrderedRing { }

/// A (commutative) `Ring` where $1 \neq 0$ and all non-zero elements are invertible
pub trait Field : Ring + MultiplicativeGroup + Powi {
  fn half() -> Self {
    Self::one() / Self::two()
  }
  fn two() -> Self {
    Self::one() + Self::one()
  }
  fn three() -> Self {
    Self::one() + Self::one() + Self::one()
  }
  fn four() -> Self {
    Self::two() * Self::two()
  }
  fn five() -> Self {
    Self::two() + Self::three()
  }
  fn six() -> Self {
    Self::two() * Self::three()
  }
  fn seven() -> Self {
    Self::three() + Self::four()
  }
  fn eight() -> Self {
    Self::two() * Self::two() * Self::two()
  }
  fn nine() -> Self {
    Self::three() * Self::three()
  }
  fn ten() -> Self {
    Self::two() * Self::five()
  }
}

/// Interface for a group with identity represented by `one`, operation defined by `*`
/// and `/`
pub trait MultiplicativeGroup : MultiplicativeMonoid +
  std::ops::Div <Self, Output=Self> + std::ops::DivAssign + num::Inv <Output=Self>
{ }

/// Ring of integers
pub trait Integer : OrderedRing + SignedExt + num::PrimInt { }

/// A Ring with an ordering
pub trait OrderedRing : Ring + SignedExt + MinMax + PartialOrd + std::ops::Rem { }

/// Interface for a ring as a combination of an additive group and a distributive
/// multiplication operation
pub trait Ring : AdditiveGroup + MultiplicativeMonoid
  + num::MulAdd <Self, Self, Output=Self> + num::MulAddAssign <Self, Self>
{ }

/// Interface for a group with identity represented by `zero`, and operation defined by
/// `+` and `-`
pub trait AdditiveGroup : AdditiveMonoid +
  std::ops::Sub <Self, Output=Self> + std::ops::SubAssign + std::ops::Neg <Output=Self>
{ }

/// Set with identity, inverses, and (associative) binary operation
pub trait Group : PartialEq + std::ops::Neg <Output=Self> {
  fn identity() -> Self;
  fn operation (a : Self, b : Self) -> Self;
}

/// Set with identity represented by `one` and (associative) binary operation defined by
/// `*`
pub trait MultiplicativeMonoid : Copy + PartialEq +
  std::ops::Mul <Self, Output=Self> + std::ops::MulAssign + num::One
{
  fn squared (self) -> Self {
    self * self
  }
  fn cubed (self) -> Self {
    self * self * self
  }
}

/// Set with identity represented by `zero` and (associative) binary operation defined
/// by `+`
pub trait AdditiveMonoid : Sized + PartialEq + std::iter::Sum
  + std::ops::Add <Self, Output=Self> + std::ops::AddAssign + num::Zero
{ }

/// Interface for angle units
pub trait Angle <S : Real> : Clone + Copy + PartialEq + PartialOrd + Sized +
  AdditiveGroup + std::ops::Div <Self, Output=S> + std::ops::Mul <S, Output=Self> +
  std::ops::Div <S, Output=Self> + std::ops::Rem <Self, Output=Self>
{
  // required
  /// Full rotation
  fn full_turn() -> Self;
  // provided
  /// Half rotation
  fn half_turn() -> Self {
    Self::full_turn() / S::two()
  }
  /// Restrict to `(-half_turn, half_turn]`
  fn wrap_signed (self) -> Self {
    if self > Self::half_turn() || self <= -Self::half_turn() {
      let out = (self + Self::half_turn()).wrap_unsigned() - Self::half_turn();
      if out == -Self::half_turn() {
        Self::half_turn()
      } else {
        out
      }
    } else {
      self
    }
  }
  /// Restrict to `[0, full_turn)`
  fn wrap_unsigned (self) -> Self {
    if self >= Self::full_turn() {
      self % Self::full_turn()
    } else if self < Self::zero() {
      self + Self::full_turn() * ((self / Self::full_turn()).trunc().abs() + S::one())
    } else {
      self
    }
  }
}

/// Unsigned integer power function
pub trait Pow {
  fn pow (self, exp : u32) -> Self;
}
/// Signed integer power function
pub trait Powi {
  fn powi (self, n : i32) -> Self;
}
/// Fractional power function
pub trait Powf {
  fn powf (self, n : Self) -> Self;
}
/// Exponential function
pub trait Exp {
  fn exp (self) -> Self;
}
/// Square root function
pub trait Sqrt {
  fn sqrt (self) -> Self;
}
/// Cube root function
pub trait Cbrt {
  fn cbrt (self) -> Self;
}
/// Trigonometric functions
pub trait Trig : Sized {
  fn sin      (self) -> Self;
  fn sin_cos  (self) -> (Self, Self);
  fn cos      (self) -> Self;
  fn tan      (self) -> Self;
  fn asin     (self) -> Self;
  fn acos     (self) -> Self;
  fn atan     (self) -> Self;
  fn atan2    (self, other : Self) -> Self;
}

/// Provides `min`, `max`, and `clamp` that are not necessarily consistent with those
/// from `Ord`. This is provided because `f32` and `f64` do not implement `Ord`, so this
/// trait is defined to give a uniform interface with `Ord` types.
pub trait MinMax {
  fn min   (self, other : Self) -> Self;
  fn max   (self, other : Self) -> Self;
  fn clamp (self, min : Self, max : Self) -> Self;
}

/// Function returning number representing sign of self
pub trait SignedExt : num::Signed {
  #[inline]
  fn sign (self) -> Sign {
    if self.is_zero() {
      Sign::Zero
    } else if self.is_positive() {
      Sign::Positive
    } else {
      debug_assert!(self.is_negative());
      Sign::Negative
    }
  }
  /// Maps `0.0` to `0.0`, otherwise equal to `S::signum` (which would otherwise map
  /// `+0.0 -> 1.0` and `-0.0 -> -1.0`)
  #[inline]
  fn signum_or_zero (self) -> Self where Self : num::Zero {
    if self.is_zero() {
      Self::zero()
    } else {
      self.signum()
    }
  }
}

/// Adds serde `Serialize` and `DeserializeOwned` constraints.
///
/// This makes it easier to conditionally add these constraints to type definitions when
/// `derive_serdes` feature is enabled.
#[cfg(not(feature = "derive_serdes"))]
pub trait MaybeSerDes { }
#[cfg(feature = "derive_serdes")]
pub trait MaybeSerDes : serde::Serialize + serde::de::DeserializeOwned { }

impl <S, T> MetricSpace <S> for T where S : Field, T : NormedVectorSpace <S> { }
impl <T> OrderedField for T where T : Field + OrderedRing { }
impl <T> Field        for T where T : Ring + MultiplicativeGroup + Powi + Powf { }
impl <T> Integer      for T where T : OrderedRing + num::PrimInt { }
impl <T> OrderedRing  for T where
  T : Ring + SignedExt + MinMax + PartialOrd + std::ops::Rem
{ }
impl <T> Ring         for T where
  T : AdditiveGroup + MultiplicativeMonoid
    + num::MulAdd<Self, Self, Output=Self> + num::MulAddAssign <Self, Self>
{ }
impl <T> AdditiveGroup for T where
  T : AdditiveMonoid + std::ops::Sub <Self, Output=Self> + std::ops::SubAssign
    + std::ops::Neg <Output=Self>
{ }
impl <T> MultiplicativeGroup for T where
  T : MultiplicativeMonoid + std::ops::Div <Self, Output=Self> + std::ops::DivAssign
    + num::Inv <Output=Self>
{ }
impl <T> AdditiveMonoid for T where
  T : Sized + PartialEq + std::iter::Sum + std::ops::Add <Self, Output=Self>
    + std::ops::AddAssign + num::Zero
{ }
impl <T> MultiplicativeMonoid for T where
  T : Copy + PartialEq + std::ops::Mul <Self, Output=Self> + std::ops::MulAssign
    + num::One
{ }
impl <T : num::Signed> SignedExt for T { }

impl <
  #[cfg(not(feature = "derive_serdes"))]
  T,
  #[cfg(feature = "derive_serdes")]
  T : serde::Serialize + serde::de::DeserializeOwned
> MaybeSerDes for T { }

macro impl_integer ($type:ty) {
  impl MinMax         for $type {
    fn min (self, other : Self) -> Self {
      Ord::min (self, other)
    }
    fn max (self, other : Self) -> Self {
      Ord::max (self, other)
    }
    fn clamp (self, min : Self, max : Self) -> Self {
      Ord::clamp (self, min, max)
    }
  }
  impl Pow for $type {
    fn pow (self, exp : u32) -> Self {
      self.pow (exp)
    }
  }
}
impl_integer!(i8);
impl_integer!(i16);
impl_integer!(i32);
impl_integer!(i64);

macro impl_real_float ($type:ident) {
  impl VectorSpace <Self> for $type {
    fn map <F> (self, mut f : F) -> Self where F : FnMut (Self) -> Self {
      f (self)
    }
  }
  impl Module <Self> for $type {
    type LinearEndo = Self;
  }
  impl LinearMap <Self, Self, Self> for $type {
    fn determinant (self) -> Self {
      self
    }
    fn transpose (self) -> Self {
      self
    }
  }
  impl Real for $type {
    fn pi() -> Self {
      std::$type::consts::PI
    }
    fn frac_pi_3() -> Self {
      std::$type::consts::FRAC_PI_3
    }
    #[expect(clippy::excessive_precision)]
    #[allow(clippy::allow_attributes)]
    #[allow(clippy::cast_possible_truncation)]
    #[allow(trivial_numeric_casts)]
    fn sqrt_3() -> Self {
      1.732050807568877293527446341505872366942805253810380628055f64 as $type
    }
    #[expect(clippy::excessive_precision)]
    #[allow(clippy::allow_attributes)]
    #[allow(clippy::cast_possible_truncation)]
    #[allow(trivial_numeric_casts)]
    fn frac_1_sqrt_3() -> Self {
      (1.0f64 / 1.732050807568877293527446341505872366942805253810380628055f64) as $type
    }
    fn floor (self) -> Self {
      self.floor()
    }
    fn ceil (self) -> Self {
      self.ceil()
    }
    fn trunc (self) -> Self {
      self.trunc()
    }
    fn fract (self) -> Self {
      self.fract()
    }
  }
  impl MinMax               for $type {
    fn min (self, other : Self) -> Self {
      self.min (other)
    }
    fn max (self, other : Self) -> Self {
      self.max (other)
    }
    fn clamp (self, min : Self, max : Self) -> Self {
      self.clamp (min, max)
    }
  }
  impl Powi for $type {
    fn powi (self, n : i32) -> Self {
      self.powi (n)
    }
  }
  impl Powf for $type {
    fn powf (self, n : Self) -> Self {
      self.powf (n)
    }
  }
  impl Exp for $type {
    fn exp (self) -> Self {
      self.exp()
    }
  }
  impl Sqrt for $type {
    fn sqrt (self) -> Self {
      self.sqrt()
    }
  }
  impl Cbrt for $type {
    fn cbrt (self) -> Self {
      self.cbrt()
    }
  }
  impl Trig for $type {
    fn sin (self) -> Self {
      self.sin()
    }
    fn sin_cos (self) -> (Self, Self) {
      self.sin_cos()
    }
    fn cos (self) -> Self {
      self.cos()
    }
    fn tan (self) -> Self {
      self.tan()
    }
    fn asin (self) -> Self {
      self.asin()
    }
    fn acos (self) -> Self {
      self.acos()
    }
    fn atan (self) -> Self {
      self.atan()
    }
    fn atan2 (self, other : Self) -> Self {
      self.atan2 (other)
    }
  }
}

impl_real_float!(f32);
impl_real_float!(f64);