nsys-math-utils 0.2.0

use std::marker::PhantomData;
use derive_more::{Deref, Display, From, Neg};

use crate::approx;
use crate::num_traits as num;

#[cfg(feature = "derive_serdes")]
use serde::{Deserialize, Serialize};

use crate::traits::*;

pub mod coordinate;

pub use vek::Vec2 as Vector2;
pub use vek::Vec3 as Vector3;
pub use vek::Vec4 as Vector4;
pub use vek::Mat2 as Matrix2;
pub use vek::Mat3 as Matrix3;
pub use vek::Mat4 as Matrix4;
pub use vek::Quaternion as Quaternion;

pub use self::enums::{Axis2, Axis3, Axis4, Octant, Quadrant, Sign, SignedAxis1,
  SignedAxis2, SignedAxis3, SignedAxis4};

////////////////////////////////////////////////////////////////////////////////////////
//  enums                                                                             //
////////////////////////////////////////////////////////////////////////////////////////

// we define enums in a separate module so that the generated iterator types don't
// pollute the namespace
mod enums {
  use strum::{EnumCount, EnumIter, FromRepr};
  #[cfg(feature = "derive_serdes")]
  use serde::{Deserialize, Serialize};

  /// Negative, zero, or positive
  #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
  #[derive(Clone, Copy, Debug, Eq, PartialEq, EnumCount, EnumIter, FromRepr)]
  #[repr(u8)]
  pub enum Sign {
    Negative,
    Zero,
    Positive
  }

  /// 2D quadrant
  #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
  #[derive(Clone, Copy, Debug, Eq, PartialEq, EnumCount, EnumIter, FromRepr)]
  #[repr(u8)]
  pub enum Quadrant {
    /// +X, +Y
    PosPos,
    /// -X, +Y
    NegPos,
    /// +X, -Y
    PosNeg,
    /// -X, -Y
    NegNeg
  }

  /// 3D octant
  #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
  #[derive(Clone, Copy, Debug, Eq, PartialEq, EnumCount, EnumIter, FromRepr)]
  #[repr(u8)]
  pub enum Octant {
    /// +X, +Y, +Z
    PosPosPos,
    /// -X, +Y, +Z
    NegPosPos,
    /// +X, -Y, +Z
    PosNegPos,
    /// -X, -Y, +Z
    NegNegPos,
    /// +X, +Y, -Z
    PosPosNeg,
    /// -X, +Y, -Z
    NegPosNeg,
    /// +X, -Y, -Z
    PosNegNeg,
    /// -X, -Y, -Z
    NegNegNeg
  }

  /// X or Y axis
  #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
  #[derive(Clone, Copy, Debug, Eq, PartialEq, EnumCount, EnumIter, FromRepr)]
  #[repr(u8)]
  pub enum Axis2 {
    X=0,
    Y
  }

  /// X, Y, or Z axis
  #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
  #[derive(Clone, Copy, Debug, Eq, PartialEq, EnumCount, EnumIter, FromRepr)]
  #[repr(u8)]
  pub enum Axis3 {
    X=0,
    Y,
    Z
  }

  /// X, Y, Z, or W axis
  #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
  #[derive(Clone, Copy, Debug, Eq, PartialEq, EnumCount, EnumIter, FromRepr)]
  #[repr(u8)]
  pub enum Axis4 {
    X=0,
    Y,
    Z,
    W
  }

  /// Positive or negative X axis
  #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
  #[derive(Clone, Copy, Debug, Eq, PartialEq, EnumCount, EnumIter, FromRepr)]
  #[repr(u8)]
  pub enum SignedAxis1 {
    NegX, PosX
  }

  /// Positive or negative X or Y axis
  #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
  #[derive(Clone, Copy, Debug, Eq, PartialEq, EnumCount, EnumIter, FromRepr)]
  #[repr(u8)]
  pub enum SignedAxis2 {
    NegX, PosX,
    NegY, PosY
  }

  /// Positive or negative X, Y, or Z axis
  #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
  #[derive(Clone, Copy, Debug, Eq, PartialEq, EnumCount, EnumIter, FromRepr)]
  #[repr(u8)]
  pub enum SignedAxis3 {
    NegX, PosX,
    NegY, PosY,
    NegZ, PosZ
  }

  /// Positive or negative X, Y, Z, or W axis
  #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
  #[derive(Clone, Copy, Debug, Eq, PartialEq, EnumCount, EnumIter, FromRepr)]
  #[repr(u8)]
  pub enum SignedAxis4 {
    NegX, PosX,
    NegY, PosY,
    NegZ, PosZ,
    NegW, PosW
  }
}

////////////////////////////////////////////////////////////////////////////////////////
//  structs                                                                           //
////////////////////////////////////////////////////////////////////////////////////////

/// Strictly positive scalars (non-zero)
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Default, Eq, PartialEq, PartialOrd)]
pub struct Positive <S> (S);

/// Non-negative scalars (may be zero)
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Default, Eq, PartialEq, PartialOrd)]
pub struct NonNegative <S> (S);

/// Non-zero scalars
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Default, PartialEq, PartialOrd)]
pub struct NonZero <S> (S);

/// Scalars in the closed unit interval [0, 1]
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Default, PartialEq, PartialOrd)]
pub struct Normalized <S> (S);

/// Scalars in the closed interval [-1, 1]
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Default, PartialEq, PartialOrd)]
pub struct NormalSigned <S> (S);

/// (Euler angles) representation of a 3D orientation. Internally the angles are wrapped
/// to $[0, 2\pi)$.
///
/// Throughout this library this is usually interpreted as intrinsic ZX'Y'' Euler
/// angles.
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Default, Eq, PartialEq)]
pub struct Angles3 <S> {
  pub yaw   : AngleWrapped <S>,
  pub pitch : AngleWrapped <S>,
  pub roll  : AngleWrapped <S>
}

/// Representation of a 3D position + orientation
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Default, Eq, PartialEq)]
pub struct Pose3 <S> {
  pub position : Point3  <S>,
  pub angles   : Angles3 <S>
}

/// Invertible linear map
#[derive(Clone, Copy, Debug, Default, Eq, PartialEq, Display)]
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[display("{}", _0)]
pub struct LinearIso <S, V, W, M> (M, PhantomData <(S, V, W)>);

/// Invertible linear endomorphism
#[derive(Clone, Copy, Debug, Default, PartialEq, Deref, Display, From)]
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[display("{}", _0)]
pub struct LinearAuto <S, V> (pub LinearIso <S, V, V, V::LinearEndo>) where
  V : Module <S>,
  V::LinearEndo : MaybeSerDes,
  S : Ring;

/// Orthogonal 2x2 matrix with determinant +1, i.e. a member of the circle group $SO(2)$
/// of special orthogonal matrices
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[cfg_attr(feature = "derive_serdes", serde(bound = "S : MaybeSerDes"))]
#[derive(Clone, Copy, Debug, Default, PartialEq)]
pub struct Rotation2 <S> (LinearAuto <S, Vector2 <S>>) where S : Ring + MaybeSerDes;

/// Orthogonal 3x3 matrix with determinant +1, i.e. a member of the 3D rotation group
/// $SO(3)$ of special orthogonal matrices
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[cfg_attr(feature = "derive_serdes", serde(bound = "S : MaybeSerDes"))]
#[derive(Clone, Copy, Debug, Default, PartialEq)]
pub struct Rotation3 <S> (LinearAuto <S, Vector3 <S>>) where S : Ring + MaybeSerDes;

/// Unit quaternion representing an orientation in $\mathbb{R}^3$
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Default, Eq, PartialEq)]
pub struct Versor <S> (Quaternion <S>) where S : num::Zero + num::One;

/// Homomorphism of affine spaces: a combination of a linear map and a translation
#[derive(Clone, Copy, Debug, Default, Eq, PartialEq, Display)]
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[display("({}, {})", linear_map, translation)]
pub struct AffineMap <S, A, B, M> where
  A : AffineSpace <S>,
  B : AffineSpace <S>,
  S : Field,
  M : LinearMap <S, A::Translation, B::Translation>
{
  pub linear_map  : M,
  pub translation : B::Translation,
  pub _phantom    : PhantomData <(A, S)>
}

/// Affine isomorphism
#[derive(Clone, Copy, Debug, Default, Eq, PartialEq, Display)]
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[display("({}, {})", linear_iso, translation)]
pub struct Affinity <S, A, B, M> where
  S : Field,
  A : AffineSpace <S>,
  B : AffineSpace <S>,
  M : LinearMap <S, A::Translation, B::Translation>
{
  pub linear_iso  : LinearIso <S, A::Translation, B::Translation, M>,
  pub translation : B::Translation
}

/// Isomorphism of projective spaces (a.k.a. homography or projective collineation)
#[derive(Clone, Copy, Debug, Default, Eq, PartialEq, Display)]
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[display("{}", _0)]
pub struct Projectivity <S, P, Q, M> (
  pub LinearIso <S, P::Translation, Q::Translation, M>
) where
  S : Field,
  P : ProjectiveSpace <S>,
  Q : ProjectiveSpace <S>;

////////////////////////////////////////////////////////////////////////////////////////
//  impls                                                                             //
////////////////////////////////////////////////////////////////////////////////////////

//
//  impl Deg
//
impl <S : Real> Angle <S> for Deg <S> {
  fn full_turn() -> Self {
    Deg (S::ten() * S::nine() * S::two() * S::two())
  }
  fn half_turn() -> Self {
    Deg (S::ten() * S::nine() * S::two())
  }
}
impl <S : Real> Trig for Deg <S> {
  fn sin     (self) -> Self {
    Rad::from (self).sin().into()
  }
  fn sin_cos (self) -> (Self, Self) {
    let (sin, cos) = Rad::from (self).sin_cos();
    (sin.into(), cos.into())
  }
  fn cos     (self) -> Self {
    Rad::from (self).cos().into()
  }
  fn tan     (self) -> Self {
    Rad::from (self).tan().into()
  }
  fn asin    (self) -> Self {
    Rad::from (self).asin().into()
  }
  fn acos    (self) -> Self {
    Rad::from (self).acos().into()
  }
  fn atan    (self) -> Self {
    Rad::from (self).atan().into()
  }
  fn atan2   (self, other : Self) -> Self {
    Rad::from (self).atan2 (Rad::from (other)).into()
  }
}
impl <S : Real> From <Rad <S>> for Deg <S> {
  fn from (radians : Rad <S>) -> Self {
    let full_turns = radians / Rad::full_turn();
    Deg::full_turn() * full_turns
  }
}
impl <S : Real> From <Turn <S>> for Deg <S> {
  fn from (turns : Turn <S>) -> Self {
    Deg (turns.0 * Deg::full_turn().0)
  }
}

//
//  impl Rad
//
impl <S : Real> Angle <S> for Rad <S> {
  fn full_turn() -> Self {
    Rad (S::pi() * S::two())
  }
}
impl <S : Real> Trig for Rad <S> {
  fn sin     (self) -> Self {
    Rad (self.0.sin())
  }
  fn sin_cos (self) -> (Self, Self) {
    let (sin, cos) = self.0.sin_cos();
    (Rad (sin), Rad (cos))
  }
  fn cos     (self) -> Self {
    Rad (self.0.cos())
  }
  fn tan     (self) -> Self {
    Rad (self.0.tan())
  }
  fn asin    (self) -> Self {
    Rad (self.0.asin())
  }
  fn acos    (self) -> Self {
    Rad (self.0.acos())
  }
  fn atan    (self) -> Self {
    Rad (self.0.atan())
  }
  fn atan2   (self, other : Self) -> Self {
    Rad (self.0.atan2 (other.0))
  }
}
impl <S : Real> From <Deg <S>> for Rad <S> {
  fn from (degrees : Deg <S>) -> Self {
    let full_turns = degrees / Deg::full_turn();
    Rad::full_turn() * full_turns
  }
}
impl <S : Real> From <Turn <S>> for Rad <S> {
  fn from (turns : Turn <S>) -> Self {
    Rad (turns.0 * Rad::full_turn().0)
  }
}

//
//  impl Turn
//
impl <S : Real> Angle <S> for Turn <S> {
  fn full_turn() -> Self {
    Turn (S::one())
  }
  fn half_turn() -> Self {
    Turn (S::half())
  }
}
impl <S : Real> Trig for Turn <S> {
  fn sin     (self) -> Self {
    Rad::from (self).sin().into()
  }
  fn sin_cos (self) -> (Self, Self) {
    let (sin, cos) = Rad::from (self).sin_cos();
    (sin.into(), cos.into())
  }
  fn cos     (self) -> Self {
    Rad::from (self).cos().into()
  }
  fn tan     (self) -> Self {
    Rad::from (self).tan().into()
  }
  fn asin    (self) -> Self {
    Rad::from (self).asin().into()
  }
  fn acos    (self) -> Self {
    Rad::from (self).acos().into()
  }
  fn atan    (self) -> Self {
    Rad::from (self).atan().into()
  }
  fn atan2   (self, other : Self) -> Self {
    Rad::from (self).atan2 (Rad::from (other)).into()
  }
}
impl <S : Real> From <Rad <S>> for Turn <S> {
  fn from (radians : Rad <S>) -> Self {
    Turn (radians.0 / Rad::full_turn().0)
  }
}
impl <S : Real> From <Deg <S>> for Turn <S> {
  fn from (degrees : Deg <S>) -> Self {
    Turn (degrees.0 / Deg::full_turn().0)
  }
}

//
//  impl AngleWrapped
//
impl <S : Real> AngleWrapped <S> {
  pub fn wrap (angle : Rad <S>) -> Self {
    AngleWrapped (angle.wrap_unsigned())
  }

  pub fn map <F> (self, f : F) -> Self where
    F : FnOnce (Rad <S>) -> Rad <S>
  {
    AngleWrapped (f (self.0).wrap_unsigned())
  }
}

//
//  impl AngleWrappedSigned
//
impl <S : Real> AngleWrappedSigned <S> {
  pub fn wrap (angle : Rad <S>) -> Self {
    AngleWrappedSigned (angle.wrap_signed())
  }

  pub fn map <F> (self, f : F) -> Self where
    F : FnOnce (Rad <S>) -> Rad <S>
  {
    AngleWrappedSigned (f (self.0).wrap_signed())
  }
}

//
//  impl Angles3
//

impl <S : Real> Angles3 <S> {
  pub const fn new (
    yaw   : AngleWrapped <S>,
    pitch : AngleWrapped <S>,
    roll  : AngleWrapped <S>
  ) -> Self {
    Angles3 { yaw, pitch, roll }
  }

  pub fn zero() -> Self {
    use num::Zero;
    Angles3::new (
      AngleWrapped::zero(),
      AngleWrapped::zero(),
      AngleWrapped::zero())
  }

  pub fn wrap (yaw : Rad <S>, pitch : Rad <S>, roll : Rad <S>) -> Self {
    let yaw   = AngleWrapped::wrap (yaw);
    let pitch = AngleWrapped::wrap (pitch);
    let roll  = AngleWrapped::wrap (roll);
    Angles3 { yaw, pitch, roll }
  }
}

impl <S> From <Rotation3 <S>> for Angles3 <S> where S : Real + MaybeSerDes {
  fn from (rotation : Rotation3 <S>) -> Self {
    let (yaw, pitch, roll) = {
      let (yaw, pitch, roll) = rotation.intrinsic_angles();
      ( AngleWrapped::wrap (yaw),
        AngleWrapped::wrap (pitch),
        AngleWrapped::wrap (roll)
      )
    };
    Angles3 { yaw, pitch, roll }
  }
}

//
//  impl Axis2
//
impl Axis2 {
  #[inline]
  pub const fn component (self) -> usize {
    self as usize
  }
}

//
//  impl Axis3
//
impl Axis3 {
  #[inline]
  pub const fn component (self) -> usize {
    self as usize
  }
}

//
//  impl SignedAxis3
//
impl SignedAxis3 {
  pub fn try_from <S : Field> (vector : &Vector3 <S>) -> Option <Self> {
    if *vector == [ S::one(),   S::zero(),  S::zero()].into() {
      Some (SignedAxis3::PosX)
    } else if *vector == [-S::one(),   S::zero(),  S::zero()].into() {
      Some (SignedAxis3::NegX)
    } else if *vector == [ S::zero(),  S::one(),   S::zero()].into() {
      Some (SignedAxis3::PosY)
    } else if *vector == [ S::zero(), -S::one(),   S::zero()].into() {
      Some (SignedAxis3::NegY)
    } else if *vector == [ S::zero(),  S::zero(),  S::one() ].into() {
      Some (SignedAxis3::PosZ)
    } else if *vector == [ S::zero(),  S::zero(), -S::one() ].into() {
      Some (SignedAxis3::NegZ)
    } else {
      None
    }
  }

  pub fn to_vec <S : Field> (self) -> Vector3 <S> {
    match self {
      SignedAxis3::NegX => [-S::one(),   S::zero(),  S::zero()].into(),
      SignedAxis3::PosX => [ S::one(),   S::zero(),  S::zero()].into(),
      SignedAxis3::NegY => [ S::zero(), -S::one(),   S::zero()].into(),
      SignedAxis3::PosY => [ S::zero(),  S::one(),   S::zero()].into(),
      SignedAxis3::NegZ => [ S::zero(),  S::zero(), -S::one() ].into(),
      SignedAxis3::PosZ => [ S::zero(),  S::zero(),  S::one() ].into()
    }
  }
  #[inline]
  pub fn to_unit <S> (self) -> Unit3 <S> where S : Real + std::fmt::Debug {
    Unit3::unchecked (self.to_vec())
  }
}

//
//  impl Quadrant
//
impl Quadrant {
  /// Returns some 'Some (quadrant)' if the vector is strictly contained within a
  /// quadrant and returns 'None' otherwise
  pub fn from_vec_strict <S> (vector : &Vector2 <S>) -> Option <Quadrant> where
    S : Ring + SignedExt
  {
    let sign_x = vector.x.sign();
    let sign_y = vector.y.sign();
    let quadrant = match (sign_x, sign_y) {
      (Sign::Zero,          _) |
      (         _, Sign::Zero) => return None,
      (Sign::Positive, Sign::Positive) => Quadrant::PosPos,
      (Sign::Negative, Sign::Positive) => Quadrant::NegPos,
      (Sign::Positive, Sign::Negative) => Quadrant::PosNeg,
      (Sign::Negative, Sign::Negative) => Quadrant::NegNeg
    };
    Some (quadrant)
  }
  /// Same as `from_vec_strict` for a point
  #[inline]
  pub fn from_point_strict <S> (point : &Point2 <S>) -> Option <Quadrant> where
    S : Ring + SignedExt
  {
    Quadrant::from_vec_strict (&point.to_vector())
  }
}
//  end impl Quadrant

//
//  impl Octant
//
impl Octant {
  /// Returns some 'Some (octant)' if the vector is strictly contained within an octant
  /// and returns 'None' otherwise
  pub fn from_vec_strict <S> (vector : &Vector3 <S>) -> Option <Octant> where
    S : Ring + SignedExt
  {
    let sign_x = vector.x.sign();
    let sign_y = vector.y.sign();
    let sign_z = vector.z.sign();
    let octant = match (sign_x, sign_y, sign_z) {
      (Sign::Zero,          _,          _) |
      (         _, Sign::Zero,          _) |
      (         _,          _, Sign::Zero) => return None,
      (Sign::Positive, Sign::Positive, Sign::Positive) => Octant::PosPosPos,
      (Sign::Negative, Sign::Positive, Sign::Positive) => Octant::NegPosPos,
      (Sign::Positive, Sign::Negative, Sign::Positive) => Octant::PosNegPos,
      (Sign::Negative, Sign::Negative, Sign::Positive) => Octant::NegNegPos,
      (Sign::Positive, Sign::Positive, Sign::Negative) => Octant::PosPosNeg,
      (Sign::Negative, Sign::Positive, Sign::Negative) => Octant::NegPosNeg,
      (Sign::Positive, Sign::Negative, Sign::Negative) => Octant::PosNegNeg,
      (Sign::Negative, Sign::Negative, Sign::Negative) => Octant::NegNegNeg
    };
    Some (octant)
  }
  /// Same as `from_vec_strict` for a point
  #[inline]
  pub fn from_point_strict <S> (point : &Point3 <S>) -> Option <Octant> where
    S : Ring + SignedExt
  {
    Octant::from_vec_strict (&point.to_vector())
  }
}
//  end impl Octant

//
//  impl Positive
//
impl <S : OrderedRing> Positive <S> {
  /// Returns 'None' when called with a negative or zero value.
  ///
  /// # Example
  ///
  /// ```
  /// # use math_utils::Positive;
  /// assert!(Positive::new (1.0).is_some());
  /// assert!(Positive::new (0.0).is_none());
  /// assert!(Positive::new (-1.0).is_none());
  /// ```
  pub fn new (value : S) -> Option <Self> {
    if value > S::zero() {
      Some (Positive (value))
    } else {
      None
    }
  }
  /// Positive infinity
  pub fn infinity() -> Self where S : num::float::FloatCore {
    Positive (S::infinity())
  }
  /// Panics if negative or zero.
  ///
  /// # Panics
  ///
  /// ```should_panic
  /// # use math_utils::Positive;
  /// let x = Positive::noisy (0.0);  // panic!
  /// ```
  pub fn noisy (value : S) -> Self {
    assert!(value > S::zero());
    Positive (value)
  }
  /// Create a new positive number without checking the value.
  ///
  /// In debug builds this will fail with a debug assertion if the value is negative:
  ///
  /// ```should_panic
  /// # use math_utils::Positive;
  /// let negative = Positive::unchecked (-1.0);  // panic!
  /// ```
  // TODO: make unsafe ?
  pub fn unchecked (value : S) -> Self {
    debug_assert!(value > S::zero());
    Positive (value)
  }
  /// Map an operation on the underlying scalar, panicking if the result is zero or
  /// negative.
  ///
  /// # Example
  ///
  /// ```
  /// # use math_utils::Positive;
  /// # use math_utils::num_traits::One;
  /// assert_eq!(*Positive::<f64>::one().map_noisy (|x| 2.0 * x), 2.0);
  /// ```
  ///
  /// # Panics
  ///
  /// Panics of the result is negative or zero:
  ///
  /// ```should_panic
  /// # use math_utils::Positive;
  /// # use math_utils::num_traits::One;
  /// let v = Positive::<f64>::one().map_noisy (|_| 0.0);  // panic!
  /// ```
  pub fn map_noisy (self, fun : fn (S) -> S) -> Self {
    Self::noisy (fun (self.0))
  }
}
impl <S : Field> num::One for Positive <S> {
  fn one() -> Self {
    Positive (S::one())
  }
}
impl <S : Ring> std::ops::Deref for Positive <S> {
  type Target = S;
  fn deref (&self) -> &S {
    &self.0
  }
}
impl <S : Ring> std::ops::Mul <Positive <S>> for Positive <S> {
  type Output = Positive <S>;
  fn mul (self, rhs : Positive <S>) -> Self::Output {
    Positive (self.0 * *rhs)
  }
}
impl <S : Ring> std::ops::MulAssign <Positive <S>> for Positive <S> {
  fn mul_assign (&mut self, rhs : Positive <S>) {
    self.0 *= *rhs
  }
}
impl <S : Ring> std::ops::Mul <NonNegative <S>> for Positive <S> {
  type Output = NonNegative <S>;
  fn mul (self, rhs : NonNegative <S>) -> Self::Output {
    NonNegative (self.0 * *rhs)
  }
}
impl <S : Ring> std::ops::MulAssign <NonNegative <S>> for Positive <S> {
  fn mul_assign (&mut self, rhs : NonNegative <S>) {
    self.0 *= *rhs
  }
}
impl <S : Field> std::ops::Div for Positive <S> {
  type Output = Self;
  fn div (self, rhs : Self) -> Self::Output {
    Positive (self.0 / rhs.0)
  }
}
impl <S : Field> std::ops::DivAssign for Positive <S> {
  fn div_assign (&mut self, rhs : Self) {
    self.0 /= rhs.0
  }
}
impl <S : Ring> std::ops::Add for Positive <S> {
  type Output = Self;
  fn add (self, rhs : Self) -> Self::Output {
    Positive (self.0 + rhs.0)
  }
}
impl <S : Ring> std::ops::AddAssign for Positive <S> {
  fn add_assign (&mut self, rhs : Self) {
    self.0 += rhs.0
  }
}
impl <S : Ring> std::ops::Add <NonNegative <S>> for Positive <S> {
  type Output = Self;
  fn add (self, rhs : NonNegative <S>) -> Self::Output {
    Positive (self.0 + rhs.0)
  }
}
impl <S : Ring> std::ops::AddAssign <NonNegative <S>> for Positive <S> {
  fn add_assign (&mut self, rhs : NonNegative <S>) {
    self.0 += rhs.0
  }
}
//  end impl Positive

//
//  impl NonNegative
//
impl <S : OrderedRing> NonNegative <S> {
  /// Returns 'None' when called with a negative value.
  ///
  /// # Example
  ///
  /// ```
  /// # use math_utils::NonNegative;
  /// assert!(NonNegative::new (1.0).is_some());
  /// assert!(NonNegative::new (0.0).is_some());
  /// assert!(NonNegative::new (-1.0).is_none());
  /// ```
  pub fn new (value : S) -> Option <Self> {
    if value >= S::zero() {
      Some (NonNegative (value))
    } else {
      None
    }
  }
  /// Converts negative input to absolute value.
  ///
  /// # Example
  ///
  /// ```
  /// # use math_utils::NonNegative;
  /// assert_eq!(*NonNegative::abs (1.0), 1.0);
  /// assert_eq!(*NonNegative::abs (-1.0), 1.0);
  /// ```
  pub fn abs (value : S) -> Self {
    NonNegative (value.abs())
  }
  /// Panics if negative.
  ///
  /// # Panics
  ///
  /// ```should_panic
  /// # use math_utils::NonNegative;
  /// let x = NonNegative::noisy (-1.0);  // panic!
  /// ```
  pub fn noisy (value : S) -> Self {
    assert!(S::zero() <= value);
    NonNegative (value)
  }
  /// Create a new non-negative number without checking the value.
  ///
  /// This method is completely unchecked for release builds. Debug builds will panic if
  /// the value is negative:
  ///
  /// ```should_panic
  /// # use math_utils::NonNegative;
  /// let negative = NonNegative::unchecked (-1.0);   // panic!
  /// ```
  pub fn unchecked (value : S) -> Self {
    debug_assert!(value >= S::zero());
    NonNegative (value)
  }
  /// Map an operation on the underlying scalar, converting negative result to an
  /// absolute value.
  ///
  /// # Example
  ///
  /// ```
  /// # use math_utils::NonNegative;
  /// assert_eq!(*NonNegative::abs (1.0).map_abs (|x| -2.0 * x), 2.0);
  /// ```
  pub fn map_abs (self, fun : fn (S) -> S) -> Self {
    Self::abs (fun (self.0))
  }
  /// Map an operation on the underlying scalar, panicking if the result is negative.
  ///
  /// # Example
  ///
  /// ```
  /// # use math_utils::NonNegative;
  /// assert_eq!(*NonNegative::abs (1.0).map_noisy (|x| 2.0 * x), 2.0);
  /// ```
  ///
  /// # Panics
  ///
  /// Panics of the result is negative:
  ///
  /// ```should_panic
  /// # use math_utils::NonNegative;
  /// let v = NonNegative::abs (1.0).map_noisy (|x| -1.0 * x);  // panic!
  /// ```
  pub fn map_noisy (self, fun : fn (S) -> S) -> Self {
    Self::noisy (fun (self.0))
  }
}
impl <S : Ring> From <Positive <S>> for NonNegative <S> {
  fn from (positive : Positive <S>) -> Self {
    NonNegative (*positive)
  }
}
impl <S : Ring> num::Zero for NonNegative <S> {
  fn zero() -> Self {
    NonNegative (S::zero())
  }
  fn is_zero (&self) -> bool {
    self.0.is_zero()
  }
}
impl <S : Field> num::One for NonNegative <S> {
  fn one() -> Self {
    NonNegative (S::one())
  }
}
impl <S : Ring> std::ops::Deref for NonNegative <S> {
  type Target = S;
  fn deref (&self) -> &S {
    &self.0
  }
}
impl <S : Ring> std::ops::Mul for NonNegative <S> {
  type Output = Self;
  fn mul (self, rhs : Self) -> Self::Output {
    NonNegative (self.0 * rhs.0)
  }
}
impl <S : Ring> std::ops::Mul <Positive <S>> for NonNegative <S> {
  type Output = Self;
  fn mul (self, rhs : Positive <S>) -> Self::Output {
    NonNegative (self.0 * rhs.0)
  }
}
impl <S : Field> std::ops::Div for NonNegative <S> {
  type Output = Self;
  fn div (self, rhs : Self) -> Self::Output {
    NonNegative (self.0 / rhs.0)
  }
}
impl <S : Ring> std::ops::Add for NonNegative <S> {
  type Output = Self;
  fn add (self, rhs : Self) -> Self::Output {
    NonNegative (self.0 + rhs.0)
  }
}
//  end impl NonNegative

//
//  impl NonZero
//
impl <S : Field> NonZero <S> {
  /// Returns 'None' when called with zero.
  ///
  /// # Example
  ///
  /// ```
  /// # use math_utils::NonZero;
  /// assert!(NonZero::new (2.0).is_some());
  /// assert!(NonZero::new (0.0).is_none());
  /// assert!(NonZero::new (-2.0).is_some());
  /// ```
  #[inline]
  pub fn new (value : S) -> Option <Self> {
    if S::zero() != value {
      Some (NonZero (value))
    } else {
      None
    }
  }
  /// Panics if called with zero.
  ///
  /// # Panics
  ///
  /// ```should_panic
  /// # use math_utils::NonZero;
  /// let x = NonZero::noisy (0.0);  // panic!
  /// ```
  #[inline]
  pub fn noisy (value : S) -> Self {
    assert!(S::zero() != value);
    NonZero (value)
  }
  /// Map an operation on the underlying scalar, panicking if the result is zero.
  ///
  /// # Example
  ///
  /// ```
  /// # use math_utils::NonZero;
  /// # use math_utils::num_traits::One;
  /// assert_eq!(*NonZero::<f64>::one().map_noisy (|x| 0.5 * x), 0.5);
  /// ```
  ///
  /// # Panics
  ///
  /// Panics of the result is zero.
  ///
  /// ```should_panic
  /// # use math_utils::NonZero;
  /// # use math_utils::num_traits::One;
  /// let v = NonZero::<f64>::one().map_noisy (|x| 0.0 * x);  // panic!
  /// ```
  #[inline]
  pub fn map_noisy (self, fun : fn (S) -> S) -> Self {
    Self::noisy (fun (self.0))
  }
}
impl <S : Field> num::One for NonZero <S> {
  fn one() -> Self {
    NonZero (S::one())
  }
}
impl <S : Field> std::ops::Deref for NonZero <S> {
  type Target = S;
  fn deref (&self) -> &S {
    &self.0
  }
}
impl <S : Field> std::ops::Mul for NonZero <S> {
  type Output = Self;
  fn mul (self, rhs : Self) -> Self::Output {
    NonZero (self.0 * rhs.0)
  }
}
impl <S : Field> Eq  for NonZero <S> { }
//  end impl NonZero

//
//  impl Normalized
//
impl <S : OrderedField> Normalized <S> {
  // NOTE: num::Zero requires std::ops::Add which Normalized does not implement
  #[inline]
  pub fn zero() -> Self {
    Normalized (S::zero())
  }
  #[inline]
  pub fn is_zero (&self) -> bool {
    self.0.is_zero()
  }
  /// Returns 'None' when called with a value outside of the closed unit interval
  /// \[0,1\].
  ///
  /// # Example
  ///
  /// ```
  /// # use math_utils::Normalized;
  /// assert!(Normalized::new (2.0).is_none());
  /// assert!(Normalized::new (1.0).is_some());
  /// assert!(Normalized::new (0.5).is_some());
  /// assert!(Normalized::new (0.0).is_some());
  /// assert!(Normalized::new (-1.0).is_none());
  /// ```
  #[inline]
  pub fn new (value : S) -> Option <Self> {
    if S::zero() <= value && value <= S::one() {
      Some (Normalized (value))
    } else {
      None
    }
  }
  /// Clamps to the closed unit interval \[0,1\].
  ///
  /// # Example
  ///
  /// ```
  /// # use math_utils::Normalized;
  /// assert_eq!(*Normalized::clamp (2.0), 1.0);
  /// assert_eq!(*Normalized::clamp (-1.0), 0.0);
  /// ```
  #[inline]
  #[expect(clippy::same_name_method)]
  pub fn clamp (value : S) -> Self {
    let value = S::max (S::zero(), S::min (value, S::one()));
    Normalized (value)
  }
  /// Panics if outside the closed unit interval \[0,1\].
  ///
  /// # Panics
  ///
  /// ```should_panic
  /// # use math_utils::Normalized;
  /// let x = Normalized::noisy (-1.0);  // panic!
  /// ```
  ///
  /// ```should_panic
  /// # use math_utils::Normalized;
  /// let x = Normalized::noisy (2.0);  // panic!
  /// ```
  #[inline]
  pub fn noisy (value : S) -> Self {
    assert!(value <= S::one());
    assert!(S::zero() <= value);
    Normalized (value)
  }
  /// Map an operation on the underlying scalar, clamping results to the closed unit
  /// interval \[0,1\].
  ///
  /// # Example
  ///
  /// ```
  /// # use math_utils::Normalized;
  /// # use math_utils::num_traits::One;
  /// assert_eq!(*Normalized::<f64>::one().map_clamp (|x|  2.0 * x), 1.0);
  /// assert_eq!(*Normalized::<f64>::one().map_clamp (|x| -2.0 * x), 0.0);
  /// ```
  #[inline]
  pub fn map_clamp (self, fun : fn (S) -> S) -> Self {
    Self::clamp (fun (self.0))
  }
  /// Map an operation on the underlying scalar, panicking if the result is outside the
  /// unit interval \[0,1\].
  ///
  /// # Example
  ///
  /// ```
  /// # use math_utils::Normalized;
  /// # use math_utils::num_traits::One;
  /// assert_eq!(*Normalized::<f64>::one().map_noisy (|x| 0.5 * x), 0.5);
  /// ```
  ///
  /// # Panics
  ///
  /// Panics of the result is outside \[0,1\]:
  ///
  /// ```should_panic
  /// # use math_utils::Normalized;
  /// # use math_utils::num_traits::One;
  /// let v = Normalized::<f64>::one().map_noisy (|x| -1.0 * x);  // panic!
  /// ```
  ///
  /// ```should_panic
  /// # use math_utils::Normalized;
  /// # use math_utils::num_traits::One;
  /// let v = Normalized::<f64>::one().map_noisy (|x|  2.0 * x);  // panic!
  /// ```
  #[inline]
  pub fn map_noisy (self, fun : fn (S) -> S) -> Self {
    Self::noisy (fun (self.0))
  }
}
impl <S : Field> num::One for Normalized <S> {
  fn one() -> Self {
    Normalized (S::one())
  }
}
impl <S : Field> std::ops::Deref for Normalized <S> {
  type Target = S;
  fn deref (&self) -> &S {
    &self.0
  }
}
impl <S : Field> std::ops::Mul for Normalized <S> {
  type Output = Self;
  fn mul (self, rhs : Self) -> Self::Output {
    Normalized (self.0 * rhs.0)
  }
}
impl <S : Field> Eq  for Normalized <S> { }
#[expect(clippy::derive_ord_xor_partial_ord)]
impl <S : OrderedField> Ord for Normalized <S> {
  fn cmp (&self, other : &Self) -> std::cmp::Ordering {
    // safe to unwrap: normalized values are never NaN
    self.partial_cmp (other).unwrap()
  }
}
//  end impl Normalized

//
//  impl NormalSigned
//
impl <S : OrderedField> NormalSigned <S> {
  // NOTE: num::Zero requires std::ops::Add which Normalized does not implement
  #[inline]
  pub fn zero() -> Self {
    NormalSigned (S::zero())
  }
  #[inline]
  pub fn is_zero (&self) -> bool {
    self.0.is_zero()
  }
  /// Returns 'None' when called with a value outside of the closed interval \[-1,1\].
  ///
  /// # Example
  ///
  /// ```
  /// # use math_utils::NormalSigned;
  /// assert!(NormalSigned::new (2.0).is_none());
  /// assert!(NormalSigned::new (1.0).is_some());
  /// assert!(NormalSigned::new (0.5).is_some());
  /// assert!(NormalSigned::new (0.0).is_some());
  /// assert!(NormalSigned::new (-1.0).is_some());
  /// assert!(NormalSigned::new (-2.0).is_none());
  /// ```
  #[inline]
  pub fn new (value : S) -> Option <Self> {
    if -S::one() <= value && value <= S::one() {
      Some (NormalSigned (value))
    } else {
      None
    }
  }
  /// Clamps to the closed interval [-1,1].
  ///
  /// # Example
  ///
  /// ```
  /// # use math_utils::NormalSigned;
  /// assert_eq!(*NormalSigned::clamp (2.0), 1.0);
  /// assert_eq!(*NormalSigned::clamp (-2.0), -1.0);
  /// ```
  #[inline]
  #[expect(clippy::same_name_method)]
  pub fn clamp (value : S) -> Self {
    let value = S::max (-S::one(), S::min (value, S::one()));
    NormalSigned (value)
  }
  /// Panics if outside the closed interval [-1,1].
  ///
  /// # Panics
  ///
  /// ```should_panic
  /// # use math_utils::NormalSigned;
  /// let x = NormalSigned::noisy (-2.0);   // panic!
  /// ```
  ///
  /// ```should_panic
  /// # use math_utils::NormalSigned;
  /// let x = NormalSigned::noisy (2.0);    // panic!
  /// ```
  #[inline]
  pub fn noisy (value : S) -> Self {
    assert!(value <= S::one());
    assert!(-S::one() <= value);
    NormalSigned (value)
  }
  /// Map an operation on the underlying scalar, clamping results to the closed interval
  /// [-1,1].
  ///
  /// # Example
  ///
  /// ```
  /// # use math_utils::NormalSigned;
  /// # use math_utils::num_traits::One;
  /// assert_eq!(*NormalSigned::<f64>::one().map_clamp (|x|  2.0 * x),  1.0);
  /// assert_eq!(*NormalSigned::<f64>::one().map_clamp (|x| -2.0 * x), -1.0);
  /// ```
  #[inline]
  pub fn map_clamp (self, fun : fn (S) -> S) -> Self {
    Self::clamp (fun (self.0))
  }
  /// Map an operation on the underlying scalar, panicking if the result is outside the
  /// interval [-1, 1].
  ///
  /// # Example
  ///
  /// ```
  /// # use math_utils::NormalSigned;
  /// # use math_utils::num_traits::One;
  /// assert_eq!(*NormalSigned::<f64>::one().map_noisy (|x| 0.5 * x), 0.5);
  /// ```
  ///
  /// # Panics
  ///
  /// Panics of the result is outside [-1, 1]:
  ///
  /// ```should_panic
  /// # use math_utils::NormalSigned;
  /// # use math_utils::num_traits::One;
  /// let v = NormalSigned::<f64>::one().map_noisy (|x| -2.0 * x);  // panic!
  /// ```
  ///
  /// ```should_panic
  /// # use math_utils::NormalSigned;
  /// # use math_utils::num_traits::One;
  /// let v = NormalSigned::<f64>::one().map_noisy (|x|  2.0 * x);  // panic!
  /// ```
  #[inline]
  pub fn map_noisy (self, fun : fn (S) -> S) -> Self {
    Self::noisy (fun (self.0))
  }
}
impl <S : OrderedField> From <Normalized <S>> for NormalSigned <S> {
  fn from (Normalized (value) : Normalized <S>) -> Self {
    debug_assert!(S::zero() <= value);
    debug_assert!(value <= S::one());
    NormalSigned (value)
  }
}
impl <S : Field> num::One for NormalSigned <S> {
  fn one() -> Self {
    NormalSigned (S::one())
  }
}
impl <S : Field> std::ops::Deref for NormalSigned <S> {
  type Target = S;
  fn deref (&self) -> &S {
    &self.0
  }
}
impl <S : Field> std::ops::Mul for NormalSigned <S> {
  type Output = Self;
  fn mul (self, rhs : Self) -> Self::Output {
    NormalSigned (self.0 * rhs.0)
  }
}
impl <S : Field> Eq  for NormalSigned <S> { }
#[expect(clippy::derive_ord_xor_partial_ord)]
impl <S : OrderedField> Ord for NormalSigned <S> {
  fn cmp (&self, other : &Self) -> std::cmp::Ordering {
    // safe to unwrap: normalized values are never NaN
    self.partial_cmp (other).unwrap()
  }
}
impl <S : Field> std::ops::Neg for NormalSigned <S> {
  type Output = Self;
  fn neg (self) -> Self::Output {
    NormalSigned (-self.0)
  }
}
//  end impl NormalSigned

//
//  impl Rotation2
//
impl <S> Rotation2 <S> where S : Ring + MaybeSerDes {
  /// Identity rotation
  pub fn identity() -> Self {
    use num::One;
    Self::one()
  }
  /// Returns `None` if called with a matrix that is not orthonormal with determinant
  /// +1.
  ///
  /// This method checks whether `mat * mat^T == I` and `mat.determinant() == 1`.
  pub fn new (mat : Matrix2 <S>) -> Option <Self> {
    let transform = LinearAuto (LinearIso::new (mat)?);
    if !transform.is_rotation() {
      None
    } else {
      Some (Rotation2 (transform))
    }
  }
  /// Returns `None` if called with a matrix that is not orthonormal with determinant
  /// +1.
  ///
  /// This method checks whether `mat * mat^T == I` and `mat.determinant() == 1`
  /// (approximately using relative equality).
  pub fn new_approx (mat : Matrix2 <S>) -> Option <Self> where
    S : approx::RelativeEq <Epsilon=S>
  {
    let four         = S::one() + S::one() + S::one() + S::one();
    let thirtytwo    = four * four * (S::one() + S::one());
    let epsilon      = S::default_epsilon() * thirtytwo;
    let max_relative = S::default_max_relative() * four;
    if approx::relative_ne!(mat * mat.transposed(), Matrix2::identity(),
      max_relative=max_relative, epsilon=epsilon
    ) || approx::relative_ne!(mat.determinant(), S::one(), max_relative=max_relative) {
      None
    } else {
      Some (Rotation2 (LinearAuto (LinearIso (mat, PhantomData))))
    }
  }
  /// Create a rotation for the given angle
  pub fn from_angle (angle : Rad <S>) -> Self where S : num::real::Real {
    Rotation2 (LinearAuto (LinearIso (Matrix2::rotation_z (angle.0), PhantomData)))
  }
  /// Panic if the given matrix is not orthogonal with determinant +1.
  ///
  /// This method checks whether `mat * mat^T == I` and `mat.determinant() == 1`.
  pub fn noisy (mat : Matrix2 <S>) -> Self where S : std::fmt::Debug {
    assert_eq!(mat * mat.transposed(), Matrix2::identity());
    assert_eq!(mat.determinant(), S::one());
    Rotation2 (LinearAuto (LinearIso (mat, PhantomData)))
  }
  /// Panic if the given matrix is not orthogonal with determinant +1.
  ///
  /// This method checks whether `mat * mat^T == I` and `mat.determinant() == 1`
  /// (approximately using relative equality).
  pub fn noisy_approx (mat : Matrix2 <S>) -> Self where
    S : approx::RelativeEq <Epsilon=S> + std::fmt::Debug
  {
    let four         = S::one() + S::one() + S::one() + S::one();
    let epsilon      = S::default_epsilon() * four;
    let max_relative = S::default_max_relative() * four;
    approx::assert_relative_eq!(mat * mat.transposed(), Matrix2::identity(),
      epsilon=epsilon, max_relative=max_relative);
    approx::assert_relative_eq!(mat.determinant(), S::one(), max_relative=max_relative);
    Rotation2 (LinearAuto (LinearIso (mat, PhantomData)))
  }
  /// It is a debug panic if the given matrix is not orthogonal with determinant +1.
  ///
  /// This method checks whether `mat * mat^T == I` and `mat.determinant() == 1`.
  pub fn unchecked (mat : Matrix2 <S>) -> Self where S : std::fmt::Debug {
    debug_assert!(mat * mat.transposed() == Matrix2::identity());
    debug_assert!(mat.determinant() == S::one());
    Rotation2 (LinearAuto (LinearIso (mat, PhantomData)))
  }
  /// It is a debug panic if the given matrix is not orthogonal with determinant +1.
  ///
  /// This method checks whether `mat * mat^T == I` and `mat.determinant() == 1`
  /// (approximately using relative equality).
  pub fn unchecked_approx (mat : Matrix2 <S>) -> Self where
    S : approx::RelativeEq <Epsilon=S> + std::fmt::Debug
  {
    if cfg!(debug_assertions) {
      let four         = S::one() + S::one() + S::one() + S::one();
      let epsilon      = S::default_epsilon() * four;
      let max_relative = S::default_max_relative() * four;
      approx::assert_relative_eq!(mat * mat.transposed(), Matrix2::identity(),
        max_relative=max_relative, epsilon=epsilon);
      approx::assert_relative_eq!(mat.determinant(), S::one(),
        max_relative=max_relative);
    }
    Rotation2 (LinearAuto (LinearIso (mat, PhantomData)))
  }
  /// Get the underlying matrix
  pub const fn mat (self) -> Matrix2 <S> {
    self.0.0.0
  }
  /// Rotates a point around the origin
  pub fn rotate (self, point : Point2 <S>) -> Point2 <S> {
    Point2::from (self * point.0)
  }
  /// Rotate around a center point
  pub fn rotate_around (self, point : Point2 <S>, center : Point2 <S>) -> Point2 <S> {
    self.rotate (point - center.0) + center.0
  }
}
impl <S> std::ops::Deref for Rotation2 <S> where S : Ring + MaybeSerDes {
  type Target = LinearAuto <S, Vector2 <S>>;
  fn deref (&self) -> &Self::Target {
    &self.0
  }
}
impl <S> num::Inv for Rotation2 <S> where S : Ring + MaybeSerDes {
  type Output = Self;
  fn inv (self) -> Self::Output {
    Rotation2 (LinearAuto (LinearIso (self.0.0.0.transpose(), PhantomData)))
  }
}
impl <S> std::ops::Mul <Vector2 <S>> for Rotation2 <S> where S : Ring + MaybeSerDes {
  type Output = Vector2 <S>;
  fn mul (self, rhs : Vector2 <S>) -> Self::Output {
    self.0 * rhs
  }
}
impl <S> std::ops::Mul <Rotation2 <S>> for Vector2 <S> where S : Ring + MaybeSerDes {
  type Output = Vector2 <S>;
  fn mul (self, rhs : Rotation2 <S>) -> Self::Output {
    self * rhs.0
  }
}
impl <S> std::ops::Mul for Rotation2 <S> where S : Ring + MaybeSerDes {
  type Output = Self;
  fn mul (self, rhs : Self) -> Self::Output {
    Rotation2 (self.0 * rhs.0)
  }
}
impl <S> std::ops::MulAssign <Self> for Rotation2 <S> where S : Ring + MaybeSerDes {
  fn mul_assign (&mut self, rhs : Self) {
    self.0 *= rhs.0
  }
}
impl <S> std::ops::Div for Rotation2 <S> where S : Ring + MaybeSerDes {
  type Output = Self;
  #[expect(clippy::suspicious_arithmetic_impl)]
  fn div (self, rhs : Self) -> Self::Output {
    use num::Inv;
    self * rhs.inv()
  }
}
impl <S> std::ops::DivAssign <Self> for Rotation2 <S> where S : Ring + MaybeSerDes {
  #[expect(clippy::suspicious_op_assign_impl)]
  fn div_assign (&mut self, rhs : Self) {
    use num::Inv;
    *self *= rhs.inv()
  }
}
impl <S> num::One for Rotation2 <S> where S : Ring + MaybeSerDes + num::Zero + num::One {
  fn one() -> Self {
    Rotation2 (LinearAuto (LinearIso (Matrix2::identity(), PhantomData)))
  }
}

//
//  impl Rotation3
//
impl <S> Rotation3 <S> where S : Ring + MaybeSerDes {
  /// Identity rotation
  pub fn identity() -> Self {
    use num::One;
    Self::one()
  }
  /// Returns `None` if called with a matrix that is not orthonormal with determinant
  /// +1.
  ///
  /// This method checks whether `mat * mat^T == I` and `mat.determinant() == 1`.
  pub fn new (mat : Matrix3 <S>) -> Option <Self> {
    if mat * mat.transposed() != Matrix3::identity() || mat.determinant() != S::one() {
      None
    } else {
      Some (Rotation3 (LinearAuto (LinearIso (mat, PhantomData))))
    }
  }
  /// Returns `None` if called with a matrix that is not orthonormal with determinant
  /// +1.
  ///
  /// This method checks whether `mat * mat^T == I` and `mat.determinant() == 1`
  /// (approximately using relative equality).
  pub fn new_approx (mat : Matrix3 <S>) -> Option <Self> where
    S : approx::RelativeEq <Epsilon=S>
  {
    let four         = S::one() + S::one() + S::one() + S::one();
    let thirtytwo    = four * four * (S::one() + S::one());
    let epsilon      = S::default_epsilon() * thirtytwo;
    let max_relative = S::default_max_relative() * four;
    if approx::relative_ne!(mat * mat.transposed(), Matrix3::identity(),
      max_relative=max_relative, epsilon=epsilon
    ) || approx::relative_ne!(mat.determinant(), S::one(), max_relative=max_relative) {
      None
    } else {
      Some (Rotation3 (LinearAuto (LinearIso (mat, PhantomData))))
    }
  }
  /// Panic if the given matrix is not orthogonal with determinant +1.
  ///
  /// This method checks whether `mat * mat^T == I` and `mat.determinant() == 1`.
  pub fn noisy (mat : Matrix3 <S>) -> Self {
    assert!(mat * mat.transposed() == Matrix3::identity());
    assert!(mat.determinant() == S::one());
    Rotation3 (LinearAuto (LinearIso (mat, PhantomData)))
  }
  /// Panic if the given matrix is not orthogonal with determinant +1.
  ///
  /// This method checks whether `mat * mat^T == I` and `mat.determinant() == 1`
  /// (approximately using relative equality).
  pub fn noisy_approx (mat : Matrix3 <S>) -> Self where
    S : approx::RelativeEq <Epsilon=S> + std::fmt::Debug
  {
    let four         = S::one() + S::one() + S::one() + S::one();
    let eight        = four + S::one() + S::one() + S::one() + S::one();
    let epsilon      = S::default_epsilon() * four;
    let max_relative = S::default_max_relative() * eight;
    approx::assert_relative_eq!(mat * mat.transposed(), Matrix3::identity(),
      max_relative=max_relative, epsilon=epsilon);
    approx::assert_relative_eq!(mat.determinant(), S::one(), max_relative=max_relative);
    Rotation3 (LinearAuto (LinearIso (mat, PhantomData)))
  }
  /// It is a debug panic if the given matrix is not orthogonal with determinant +1.
  ///
  /// This method checks whether `mat * mat^T == I` and `mat.determinant() == 1`.
  pub fn unchecked (mat : Matrix3 <S>) -> Self where S : std::fmt::Debug {
    debug_assert!(mat * mat.transposed() == Matrix3::identity());
    debug_assert!(mat.determinant() == S::one());
    Rotation3 (LinearAuto (LinearIso (mat, PhantomData)))
  }
  /// It is a debug panic if the given matrix is not orthogonal with determinant +1.
  ///
  /// This method checks whether `mat * mat^T == I` and `mat.determinant() == 1`
  /// (approximately using relative equality).
  pub fn unchecked_approx (mat : Matrix3 <S>) -> Self where
    S : approx::RelativeEq <Epsilon=S> + std::fmt::Debug
  {
    if cfg!(debug_assertions) {
      let four         = S::one() + S::one() + S::one() + S::one();
      let epsilon      = S::default_epsilon() * four;
      let max_relative = S::default_max_relative() * four;
      approx::assert_relative_eq!(mat * mat.transposed(), Matrix3::identity(),
        max_relative=max_relative, epsilon=epsilon);
      approx::assert_relative_eq!(mat.determinant(), S::one(),
        max_relative=max_relative);
    }
    Rotation3 (LinearAuto (LinearIso (mat, PhantomData)))
  }
  /// Construct a rotation around the X axis.
  ///
  /// Positive angles are counter-clockwise.
  pub fn from_angle_x (angle : Rad <S>) -> Self where S : num::real::Real {
    Rotation3 (LinearAuto (LinearIso (Matrix3::rotation_x (angle.0), PhantomData)))
  }
  /// Construct a rotation around the Y axis
  ///
  /// Positive angles are counter-clockwise.
  pub fn from_angle_y (angle : Rad <S>) -> Self where S : num::real::Real {
    Rotation3 (LinearAuto (LinearIso (Matrix3::rotation_y (angle.0), PhantomData)))
  }
  /// Construct a rotation around the Z axis
  ///
  /// Positive angles are counter-clockwise.
  pub fn from_angle_z (angle : Rad <S>) -> Self where S : num::real::Real {
    Rotation3 (LinearAuto (LinearIso (Matrix3::rotation_z (angle.0), PhantomData)))
  }
  /// Construct a rotation from intrinsic ZX'Y'' Euler angles:
  ///
  /// - yaw: rotation around Z axis
  /// - pitch: rotation around X' axis
  /// - roll: rotation around Y'' axis
  pub fn from_angles_intrinsic (yaw : Rad <S>, pitch : Rad <S>, roll : Rad <S>)
    -> Self where S : num::real::Real
  {
    let rotation1 = Rotation3::from_angle_z (yaw);
    let rotation2 = Rotation3::from_axis_angle (rotation1.cols.x, pitch) * rotation1;
    Rotation3::from_axis_angle (rotation2.cols.y, roll) * rotation2
  }
  /// Construct a rotation around the given axis vector with the given angle
  pub fn from_axis_angle (axis : Vector3 <S>, angle : Rad <S>) -> Self where
    S : num::real::Real
  {
    Rotation3 (LinearAuto (LinearIso (
      Matrix3::rotation_3d (angle.0, axis), PhantomData)))
  }
  /// Construct an orthonormal matrix from a set of linearly-independent vectors using
  /// the Gram-Schmidt process
  pub fn orthonormalize (v1 : Vector3 <S>, v2 : Vector3 <S>, v3 : Vector3 <S>)
    -> Option <Self> where S : Field + Sqrt
  {
    if Matrix3::from_col_arrays ([v1, v2, v3].map (Vector3::into_array)).determinant()
      == S::zero()
    {
      None
    } else {
      let project = |u : Vector3 <S>, v : Vector3 <S>| u * (u.dot (v) / u.self_dot());
      let u1 = v1;
      let u2 = v2 - project (u1, v2);
      let u3 = v3 - project (u1, v3) - project (u2, v3);
      Some (Rotation3 (LinearAuto (LinearIso (Matrix3::from_col_arrays ([
        u1.normalize().into_array(),
        u2.normalize().into_array(),
        u3.normalize().into_array()
      ]), PhantomData))))
    }
  }

  /// Returns a rotation with zero roll and oriented towards the target point.
  ///
  /// Returns the identity rotation if `point` is the origin.
  pub fn look_at (point : Point3 <S>) -> Self where
    S : num::Float + num::FloatConst + approx::RelativeEq <Epsilon=S> + std::fmt::Debug
  {
    if point == Point::origin() {
      Self::identity()
    } else {
      use num::Zero;
      let forward = point.0.normalized();
      let right = {
        let projected = forward.with_z (S::zero());
        if !projected.is_zero() {
          Self::from_angle_z (-Rad (S::FRAC_PI_2())) * projected.normalized()
        } else {
          Vector3::unit_x()
        }
      };
      let up = right.cross (forward);
      let mat =
        Matrix3::from_col_arrays ([right, forward, up].map (Vector3::into_array));
      Self::noisy_approx (mat)
    }
  }

  /// Return intrinsic yaw (Z), pitch (X'), roll (Y'') angles.
  ///
  /// NOTE: this function returns the raw output of the `atan2` operations used to
  /// compute the angles. This function returns values in the range `[-\pi, \pi]`.
  pub fn intrinsic_angles (&self) -> (Rad <S>, Rad <S>, Rad <S>) where S : Real {
    let cols  = &self.cols;
    let yaw   = Rad (S::atan2 (-cols.y.x, cols.y.y));
    let pitch = Rad (S::atan2 (cols.y.z, S::sqrt (S::one() - cols.y.z * cols.y.z)));
    let roll  = Rad (S::atan2 (-cols.x.z, cols.z.z));
    (yaw, pitch, roll)
  }

  /// Convert to versor (unit quaternion) representation.
  ///
  /// <https://stackoverflow.com/a/63745757>
  pub fn versor (self) -> Versor <S> where
    S : Real + MaybeSerDes + num::real::Real + approx::RelativeEq <Epsilon=S>
  {
    let diagonal = self.diagonal();
    let t     = diagonal.sum();
    let m     = MinMax::max (MinMax::max (MinMax::max (
      diagonal.x, diagonal.y), diagonal.z), t);
    let qmax  = S::half() * Sqrt::sqrt (S::one() - t + S::two() * m);
    let qmax4_recip = (S::four() * qmax).recip();
    let [qx, qy, qz, qw] : [S; 4];
    let cols = self.cols;
    if m == diagonal.x {
      qx = qmax;
      qy = qmax4_recip * (cols.x.y + cols.y.x);
      qz = qmax4_recip * (cols.x.z + cols.z.x);
      qw = -qmax4_recip * (cols.z.y - cols.y.z);
    } else if m == diagonal.y {
      qx = qmax4_recip * (cols.x.y + cols.y.x);
      qy = qmax;
      qz = qmax4_recip * (cols.y.z + cols.z.y);
      qw = -qmax4_recip * (cols.x.z - cols.z.x);
    } else if m == diagonal.z {
      qx = qmax4_recip * (cols.x.z + cols.z.x);
      qy = qmax4_recip * (cols.y.z + cols.z.y);
      qz = qmax;
      qw = -qmax4_recip * (cols.y.x - cols.x.y);
    } else {
      qx = -qmax4_recip * (cols.z.y - cols.y.z);
      qy = -qmax4_recip * (cols.x.z - cols.z.x);
      qz = -qmax4_recip * (cols.y.x - cols.x.y);
      qw = qmax;
    }
    let quat = Quaternion::from_xyzw (qx, qy, qz, qw);
    Versor::unchecked_approx (quat)
  }
  /// Get the underlying matrix
  pub const fn mat (self) -> Matrix3 <S> {
    self.0.0.0
  }
  /// Rotates a point around the origin
  pub fn rotate (self, point : Point3 <S>) -> Point3 <S> {
    Point3::from (self * point.0)
  }
  /// Rotate around a center point
  pub fn rotate_around (self, point : Point3 <S>, center : Point3 <S>) -> Point3 <S> {
    self.rotate (point - center.0) + center.0
  }
}
impl <S> std::ops::Deref for Rotation3 <S> where S : Ring + MaybeSerDes {
  type Target = LinearAuto <S, Vector3 <S>>;
  fn deref (&self) -> &Self::Target {
    &self.0
  }
}
impl <S> num::Inv for Rotation3 <S> where S : Ring + MaybeSerDes {
  type Output = Self;
  fn inv (self) -> Self::Output {
    Rotation3 (LinearAuto (LinearIso (self.0.0.0.transpose(), PhantomData)))
  }
}
impl <S> std::ops::Mul <Vector3 <S>> for Rotation3 <S> where S : Ring + MaybeSerDes {
  type Output = Vector3 <S>;
  fn mul (self, rhs : Vector3 <S>) -> Self::Output {
    self.0 * rhs
  }
}
impl <S> std::ops::Mul <Rotation3 <S>> for Vector3 <S> where S : Ring + MaybeSerDes {
  type Output = Vector3 <S>;
  fn mul (self, rhs : Rotation3 <S>) -> Self::Output {
    self * rhs.0
  }
}
impl <S> std::ops::Mul for Rotation3 <S> where S : Ring + MaybeSerDes {
  type Output = Self;
  fn mul (self, rhs : Self) -> Self::Output {
    Rotation3 (self.0 * rhs.0)
  }
}
impl <S> std::ops::MulAssign <Self> for Rotation3 <S> where S : Ring + MaybeSerDes {
  fn mul_assign (&mut self, rhs : Self) {
    self.0 *= rhs.0
  }
}
impl <S> std::ops::Div for Rotation3 <S> where S : Ring + MaybeSerDes {
  type Output = Self;
  #[expect(clippy::suspicious_arithmetic_impl)]
  fn div (self, rhs : Self) -> Self::Output {
    use num::Inv;
    self * rhs.inv()
  }
}
impl <S> std::ops::DivAssign <Self> for Rotation3 <S> where S : Ring + MaybeSerDes {
  #[expect(clippy::suspicious_op_assign_impl)]
  fn div_assign (&mut self, rhs : Self) {
    use num::Inv;
    *self *= rhs.inv()
  }
}
impl <S> num::One for Rotation3 <S> where S : Ring + MaybeSerDes + num::Zero + num::One {
  fn one() -> Self {
    Rotation3 (LinearAuto (LinearIso (Matrix3::identity(), PhantomData)))
  }
}
impl <S> From <Angles3 <S>> for Rotation3 <S> where
  S : Real + num::real::Real + MaybeSerDes
{
  fn from (angles : Angles3 <S>) -> Self {
    Rotation3::from_angles_intrinsic (
      angles.yaw.angle(), angles.pitch.angle(), angles.roll.angle())
  }
}

//
//  impl Versor
//
impl <S : Real> Versor <S> {
  /// Normalizes the given quaternion.
  ///
  /// Panics if the zero quaternion is given.
  pub fn normalize (quaternion : Quaternion <S>) -> Self where S : std::fmt::Debug {
    assert_ne!(quaternion, Quaternion::zero());
    Versor (normalize_quaternion (quaternion))
  }
  /// Panic if the given quaternion is not a unit quaternion.
  ///
  /// This method checks whether `quaternion == quaternion.normalized()`.
  pub fn noisy (unit_quaternion : Quaternion <S>) -> Self where
    S : num::real::Real + std::fmt::Debug
  {
    assert_eq!(unit_quaternion, normalize_quaternion (unit_quaternion));
    Versor (unit_quaternion)
  }
  /// Panic if the given quaternion is not a unit quaternion.
  ///
  /// Checks if the magnitude is approximately one.
  pub fn noisy_approx (unit_quaternion : Quaternion <S>) -> Self where
    S : num::real::Real + approx::RelativeEq <Epsilon=S>
  {
    assert!(unit_quaternion.into_vec4().is_normalized());
    Versor (unit_quaternion)
  }
  /// It is a debug panic if the given quaternion is not a unit quaternion.
  ///
  /// This method checks whether `quaternion == quaternion.normalized()`.
  pub fn unchecked (unit_quaternion : Quaternion <S>) -> Self where S : std::fmt::Debug {
    debug_assert_eq!(unit_quaternion, normalize_quaternion (unit_quaternion));
    Versor (unit_quaternion)
  }
  /// It is a debug panic if the given quaternion is not a unit quaternion.
  ///
  /// Checks if the magnitude is approximately one.
  pub fn unchecked_approx (unit_quaternion : Quaternion <S>) -> Self where
    S : num::real::Real + approx::RelativeEq <Epsilon=S>
  {
    debug_assert!(unit_quaternion.into_vec4().is_normalized());
    Versor (unit_quaternion)
  }
  /// Constructs the rotation using the direction and magnitude of the given vector
  // TODO: work around rotation_3d Float constraint
  pub fn from_scaled_axis (scaled_axis : Vector3 <S>) -> Self where
    S : std::fmt::Debug + num::Float
  {
    if scaled_axis == Vector3::zero() {
      let versor = Quaternion::rotation_3d (S::zero(), scaled_axis);
      debug_assert_eq!(versor.magnitude(), S::one());
      Versor (versor)
    } else {
      let angle = scaled_axis.magnitude();
      debug_assert!(S::zero() < angle);
      let axis  = scaled_axis / angle;
      Versor (Quaternion::rotation_3d (angle, axis))
    }
  }
}
impl <S> std::ops::Mul <Vector3 <S>> for Versor <S> where S : Real + num::Float {
  type Output = Vector3 <S>;
  fn mul (self, rhs : Vector3 <S>) -> Self::Output {
    self.0 * rhs
  }
}
impl <S> From <Rotation3 <S>> for Versor <S> where
  S : Real + MaybeSerDes + num::real::Real + approx::RelativeEq <Epsilon=S>
{
  fn from (rot : Rotation3 <S>) -> Self {
    rot.versor()
  }
}
impl <S : Real> std::ops::Deref for Versor <S> {
  type Target = Quaternion <S>;
  fn deref (&self) -> &Quaternion <S> {
    &self.0
  }
}
//  end impl Versor

//
//  impl LinearIso
//
impl <S, V, W, M> LinearIso <S, V, W, M> where
  M : LinearMap <S, V, W> + Copy,
  V : Module <S>,
  W : Module <S>,
  S : Ring
{
  pub fn mat (self) -> M {
    *self
  }
  /// Checks whether the determinant is non-zero
  pub fn is_invertible (linear_map : M) -> bool {
    linear_map.determinant() != S::zero()
  }
  /// Checks whether the determinant is non-zero with approximate equality
  pub fn is_invertible_approx (linear_map : M, epsilon : Option <S>) -> bool where
    S : approx::AbsDiffEq <Epsilon=S>
  {
    let epsilon = epsilon.unwrap_or_else (||{
      let four = S::one() + S::one() + S::one() + S::one();
      S::default_epsilon() * four
    });
    approx::abs_diff_ne!(linear_map.determinant(), S::zero(), epsilon=epsilon)
  }
  /// Checks if map is invertible
  pub fn new (linear_map : M) -> Option <Self> {
    if Self::is_invertible (linear_map) {
      Some (LinearIso (linear_map, PhantomData))
    } else {
      None
    }
  }
  /// Checks if map is invertible with approximate equality
  pub fn new_approx (linear_map : M, epsilon : Option <S>) -> Option <Self> where
    S : approx::RelativeEq <Epsilon=S>
  {
    if Self::is_invertible_approx (linear_map, epsilon) {
      Some (LinearIso (linear_map, PhantomData))
    } else {
      None
    }
  }
}
impl <S, V, W, M> std::ops::Deref for LinearIso <S, V, W, M> where
  M : LinearMap <S, V, W>,
  V : Module <S>,
  W : Module <S>,
  S : Ring
{
  type Target = M;
  fn deref (&self) -> &Self::Target {
    &self.0
  }
}
impl <S, V> num::One for LinearIso <S, V, V, V::LinearEndo> where
  V : Module <S>,
  S : Ring
{
  fn one () -> Self {
    LinearIso (V::LinearEndo::one(), PhantomData)
  }
}
impl <S, V, W, M> std::ops::Mul <V> for LinearIso <S, V, W, M> where
  M : LinearMap <S, V, W>,
  V : Module <S>,
  W : Module <S>,
  S : Ring
{
  type Output = W;
  fn mul (self, rhs : V) -> Self::Output {
    self.0 * rhs
  }
}
impl <S, V> std::ops::Mul for LinearIso <S, V, V, V::LinearEndo> where
  V : Module <S>,
  S : Ring
{
  type Output = Self;
  fn mul (self, rhs : Self) -> Self::Output {
    LinearIso (self.0 * rhs.0, PhantomData)
  }
}
impl <S, V> std::ops::MulAssign for LinearIso <S, V, V, V::LinearEndo> where
  V : Module <S>,
  S : Ring
{
  fn mul_assign (&mut self, rhs : Self) {
    self.0 *= rhs.0
  }
}

//
//  impl LinearAuto
//
impl <S, V> LinearAuto <S, V> where
  V : Module <S>,
  V::LinearEndo : MaybeSerDes,
  S : Ring
{
  pub fn mat (self) -> V::LinearEndo {
    **self
  }
  /// Returns false if called with a matrix that is not orthogonal with determinant
  /// +/-1.
  ///
  /// This method checks whether `mat * mat^T == I` and `mat.determinant() == +/-1`.
  pub fn is_orthonormal (&self) -> bool {
    use num::One;
    let determinant = self.0.0.determinant();
    self.0.0 * self.0.0.transpose() == V::LinearEndo::one() &&
    (determinant == S::one() || determinant == -S::one())
  }
  /// Returns `None` if called with a matrix that is not orthogonal with determinant
  /// +/-1.
  ///
  /// This method checks whether `mat * mat^T == I` and `mat.determinant() == +/-1`
  /// (approximately using relative equality).
  pub fn is_orthonormal_approx (&self) -> bool where
    S : approx::RelativeEq <Epsilon=S>,
    V::LinearEndo : approx::RelativeEq <Epsilon=S>
  {
    use num::One;
    let four         = S::one() + S::one() + S::one() + S::one();
    let epsilon      = S::default_epsilon() * four;
    let max_relative = S::default_max_relative() * four;
    let determinant  = self.0.0.determinant();
    approx::relative_eq!(self.0.0 * self.0.0.transpose(), V::LinearEndo::one(),
      max_relative=max_relative, epsilon=epsilon) &&
    ( approx::relative_eq!(determinant,  S::one(), max_relative=max_relative) ||
      approx::relative_eq!(determinant, -S::one(), max_relative=max_relative) )
  }
  /// Checks whether the transformation is a special orthogonal matrix.
  ///
  /// This method checks whether `mat * mat^T == I` and `mat.determinant() == 1`.
  pub fn is_rotation (&self) -> bool {
    use num::One;
    let determinant = self.0.0.determinant();
    self.0.0 * self.0.0.transpose() == V::LinearEndo::one() &&
    determinant == S::one()
  }
  /// Checks whether the transformation is a special orthogonal matrix.
  ///
  /// This method checks whether `mat * mat^T == I` and `mat.determinant() == 1`
  /// (approximately using relative equality).
  pub fn is_rotation_approx (&self) -> bool where
    S : approx::RelativeEq <Epsilon=S>,
    V::LinearEndo : approx::RelativeEq <Epsilon=S>
  {
    use num::One;
    let four         = S::one() + S::one() + S::one() + S::one();
    let epsilon      = S::default_epsilon() * four;
    let max_relative = S::default_max_relative() * four;
    let determinant  = self.0.0.determinant();
    approx::relative_eq!(self.0.0 * self.0.0.transpose(), V::LinearEndo::one(),
      max_relative=max_relative, epsilon=epsilon) &&
    approx::relative_eq!(determinant,  S::one(), max_relative=max_relative)
  }
}
impl <S, V> std::ops::Mul <V> for LinearAuto <S, V> where
  V : Module <S>,
  V::LinearEndo : MaybeSerDes,
  S : Ring
{
  type Output = V;
  fn mul (self, rhs : V) -> Self::Output {
    self.0 * rhs
  }
}
impl <S, V> std::ops::Mul for LinearAuto <S, V> where
  V : Module <S>,
  V::LinearEndo : MaybeSerDes,
  S : Ring
{
  type Output = Self;
  fn mul (self, rhs : Self) -> Self::Output {
    LinearAuto (self.0 * rhs.0)
  }
}
impl <S, V> std::ops::MulAssign <Self> for LinearAuto <S, V> where
  V : Module <S>,
  V::LinearEndo : MaybeSerDes,
  S : Ring
{
  fn mul_assign (&mut self, rhs : Self) {
    self.0 *= rhs.0
  }
}
impl <S, V> num::One for LinearAuto <S, V> where
  V : Module <S>,
  V::LinearEndo : MaybeSerDes,
  S : Ring
{
  fn one() -> Self {
    LinearAuto (LinearIso::one())
  }
}

//
//  impl AffineMap
//
impl <S, A, B, M> AffineMap <S, A, B, M> where
  A : AffineSpace <S>,
  B : AffineSpace <S>,
  M : LinearMap <S, A::Translation, B::Translation>,
  S : Field
{
  pub const fn new (linear_map : M, translation : B::Translation) -> Self {
    AffineMap { linear_map, translation, _phantom: PhantomData }
  }
  pub fn map (self, point : A) -> B {
    B::from_vector (self.linear_map * point.to_vector() + self.translation)
  }
  pub fn transform (self, translation : A::Translation) -> B::Translation {
    self.linear_map * translation
  }
}

//
//  impl Affinity
//
impl <S, A, B, M> Affinity <S, A, B, M> where
  M : LinearMap <S, A::Translation, B::Translation>,
  A : AffineSpace <S>,
  B : AffineSpace <S>,
  S : Field
{
  pub const fn new (
    linear_iso  : LinearIso <S, A::Translation, B::Translation, M>,
    translation : B::Translation
  ) -> Self {
    Affinity { linear_iso, translation }
  }
  pub fn mat (self) -> M {
    *self.linear_iso
  }
  pub fn map (self, point : A) -> B {
    B::from_vector (self.linear_iso * point.to_vector() + self.translation)
  }
  pub fn transform (self, translation : A::Translation) -> B::Translation {
    self.linear_iso * translation
  }
}

//
//  impl Projectivity
//
impl <S, P, Q, M> Projectivity <S, P, Q, M> where
  M : LinearMap <S, P::Translation, Q::Translation>,
  P : ProjectiveSpace <S>,
  Q : ProjectiveSpace <S>,
  S : Field
{
  pub const fn new (linear_iso : LinearIso <S, P::Translation, Q::Translation, M>)
    -> Self
  {
    Projectivity (linear_iso)
  }
  pub fn mat (self) -> M {
    **self
  }
}
impl <S, P, Q, M> std::ops::Deref for Projectivity <S, P, Q, M> where
  M : LinearMap <S, P::Translation, Q::Translation>,
  P : ProjectiveSpace <S>,
  Q : ProjectiveSpace <S>,
  S : Field
{
  type Target = LinearIso <S, P::Translation, Q::Translation, M>;
  fn deref (&self) -> &Self::Target {
    &self.0
  }
}
impl <S, P, Q, M> std::ops::Mul <P> for Projectivity <S, P, Q, M> where
  M : LinearMap <S, P::Translation, Q::Translation>,
  P : ProjectiveSpace <S>,
  Q : ProjectiveSpace <S>,
  S : Field
{
  type Output = Q;
  fn mul (self, rhs : P) -> Q {
    Q::from_vector (self.0 * rhs.to_vector())
  }
}

pub macro vector {
  ($x:expr, $y:expr) => {
    $crate::vector2 ($x, $y)
  },
  ($x:expr, $y:expr, $z:expr) => {
    $crate::vector3 ($x, $y, $z)
  },
  ($x:expr, $y:expr, $z:expr, $w:expr) => {
    $crate::vector4 ($x, $y, $z, $w)
  }
}

pub macro matrix {
  ($v00:expr, $v01:expr; $v10:expr, $v11:expr$(;)?) => {
    $crate::Matrix2::new ($v00, $v01, $v10, $v11)
  },
  ( $v00:expr, $v01:expr, $v02:expr;
    $v10:expr, $v11:expr, $v12:expr;
    $v20:expr, $v21:expr, $v22:expr$(;)?
  ) => {
    $crate::Matrix3::new ($v00, $v01, $v02, $v10, $v11, $v12, $v20, $v21, $v22)
  },
  ( $v00:expr, $v01:expr, $v02:expr, $v03:expr;
    $v10:expr, $v11:expr, $v12:expr, $v13:expr;
    $v20:expr, $v21:expr, $v22:expr, $v23:expr;
    $v30:expr, $v31:expr, $v32:expr, $v33:expr$(;)?
  ) => {
    $crate::Matrix4::new (
      $v00, $v01, $v02, $v03,
      $v10, $v11, $v12, $v13,
      $v20, $v21, $v22, $v23,
      $v30, $v31, $v32, $v33)
  }
}

macro_rules! impl_angle {
  ( $angle:ident, $docname:expr ) => {
    #[doc = $docname]
    #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
    // TODO: derive From here causes a test failure in the intrinsic angles test
    #[derive(Clone, Copy, Debug, Default, Eq, PartialEq, PartialOrd)]
    pub struct $angle <S> (pub S);
    impl_numcast_primitive!($angle);
    impl <S : Ring> std::iter::Sum for $angle <S> {
      fn sum <I> (iter : I) -> Self where I : Iterator <Item = Self> {
        use num::Zero;
        iter.fold ($angle::zero(), |acc, item| acc + item)
      }
    }
    impl <S : Ring> std::ops::Neg for $angle <S> {
      type Output = Self;
      fn neg (self) -> Self::Output {
        $angle (-self.0)
      }
    }
    impl <S : Ring> std::ops::Add for $angle <S> {
      type Output = Self;
      fn add (self, other : Self) -> Self::Output {
        $angle (self.0 + other.0)
      }
    }
    impl <S : Ring> std::ops::AddAssign for $angle <S> {
      fn add_assign (&mut self, other : Self) {
        self.0 += other.0
      }
    }
    impl <S : Ring> std::ops::Sub for $angle <S> {
      type Output = Self;
      fn sub (self, other : Self) -> Self::Output {
        $angle (self.0 - other.0)
      }
    }
    impl <S : Ring> std::ops::SubAssign for $angle <S> {
      fn sub_assign (&mut self, other : Self) {
        self.0 -= other.0
      }
    }
    impl <S : Ring> std::ops::Mul <S> for $angle <S> {
      type Output = Self;
      fn mul (self, scalar : S) -> Self::Output {
        $angle (scalar * self.0)
      }
    }
    impl <S : Ring> std::ops::MulAssign <S> for $angle <S> {
      fn mul_assign (&mut self, scalar : S) {
        self.0 *= scalar
      }
    }
    impl <S : Field> std::ops::Div for $angle <S> {
      type Output = S;
      fn div (self, other : Self) -> Self::Output {
        self.0 / other.0
      }
    }
    impl <S : Field> std::ops::Div <S> for $angle <S> {
      type Output = Self;
      fn div (self, scalar : S) -> Self::Output {
        $angle (self.0 / scalar)
      }
    }
    impl <S : Field> std::ops::DivAssign <S> for $angle <S> {
      fn div_assign (&mut self, scalar : S) {
        self.0 /= scalar
      }
    }
    impl <S : OrderedField> std::ops::Rem for $angle <S> {
      type Output = Self;
      fn rem (self, other : Self) -> Self::Output {
        $angle (self.0 % other.0)
      }
    }
    impl <S : Ring> num::Zero for $angle <S> {
      fn zero() -> Self {
        $angle (S::zero())
      }
      fn is_zero (&self) -> bool {
        self.0.is_zero()
      }
    }
    impl <S : Ring> From <S> for $angle <S> {
      fn from (s : S) -> Self {
        $angle (s)
      }
    }
  }
}

macro_rules! impl_angle_wrapped {
  ( $angle:ident, $comment:expr ) => {
    #[doc = $comment]
    #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
    #[derive(Clone, Copy, Debug, Default, Eq, PartialEq, PartialOrd)]
    pub struct $angle <S> (Rad <S>);

    impl_numcast!($angle);
    impl <S : Real> $angle <S> {
      pub const fn angle (&self) -> Rad <S> {
        self.0
      }
    }
    impl <S : Real> std::iter::Sum for $angle <S> {
      fn sum <I> (iter : I) -> Self where I : Iterator <Item = Self> {
        use num::Zero;
        iter.fold ($angle::zero(), |acc, item| acc + item)
      }
    }
    impl <S : Real> std::ops::Neg for $angle <S> {
      type Output = Self;
      fn neg (self) -> Self::Output {
        self.map (Rad::neg)
      }
    }
    impl <S : Real> std::ops::Add for $angle <S> {
      type Output = Self;
      fn add (self, other : Self) -> Self::Output {
        self.map (|angle| angle + other.0)
      }
    }
    impl <S : Real> std::ops::AddAssign for $angle <S> {
      fn add_assign (&mut self, other : Self) {
        *self = *self + other
      }
    }
    impl <S : Real> std::ops::Add <Rad <S>> for $angle <S> {
      type Output = Self;
      fn add (self, angle : Rad <S>) -> Self::Output {
        self.map (|a| a + angle)
      }
    }
    impl <S : Real> std::ops::AddAssign <Rad <S>> for $angle <S> {
      fn add_assign (&mut self, angle : Rad <S>) {
        *self = *self + angle
      }
    }
    impl <S : Real> std::ops::Sub for $angle <S> {
      type Output = Self;
      fn sub (self, other : Self) -> Self::Output {
        self.map (|angle| angle - other.0)
      }
    }
    impl <S : Real> std::ops::SubAssign for $angle <S> {
      fn sub_assign (&mut self, other : Self) {
        *self = *self - other
      }
    }
    impl <S : Real> std::ops::Sub <Rad <S>> for $angle <S> {
      type Output = Self;
      fn sub (self, angle : Rad <S>) -> Self::Output {
        self.map (|a| a - angle)
      }
    }
    impl <S : Real> std::ops::SubAssign <Rad <S>> for $angle <S> {
      fn sub_assign (&mut self, angle : Rad <S>) {
        *self = *self - angle
      }
    }
    impl <S : Real> std::ops::Mul <S> for $angle <S> {
      type Output = Self;
      fn mul (self, scalar : S) -> Self::Output {
        self.map (|angle| angle * scalar)
      }
    }
    impl <S : Real> std::ops::MulAssign <S> for $angle <S> {
      fn mul_assign (&mut self, scalar : S) {
        *self = *self * scalar
      }
    }
    impl <S : Real> std::ops::Div for $angle <S> {
      type Output = S;
      fn div (self, other : Self) -> Self::Output {
        self.0 / other.0
      }
    }
    impl <S : Real> std::ops::Div <S> for $angle <S> {
      type Output = Self;
      fn div (self, scalar : S) -> Self::Output {
        self.map (|angle| angle / scalar)
      }
    }
    impl <S : Real> std::ops::DivAssign <S> for $angle <S> {
      fn div_assign (&mut self, scalar : S) {
        *self = *self / scalar
      }
    }
    impl <S : Real> std::ops::Rem for $angle <S> {
      type Output = Self;
      fn rem (self, other : Self) -> Self::Output {
        self.map (|angle| angle % other.0)
      }
    }
    impl <S : Real> num::Zero for $angle <S> {
      fn zero() -> Self {
        $angle (Rad::zero())
      }
      fn is_zero (&self) -> bool {
        self.0.is_zero()
      }
    }
  }
}

macro_rules! impl_dimension {
  ( $point:ident, $vector:ident, $nonzero:ident, $unit:ident,
    $point_fun:ident, $vector_fun:ident, $matrix_fun:ident,
    [$($component:ident),+],
    [$(($axis_method:ident, $unit_method:ident)),+], [$($patch:ident)?],
    [$($projective:ident)?],
    $matrix:ident, [$($submatrix:ident)?], $ndims:expr, $dimension:expr
  ) => {
    #[doc = $dimension]
    #[doc = "position"]
    #[derive(Clone, Copy, Debug, Default, Eq, PartialEq, Display, From, Neg)]
    #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
    #[display("{}", _0)]
    #[repr(C)]
    pub struct $point <S> (pub $vector <S>);

    #[doc = $dimension]
    #[doc = "non-zero vector"]
    #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
    #[derive(Clone, Copy, Debug, Eq, PartialEq, Display)]
    #[display("{}", _0)]
    #[repr(C)]
    pub struct $nonzero <S> ($vector <S>);

    #[doc = $dimension]
    #[doc = "unit vector"]
    #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
    #[derive(Clone, Copy, Debug, Eq, PartialEq, Display)]
    #[display("{}", _0)]
    #[repr(C)]
    pub struct $unit <S> ($vector <S>);

    //
    //  impl PointN
    //
    #[doc = "Construct a"]
    #[doc = $dimension]
    #[doc = "point"]
    pub const fn $point_fun <S> ($($component : S),+) -> $point <S> {
      $point::new ($($component),+)
    }
    impl_numcast!($point);
    impl <S> $point <S> {
      pub const fn new ($($component : S),+) -> Self {
        $point ($vector::new ($($component),+))
      }
      $(
      pub const fn $component (&self) -> S where S : Copy {
        self.0.$component
      }
      )+
      pub fn distance (self, other : $point <S>) -> S where S : num::real::Real {
        self.0.distance (other.0)
      }
      pub fn distance_squared (self, other : $point <S>) -> S where S : Field {
        self.0.distance_squared (other.0)
      }
    }
    // projective completion
    $(
    impl <S : Field> ProjectiveSpace <S> for $projective <S> {
      type Patch = $point <S>;
    }
    )?
    impl <S : Field> AffineSpace <S> for $point <S> {
      type Translation = $vector <S>;
    }
    impl <S : Ring> Point <$vector <S>> for $point <S> {
      fn to_vector (self) -> $vector <S> {
        self.0
      }
      fn from_vector (vector : $vector <S>) -> Self {
        $point (vector)
      }
    }
    impl <S> From <[S; $ndims]> for $point <S> {
      fn from (array : [S; $ndims]) -> Self {
        $point (array.into())
      }
    }
    $(
    impl <S> From <($patch <S>, S)> for $point <S> {
      fn from ((patch, scale) : ($patch <S>, S)) -> Self {
        $point ((patch.0, scale).into())
      }
    }
    )?
    impl <S> std::cmp::PartialOrd for $point <S> where S : PartialOrd {
      fn partial_cmp (&self, other : &Self) -> Option <std::cmp::Ordering> {
        self.0.as_slice().partial_cmp (&other.0.as_slice())
      }
    }
    impl <S> std::ops::Add <$vector <S>> for $point <S> where S : AdditiveGroup {
      type Output = Self;
      fn add (self, displacement : $vector <S>) -> Self::Output {
        $point (self.0 + displacement)
      }
    }
    impl <S> std::ops::Add <&$vector <S>> for $point <S> where S : AdditiveGroup + Copy {
      type Output = Self;
      fn add (self, displacement : &$vector <S>) -> Self::Output {
        $point (self.0 + *displacement)
      }
    }
    impl <S> std::ops::Add <$vector <S>> for &$point <S> where S : AdditiveGroup + Copy {
      type Output = $point <S>;
      fn add (self, displacement : $vector <S>) -> Self::Output {
        $point (self.0 + displacement)
      }
    }
    impl <S> std::ops::Add <&$vector <S>> for &$point <S> where
      S : AdditiveGroup + Copy
    {
      type Output = $point <S>;
      fn add (self, displacement : &$vector <S>) -> Self::Output {
        $point (self.0 + *displacement)
      }
    }
    impl <S> std::ops::AddAssign <$vector <S>> for $point <S> where S : AdditiveGroup {
      fn add_assign (&mut self, displacement : $vector <S>) {
        self.0 += displacement
      }
    }
    impl <S> std::ops::AddAssign <&$vector <S>> for $point <S> where
      S : AdditiveGroup + Copy
    {
      fn add_assign (&mut self, displacement : &$vector <S>) {
        self.0 += *displacement
      }
    }
    impl <S> std::ops::Sub <$vector <S>> for $point <S> where S : AdditiveGroup {
      type Output = Self;
      fn sub (self, displacement : $vector <S>) -> Self::Output {
        $point (self.0 - displacement)
      }
    }
    impl <S> std::ops::Sub <&$vector <S>> for $point <S> where S : AdditiveGroup + Copy {
      type Output = Self;
      fn sub (self, displacement : &$vector <S>) -> Self::Output {
        $point (self.0 - *displacement)
      }
    }
    impl <S> std::ops::Sub <$point <S>> for $point <S> where S : AdditiveGroup {
      type Output = $vector <S>;
      fn sub (self, other : Self) -> Self::Output {
        self.0 - other.0
      }
    }
    impl <S> std::ops::Sub <&$point <S>> for $point <S> where S : AdditiveGroup + Copy {
      type Output = $vector <S>;
      fn sub (self, other : &Self) -> Self::Output {
        self.0 - other.0
      }
    }
    impl <S> std::ops::Sub <$vector <S>> for &$point <S> where S : AdditiveGroup + Copy {
      type Output = $point <S>;
      fn sub (self, displacement : $vector <S>) -> Self::Output {
        $point (self.0 - displacement)
      }
    }
    impl <S> std::ops::Sub <&$vector <S>> for &$point <S> where
      S : AdditiveGroup + Copy
    {
      type Output = $point <S>;
      fn sub (self, displacement : &$vector <S>) -> Self::Output {
        $point (self.0 - *displacement)
      }
    }
    impl <S> std::ops::Sub <$point <S>> for &$point <S> where S : AdditiveGroup + Copy {
      type Output = $vector <S>;
      fn sub (self, other : $point <S>) -> Self::Output {
        self.0 - other.0
      }
    }
    impl <S> std::ops::Sub <&$point <S>> for &$point <S> where S : AdditiveGroup + Copy {
      type Output = $vector <S>;
      fn sub (self, other : &$point <S>) -> Self::Output {
        self.0 - other.0
      }
    }
    impl <S> std::ops::SubAssign <$vector <S>> for $point <S> where
      S : AdditiveGroup
    {
      fn sub_assign (&mut self, displacement : $vector <S>) {
        self.0 -= displacement
      }
    }
    impl <S> std::ops::SubAssign <&$vector <S>> for $point <S> where
      S : AdditiveGroup + Copy
    {
      fn sub_assign (&mut self, displacement : &$vector <S>) {
        self.0 -= *displacement
      }
    }
    impl <S> approx::AbsDiffEq for $point <S> where
      S : approx::AbsDiffEq,
      S::Epsilon : Copy
    {
      type Epsilon = S::Epsilon;
      fn default_epsilon() -> Self::Epsilon {
        S::default_epsilon()
      }
      #[inline]
      fn abs_diff_eq (&self, other : &Self, epsilon : Self::Epsilon) -> bool {
        self.0.abs_diff_eq (&other.0, epsilon)
      }
    }
    impl <S> approx::RelativeEq for $point <S> where
      S : approx::RelativeEq,
      S::Epsilon : Copy
    {
      fn default_max_relative() -> Self::Epsilon {
        S::default_max_relative()
      }
      #[inline]
      fn relative_eq (&self,
        other : &Self, epsilon : Self::Epsilon, max_relative : Self::Epsilon
      ) -> bool {
        self.0.relative_eq (&other.0, epsilon, max_relative)
      }
    }
    impl <S> approx::UlpsEq for $point <S> where
      S : approx::UlpsEq,
      S::Epsilon : Copy
    {
      fn default_max_ulps() -> u32 {
        S::default_max_ulps()
      }
      #[inline]
      fn ulps_eq (&self, other : &Self, epsilon : Self::Epsilon, max_ulps : u32)
        -> bool
      {
        self.0.ulps_eq (&other.0, epsilon, max_ulps)
      }
    }
    //
    //  impl VectorN
    //
    #[doc = "Construct a"]
    #[doc = $dimension]
    #[doc = "vector"]
    pub const fn $vector_fun <S> ($($component : S),+) -> $vector <S> {
      $vector::new ($($component),+)
    }
    impl <S> InnerProductSpace <S> for $vector <S> where S : Field {
      fn outer_product (self, rhs : Self) -> $matrix <S> {
        let mut cols : $vector <$vector <S>> = [$vector::zero(); $ndims].into();
        for col in 0..$ndims {
          for row in 0..$ndims {
            cols[col][row] = self[row] * rhs[col];
          }
        }
        $matrix { cols }
      }
    }
    impl <S : Ring> Dot <S> for $vector <S> {
      fn dot (self, other : Self) -> S {
        self.dot (other)
      }
    }
    impl <S : Field> VectorSpace <S> for $vector <S> {
      fn map <F> (self, f : F) -> Self where F : FnMut (S) -> S {
        self.map (f)
      }
    }
    impl <S> Module <S> for $vector <S> where S : Ring + num::MulAdd {
      type LinearEndo = $matrix <S>;
    }
    impl <S> GroupAction <$point <S>> for $vector <S> where S : AdditiveGroup {
      fn action (self, point : $point <S>) -> $point <S> {
        point + self
      }
    }
    impl <S : AdditiveGroup> Group for $vector <S> {
      fn identity() -> Self {
        Self::zero()
      }
      fn operation (a : Self, b : Self) -> Self {
        a + b
      }
    }
    impl <S> std::ops::Mul <LinearAuto <S, Self>> for $vector <S> where
      S : Ring + MaybeSerDes
    {
      type Output = Self;
      fn mul (self, rhs : LinearAuto <S, Self>) -> Self::Output {
        self * rhs.0
      }
    }
    impl <S> std::ops::Mul <LinearIso <S, Self, Self, $matrix <S>>> for $vector <S> where
      S : Ring + MaybeSerDes
    {
      type Output = Self;
      fn mul (self, rhs : LinearIso <S, Self, Self, $matrix <S>>) -> Self::Output {
        self * rhs.0
      }
    }
    impl <S> num::Inv for LinearAuto <S, $vector <S>> where
      S : Ring + num::real::Real + MaybeSerDes
    {
      type Output = Self;
      fn inv (self) -> Self::Output {
        LinearAuto (self.0.inv())
      }
    }
    impl <S> std::ops::Div for LinearAuto <S, $vector <S>> where
      S : Ring + num::Float + MaybeSerDes
    {
      type Output = Self;
      fn div (self, rhs : Self) -> Self::Output {
        LinearAuto (self.0 / rhs.0)
      }
    }
    impl <S> std::ops::DivAssign <Self> for LinearAuto <S, $vector <S>> where
      S : Ring + num::Float + MaybeSerDes
    {
      fn div_assign (&mut self, rhs : Self) {
        self.0 /= rhs.0
      }
    }
    //
    //  impl MatrixN
    //
    #[doc = "Construct a"]
    #[doc = $dimension]
    #[doc = "matrix"]
    pub const fn $matrix_fun <S> ($($component : $vector <S>),+) -> $matrix <S> {
      $matrix {
        cols: $vector_fun ($($component),+)
      }
    }
    impl <S, V, W> num::Inv for LinearIso <S, V, W, $matrix <S>> where
      V : Module <S>,
      W : Module <S>,
      S : Ring + num::real::Real
    {
      type Output = LinearIso <S, W, V, $matrix <S>>;
      fn inv (self) -> Self::Output {
        // TODO: currently the vek crate only implements invert for Mat4
        let mat4 = Matrix4::from (self.0);
        LinearIso (mat4.inverted().into(), PhantomData::default())
      }
    }
    impl <S, V> std::ops::Div for LinearIso <S, V, V, $matrix <S>> where
      V : Module <S, LinearEndo=$matrix <S>>,
      S : Ring + num::real::Real
    {
      type Output = Self;
      #[expect(clippy::suspicious_arithmetic_impl)]
      fn div (self, rhs : Self) -> Self::Output {
        use num::Inv;
        self * rhs.inv()
      }
    }
    impl <S, V> std::ops::DivAssign for LinearIso <S, V, V, $matrix <S>> where
      V : Module <S, LinearEndo=$matrix <S>>,
      S : Ring + num::real::Real
    {
      #[expect(clippy::suspicious_op_assign_impl)]
      fn div_assign (&mut self, rhs : Self) {
        use num::Inv;
        *self *= rhs.inv()
      }
    }
    $(
    impl <S : Copy> Matrix <S> for $matrix <S> {
      type Rows = $vector <$vector <S>>;
      type Submatrix = $submatrix <S>;
      fn rows (self) -> $vector <$vector <S>> {
        self.into_row_arrays().map ($vector::from).into()
      }
      fn submatrix (self, i : usize, j : usize) -> $submatrix <S> {
        let a : [S; ($ndims-1) * ($ndims-1)] = std::array::from_fn (|index|{
          let i = (i + 1 + index / ($ndims-1)) % $ndims;
          let j = (j + 1 + index % ($ndims-1)) % $ndims;
          self.cols[i][j]
        });
        $submatrix::from_col_array (a)
      }
      /// Fill additional row and column with zeros
      fn fill_zeros (submatrix : $submatrix <S>) -> Self where S : num::Zero{
        let cols = $vector::from (
          std::array::from_fn (|i| if i < $ndims-1 {
            $vector::from (submatrix.cols[i])
          } else {
            $vector::zero()
          })
        );
        $matrix { cols }
      }
    }
    )?
    impl <S : Ring> LinearMap <S, $vector <S>, $vector <S>> for $matrix <S> {
      fn determinant (self) -> S {
        self.determinant()
      }
      fn transpose (self) -> Self {
        self.transposed()
      }
    }
    //
    //  impl NonZeroN
    //
    impl_numcast!($nonzero);
    impl <S : Ring> $nonzero <S> {
      /// Returns 'None' if called with the zero vector
      // TODO: move broken doctests to unit tests
      // # Example
      //
      // ```
      // # use math_utils::$nonzero;
      // assert!($nonzero::new ([1.0, 0.0].into()).is_some());
      // assert!($nonzero::new ([0.0, 0.0].into()).is_none());
      // ```
      pub fn new (vector : $vector <S>) -> Option <Self> {
        use num::Zero;
        if !vector.is_zero() {
          Some ($nonzero (vector))
        } else {
          None
        }
      }
      /// Panics if zero vector is given.
      // TODO: move broken doctests to unit tests
      // # Panics
      //
      // ```should_panic
      // # use math_utils::$nonzero;
      // let x = $nonzero::noisy ([0.0, 0.0].into());  // panic!
      // ```
      pub fn noisy (vector : $vector <S>) -> Self {
        use num::Zero;
        assert!(!vector.is_zero());
        $nonzero (vector)
      }
      /// Map an operation on the underlying vector, panicking if the result is zero
      // TODO: move broken doctests to unit tests
      // # Panics
      //
      // Panics of the result is zero:
      //
      // ```should_panic
      // # use math_utils::$nonzero;
      // let v = $nonzero::noisy ([1.0, 1.0].into())
      //   .map_noisy (|x| 0.0 * x);  // panic!
      // ```
      pub fn map_noisy (self, fun : fn ($vector <S>) -> $vector <S>) -> Self {
        Self::noisy (fun (self.0))
      }
    }
    impl <S : Ring> std::ops::Deref for $nonzero <S> {
      type Target = $vector <S>;
      fn deref (&self) -> &$vector <S> {
        &self.0
      }
    }
    //
    //  impl UnitN
    //
    impl_numcast!($unit);
    impl <S : Real> $unit <S> {
      // TODO: move broken doc tests to unit tests
      // ```
      // # use math_utils::$unit;
      // assert_eq!($unit::axis_x(), $unit::new ([1.0, 0.0].into()).unwrap());
      // ```
      $(
      #[inline]
      pub fn $axis_method() -> Self {
        $unit ($vector::$unit_method())
      }
      )+
      /// Returns 'None' if called with a non-normalized vector.
      ///
      /// This method checks whether `vector == vector.normalize()`.
      pub fn new (vector : $vector <S>) -> Option <Self> {
        let normalized = vector.normalize();
        if vector == normalized {
          Some ($unit (vector))
        } else {
          None
        }
      }
      /// Returns 'None' if called with an non-normalized vector determined by checking
      /// if the magnitude is approximately one.
      ///
      /// Note the required `RelativeEq` trait is only implemented for floating-point
      /// types.
      pub fn new_approx (vector : $vector <S>) -> Option <Self> where
        S : num::real::Real + approx::RelativeEq <Epsilon=S>
      {
        if vector.is_normalized() {
          Some ($unit (vector))
        } else {
          None
        }
      }
      /// Normalizes a given non-zero vector.
      // TODO: move broken doc tests to unit tests
      // # Example
      //
      // ```
      // # use math_utils::$unit;
      // assert_eq!(
      //   *$unit::normalize ([2.0, 0.0].into()),
      //   [1.0, 0.0].into()
      // );
      // ```
      //
      // # Panics
      //
      // Panics if the zero vector is given:
      //
      // ```should_panic
      // # use math_utils::$unit;
      // let x = $unit::normalize ([0.0, 0.0].into());  // panic!
      // ```
      pub fn normalize (vector : $vector <S>) -> Self {
        use num::Zero;
        assert!(!vector.is_zero());
        $unit (vector.normalize())
      }
      /// Normalizes a given non-zero vector only if the vector is not already
      /// normalized.
      pub fn normalize_approx (vector : $vector <S>) -> Self where
        S : num::real::Real + approx::RelativeEq <Epsilon=S>
      {
        use num::Zero;
        assert!(!vector.is_zero());
        let vector = if vector.is_normalized() {
          vector
        } else {
          vector.normalize()
        };
        $unit (vector)
      }
      /// Panics if a non-normalized vector is given.
      ///
      /// This method checks whether `vector == vector.normalize()`.
      // TODO: move broken doc tests to unit tests
      // ```should_panic
      // # use math_utils::$unit;
      // let x = $unit::noisy ([2.0, 0.0].into());  // panic!
      // ```
      pub fn noisy (vector : $vector <S>) -> Self where S : std::fmt::Debug {
        assert_eq!(vector, vector.normalize());
        $unit (vector)
      }
      /// Panics if a non-normalized vector is given.
      ///
      /// Checks if the magnitude is approximately one.
      // TODO: move broken doc tests to unit tests
      // ```should_panic
      // # use math_utils::$unit;
      // let x = $unit::noisy ([2.0, 0.0].into());  // panic!
      // ```
      pub fn noisy_approx (vector : $vector <S>) -> Self where
        S : num::real::Real + approx::RelativeEq <Epsilon=S>
      {
        assert!(vector.is_normalized());
        $unit (vector)
      }
      /// It is a debug assertion if the given vector is not normalized.
      ///
      /// In debug builds, method checks that `vector == vector.normalize()`.
      #[inline]
      pub fn unchecked (vector : $vector <S>) -> Self where S : std::fmt::Debug {
        debug_assert_eq!(vector, vector.normalize());
        $unit (vector)
      }
      /// It is a debug assertion if the given vector is not normalized.
      ///
      /// In debug builds, checks if the magnitude is approximately one.
      // TODO: move broken doc tests to unit tests
      // ```should_panic
      // # use math_utils::$unit;
      // let n = $unit::unchecked ([2.0, 0.0].into());  // panic!
      // ```
      #[inline]
      pub fn unchecked_approx (vector : $vector <S>) -> Self where
        S : num::real::Real + approx::RelativeEq <Epsilon=S>
      {
        debug_assert!(vector.is_normalized());
        $unit (vector)
      }
      /// Generate a random unit vector. Uses naive algorithm.
      pub fn random_unit <R : rand::Rng> (rng : &mut R) -> Self where
        S : rand::distr::uniform::SampleUniform
      {
        let vector = $vector {
          $($component: rng.random_range (-S::one()..S::one())),+
        };
        $unit (vector.normalize())
      }
      /// Return the unit vector pointing in the opposite direction
      // TODO: implement num::Inv ?
      pub fn invert (mut self) -> Self {
        self.0 = -self.0;
        self
      }
      /// Map an operation on the underlying vector, panicking if the result is zero.
      ///
      /// # Panics
      ///
      /// Panics of the result is not normalized.
      // TODO: move broken doc tests to unit tests
      // ```should_panic
      // # use math_utils::$unit;
      // // panic!
      // let v = $unit::noisy ([1.0, 0.0].into()).map_noisy (|x| 2.0 * x);
      // ```
      pub fn map_noisy (self, fun : fn ($vector <S>) -> $vector <S>) -> Self where
        S : std::fmt::Debug
      {
        Self::noisy (fun (self.0))
      }
      /// Map an operation on the underlying vector, panicking if the result is zero.
      ///
      /// # Panics
      ///
      /// Panics of the result is not normalized.
      // TODO: move broken doc tests to unit tests
      // ```should_panic
      // # use math_utils::$unit;
      // // panic!
      // let v = $unit::noisy ([1.0, 0.0].into()).map_noisy (|x| 2.0 * x);
      // ```
      pub fn map_noisy_approx (self, fun : fn ($vector <S>) -> $vector <S>) -> Self where
        S : num::real::Real + approx::RelativeEq <Epsilon=S>
      {
        Self::noisy_approx (fun (self.0))
      }
    }
    impl <S : Real> std::ops::Deref for $unit <S> {
      type Target = $vector <S>;
      fn deref (&self) -> &$vector <S> {
        &self.0
      }
    }
    impl <S : Real> std::ops::Neg for $unit <S> {
      type Output = Self;
      fn neg (self) -> Self::Output {
        Self (-self.0)
      }
    }
    //  end UnitN
  }
}

macro_rules! impl_numcast {
  ($type:ident) => {
    impl <S> $type <S> {
      #[inline]
      pub fn numcast <T> (self) -> Option <$type <T>> where
        S : num::NumCast,
        T : num::NumCast
      {
        self.0.numcast().map ($type)
      }
    }
  }
}

macro_rules! impl_numcast_primitive {
  ($type:ident) => {
    impl <S> $type <S> {
      #[inline]
      pub fn numcast <T> (self) -> Option <$type <T>> where
        S : num::NumCast,
        T : num::NumCast
      {
        T::from (self.0).map ($type)
      }
    }
  }
}

#[doc(hidden)]
#[macro_export]
macro_rules! impl_numcast_fields {
  ($type:ident, $($field:ident),+) => {
    impl <S> $type <S> {
      #[inline]
      pub fn numcast <T> (self) -> Option <$type <T>> where
        S : num::NumCast,
        T : num::NumCast
      {
        Some ($type {
          $($field: self.$field.numcast()?),+
        })
      }
    }
  }
}

impl_angle!(Deg, "Degrees");
impl_angle!(Rad, "Radians");
impl_angle!(Turn, "Turns");
impl_angle_wrapped!(AngleWrapped, "Unsigned wrapped angle restricted to $[0, 2\\pi)$");
impl_angle_wrapped!(AngleWrappedSigned,
  "Signed wrapped angle restricted to $(-\\pi, \\pi]$");

impl_dimension!(Point2, Vector2, NonZero2, Unit2, point2, vector2, matrix2, [x, y],
  [(axis_x, unit_x), (axis_y, unit_y)],
  [], [Point3], Matrix2, [], 2, "2D");
impl_dimension!(Point3, Vector3, NonZero3, Unit3, point3, vector3, matrix3, [x, y, z],
  [(axis_x, unit_x), (axis_y, unit_y), (axis_z, unit_z)],
  [Point2], [Point4], Matrix3, [Matrix2], 3, "3D");
impl_dimension!(Point4, Vector4, NonZero4, Unit4, point4, vector4, matrix4, [x, y, z, w],
  [(axis_x, unit_x), (axis_y, unit_y), (axis_z, unit_z), (axis_w, unit_w)],
  [Point3], [], Matrix4, [Matrix3], 4, "4D");

impl_numcast_primitive!(Positive);
impl_numcast_primitive!(NonNegative);
impl_numcast_primitive!(NonZero);
impl_numcast_primitive!(Normalized);
impl_numcast_primitive!(NormalSigned);
impl_numcast_fields!(Angles3, yaw, pitch, roll);

//
//  private
//
#[inline]
fn normalize_quaternion <S : Real> (quaternion : Quaternion <S>) -> Quaternion <S> {
  Quaternion::from_vec4 (quaternion.into_vec4().normalize())
}

#[cfg(test)]
mod tests {
  use super::*;
  use crate::num_traits;
  use approx;
  use rand;
  use rand_xorshift;

  #[test]
  fn defaults() {
    use num_traits::Zero;

    assert_eq!(Positive::<f32>::default().0, 0.0);
    assert_eq!(NonNegative::<f32>::default().0, 0.0);
    assert_eq!(Normalized::<f32>::default().0, 0.0);
    assert_eq!(NormalSigned::<f32>::default().0, 0.0);
    assert_eq!(Deg::<f32>::default().0, 0.0);
    assert_eq!(Rad::<f32>::default().0, 0.0);
    assert_eq!(Turn::<f32>::default().0, 0.0);
    assert_eq!(
      Angles3::<f32>::default(),
      Angles3::new (AngleWrapped::zero(), AngleWrapped::zero(), AngleWrapped::zero()));
    assert_eq!(
      Pose3::<f32>::default(),
      Pose3 { position: Point3::origin(), angles: Angles3::default() });
    assert_eq!(
      LinearIso::<f32, Vector3 <f32>, Vector3 <f32>, Matrix3 <f32>>::default().0,
      Matrix3::identity());
    assert_eq!(LinearAuto::<f32, Vector3 <f32>>::default().0.0, Matrix3::identity());
    assert_eq!(Rotation2::<f32>::default().0.0.0, Matrix2::identity());
    assert_eq!(Rotation3::<f32>::default().0.0.0, Matrix3::identity());
    assert_eq!(Versor::<f32>::default().0, Quaternion::from_xyzw (0.0, 0.0, 0.0, 1.0));
    assert_eq!(
      AffineMap::<f32, Point3 <f32>, Point3 <f32>, Matrix3 <f32>>::default(),
      AffineMap::new (Matrix3::identity(), Vector3::zero()));
    assert_eq!(
      Affinity::<f32, Point3 <f32>, Point3 <f32>, Matrix3 <f32>>::default(),
      Affinity::new (LinearIso::new (Matrix3::identity()).unwrap(), Vector3::zero()));
    assert_eq!(
      Projectivity::<f32, Point4 <f32>, Point4 <f32>, Matrix4 <f32>>::default().0,
      LinearIso::new (Matrix4::identity()).unwrap());
  }

  #[test]
  fn vector_sum() {
    let s = [
      [0.0, 1.0],
      [1.0, 1.0],
      [0.5, 0.5]
    ].map (Vector2::<f32>::from);
    let _sum = s.iter().copied().sum::<Vector2 <f32>>();
  }

  #[test]
  fn vector_outer_product() {
    let v1 = vector3 (1.0, 2.0, 3.0);
    let v2 = vector3 (2.0, 3.0, 4.0);
    assert_eq!(v1.outer_product (v2).into_row_arrays(), [
      [2.0, 3.0,  4.0],
      [4.0, 6.0,  8.0],
      [6.0, 9.0, 12.0]
    ]);
  }

  #[test]
  fn rotation3_noisy_approx() {
    use rand::{Rng, SeedableRng};
    // construct orthonormal matrices and verify they are orthonormal using noisy_approx
    // constructor
    let mut rng = rand_xorshift::XorShiftRng::seed_from_u64 (0);
    let up = Vector3::unit_z();
    for _ in 0..1000 {
      let forward = vector3 (
        rng.random_range (-1000.0f32..1000.0),
        rng.random_range (-1000.0f32..1000.0),
        0.0
      ).normalized();
      let right = forward.cross (up);
      let mat = Matrix3::from_col_arrays ([
        right.into_array(),
        forward.into_array(),
        up.into_array()
      ]);
      let _ = Rotation3::noisy_approx (mat);
    }
  }

  #[test]
  fn rotation3_intrinsic_angles() {
    use std::f32::consts::PI;
    use num_traits::Zero;
    use rand::{Rng, SeedableRng};
    // yaw: CCW quarter turn around Z (world up) axis
    let rot = Rotation3::from_angles_intrinsic (
      Turn (0.25).into(), Rad::zero(), Rad::zero());
    approx::assert_relative_eq!(rot.mat(),
      Matrix3::from_col_arrays ([
        [ 0.0, 1.0, 0.0],
        [-1.0, 0.0, 0.0],
        [ 0.0, 0.0, 1.0]
      ])
    );
    assert_eq!((Turn (0.25).into(), Rad::zero(), Rad::zero()), rot.intrinsic_angles());
    let rot = Rotation3::from_angles_intrinsic (
      Turn (0.5).into(), Rad::zero(), Rad::zero());
    approx::assert_relative_eq!(rot.mat(),
      Matrix3::from_col_arrays ([
        [-1.0,  0.0, 0.0],
        [ 0.0, -1.0, 0.0],
        [ 0.0,  0.0, 1.0]
      ])
    );
    assert_eq!((Turn (0.5).into(), Rad::zero(), Rad::zero()), rot.intrinsic_angles());
    // pitch: CCW quarter turn around X (world right) axis
    let rot = Rotation3::from_angles_intrinsic (
      Rad::zero(), Turn (0.25).into(), Rad::zero());
    approx::assert_relative_eq!(rot.mat(),
      Matrix3::from_col_arrays ([
        [1.0,  0.0, 0.0],
        [0.0,  0.0, 1.0],
        [0.0, -1.0, 0.0]
      ])
    );
    assert_eq!((Rad::zero(), Turn (0.25).into(), Rad::zero()), rot.intrinsic_angles());
    // roll: CCW quarter turn around Y (world forward) axis
    let rot = Rotation3::from_angles_intrinsic (
      Rad::zero(), Rad::zero(), Turn (0.25).into());
    approx::assert_relative_eq!(rot.mat(),
      Matrix3::from_col_arrays ([
        [0.0, 0.0, -1.0],
        [0.0, 1.0,  0.0],
        [1.0, 0.0,  0.0]
      ])
    );
    assert_eq!((Rad::zero(), Rad::zero(), Turn (0.25).into()), rot.intrinsic_angles());
    // pitch after yaw
    let rot = Rotation3::from_angles_intrinsic (
      Turn (0.25).into(), Turn (0.25).into(), Rad::zero());
    approx::assert_relative_eq!(rot.mat(),
      Matrix3::from_col_arrays ([
        [0.0, 1.0, 0.0],
        [0.0, 0.0, 1.0],
        [1.0, 0.0, 0.0]
      ])
    );
    let angles = rot.intrinsic_angles();
    approx::assert_relative_eq!(Rad::from (Turn (0.25)).0, angles.0.0);
    approx::assert_relative_eq!(Rad::from (Turn (0.25)).0, angles.1.0);
    approx::assert_relative_eq!(0.0, angles.2.0);
    // roll after pitch after yaw
    let rot = Rotation3::from_angles_intrinsic (
      Turn (0.25).into(), Turn (0.25).into(), Turn (0.25).into());
    approx::assert_relative_eq!(rot.mat(),
      Matrix3::from_col_arrays ([
        [-1.0, 0.0, 0.0],
        [ 0.0, 0.0, 1.0],
        [ 0.0, 1.0, 0.0]
      ])
    );
    let angles = rot.intrinsic_angles();
    approx::assert_relative_eq!(Rad::from (Turn (0.25)).0, angles.0.0);
    approx::assert_relative_eq!(Rad::from (Turn (0.25)).0, angles.1.0);
    approx::assert_relative_eq!(Rad::from (Turn (0.25)).0, angles.2.0);

    // construct random orthonormal matrices and verify the intrinsic angles are in the
    // range [-pi, pi]
    let mut rng = rand_xorshift::XorShiftRng::seed_from_u64 (0);
    let mut random_vector = || vector3 (
      rng.random_range (-100.0f32..100.0),
      rng.random_range (-100.0f32..100.0),
      rng.random_range (-100.0f32..100.0));
    for _ in 0..1000 {
      let rot = Rotation3::orthonormalize (
        random_vector(), random_vector(), random_vector()
      ).unwrap();
      let (yaw, pitch, roll) = rot.intrinsic_angles();
      assert!(yaw.0   <  PI);
      assert!(yaw.0   > -PI);
      assert!(pitch.0 <  PI);
      assert!(pitch.0 > -PI);
      assert!(roll.0  <  PI);
      assert!(roll.0  > -PI);
    }
  }

  #[test]
  fn rotation3_into_versor() {
    use std::f32::consts::PI;
    let rot  = Rotation3::from_angle_x (Rad (0.01));
    let quat = rot.versor().0;
    approx::assert_relative_eq!(rot.mat(), quat.into());
    let rot  = Rotation3::from_angle_x (Rad (-0.01));
    let quat = rot.versor().0;
    approx::assert_relative_eq!(rot.mat(), quat.into());

    let rot  = Rotation3::from_angle_y (Rad (0.01));
    let quat = rot.versor().0;
    approx::assert_relative_eq!(rot.mat(), quat.into());
    let rot  = Rotation3::from_angle_y (Rad (-0.01));
    let quat = rot.versor().0;
    approx::assert_relative_eq!(rot.mat(), quat.into());

    let rot  = Rotation3::from_angle_z (Rad (0.01));
    let quat = rot.versor().0;
    approx::assert_relative_eq!(rot.mat(), quat.into());
    let rot  = Rotation3::from_angle_z (Rad (-0.01));
    let quat = rot.versor().0;
    approx::assert_relative_eq!(rot.mat(), quat.into());

    let rot  = Rotation3::from_angle_x (Rad (PI / 2.0));
    let quat = rot.versor().0;
    approx::assert_relative_eq!(rot.mat(), quat.into());
    let rot  = Rotation3::from_angle_x (Rad (-PI / 2.0));
    let quat = rot.versor().0;
    approx::assert_relative_eq!(rot.mat(), quat.into());

    let rot  = Rotation3::from_angle_y (Rad (PI / 2.0));
    let quat = rot.versor().0;
    approx::assert_relative_eq!(rot.mat(), quat.into());
    let rot  = Rotation3::from_angle_y (Rad (-PI / 2.0));
    let quat = rot.versor().0;
    approx::assert_relative_eq!(rot.mat(), quat.into());

    let rot  = Rotation3::from_angle_z (Rad (PI / 2.0));
    let quat = rot.versor().0;
    approx::assert_relative_eq!(rot.mat(), quat.into());
    let rot  = Rotation3::from_angle_z (Rad (-PI / 2.0));
    let quat = rot.versor().0;
    approx::assert_relative_eq!(rot.mat(), quat.into());
  }
}