math_utils/types/

mod.rs

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

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::Mat2 as Matrix2;
pub use vek::Vec3 as Vector3;
pub use vek::Mat3 as Matrix3;
pub use vek::Vec4 as Vector4;
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, 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, 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, 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, PartialEq)]
pub struct Pose3 <S> {
  pub position : Point3  <S>,
  pub angles   : Angles3 <S>
}

/// Invertible linear map
#[derive(Clone, Copy, Debug, Default, 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, PartialEq)]
pub struct Versor <S> (Quaternion <S>) where S : num::Zero + num::One;

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

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

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

////////////////////////////////////////////////////////////////////////////////
//  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::one() / S::two())
  }
}
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 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 fn component (self) -> usize {
    self as usize
  }
}

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

//
//  impl SignedAxis3
//
impl SignedAxis3 {
  pub fn try_from <S : Field> (vector : &Vector3 <S>) -> Option <Self> {
    return 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 : Real> (self) -> Unit3 <S> where S : std::fmt::Debug {
    Unit3::unchecked (self.to_vec())
  }
}

//
//  impl Quadrant
//
impl Quadrant {
  /// Returns some 'Some (quadrant)' if the vector is strictly contained with an
  /// 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 with 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 : Ring> 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!
  /// ```
  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 for Positive <S> {
  type Output = Self;
  fn mul (self, rhs : Self) -> Self::Output {
    Positive (self.0 * rhs.0)
  }
}
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 : Field> std::ops::Div for Positive <S> {
  type Output = Self;
  fn div (self, rhs : Self) -> Self::Output {
    Positive (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::Add <NonNegative <S>> for Positive <S> {
  type Output = Self;
  fn add (self, rhs : NonNegative <S>) -> Self::Output {
    Positive (self.0 + rhs.0)
  }
}
//  end impl Positive

//
//  impl NonNegative
//
impl <S : Ring> 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 : Field> 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]
  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> { }
impl <S : Field> Ord for Normalized <S> {
  fn cmp (&self, rhs : &Self) -> std::cmp::Ordering {
    // safe to unwrap: normalized values are never NaN
    self.partial_cmp (rhs).unwrap()
  }
}
//  end impl Normalized

//
//  impl NormalSigned
//
impl <S : Field> 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]
  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 : Field> 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> { }
impl <S : Field> Ord for NormalSigned <S> {
  fn cmp (&self, rhs : &Self) -> std::cmp::Ordering {
    // safe to unwrap: normalized values are never NaN
    self.partial_cmp (rhs).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 epsilon      = S::default_epsilon() * four;
    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::default()))))
    }
  }
  /// 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::default())))
  }
  /// 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::default())))
  }
  /// 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::default())))
  }
  /// 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::default())))
  }
  /// 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::default())))
  }
}
impl <S> std::ops::Deref for Rotation2 <S> where S : Ring + MaybeSerDes {
  type Target = Matrix2 <S>;
  fn deref (&self) -> &Self::Target {
    &self.0
  }
}
impl <S> MultiplicativeGroup for Rotation2 <S> where
  S : Ring + MaybeSerDes
{ }
impl <S> MultiplicativeMonoid for Rotation2 <S> where
  S : Ring + MaybeSerDes
{ }
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::default()
    )))
  }
}
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;
  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
{
  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::default())))
  }
}

//
//  impl Rotation3
//
impl <S> Rotation3 <S> where S : Ring + MaybeSerDes {
  /// Identity rotation
  pub fn identity() -> Self {
    use num::One;
    Self::one()
  }
  /// 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::default())))
  }
  /// 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::default())))
  }
  /// 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::default())))
  }
  /// 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.into(), pitch) * rotation1;
    Rotation3::from_axis_angle (rotation2.cols.y.into(), 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::default())))
  }
  /// 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::default()))))
    }
  }

  /// 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)
    }
  }

  /// 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::default()))))
    }
  }
  /// 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 epsilon      = S::default_epsilon() * four;
    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::default()))))
    }
  }
  /// 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::default())))
  }
  /// 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::default())))
  }
  /// 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::default())))
  }
  /// 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::default())))
  }
  /// 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::one() / S::two()) * 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)
  }
}
impl <S> std::ops::Deref for Rotation3 <S> where S : Ring + MaybeSerDes {
  type Target = Matrix3 <S>;
  fn deref (&self) -> &Self::Target {
    &self.0
  }
}
impl <S> MultiplicativeGroup  for Rotation3 <S> where S : Ring + MaybeSerDes { }
impl <S> MultiplicativeMonoid for Rotation3 <S> where S : Ring + MaybeSerDes { }
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::default()
    )))
  }
}
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;
  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
{
  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::default())))
  }
}
impl <S> From <Angles3 <S>> for Rotation3 <S> where
  S : Real + MaybeSerDes + num::real::Real
{
  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 : 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> 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 AffineMap
//
impl <S, A, B, M> AffineMap <S, A, B, M> where
  A : AffineSpace <S>,
  B : AffineSpace <S>,
  M : LinearMap <S, A::Vector, B::Vector>,
  S : Field + std::fmt::Display,
  B::Vector : std::fmt::Display
{
  pub fn new (linear_map : M, translation : B::Vector) -> Self {
    AffineMap { linear_map, translation, _phantom: PhantomData::default() }
  }
  pub fn transform (self, point : A::Point) -> B::Point {
    B::Point::from_vector (
      self.linear_map * point.to_vector() + self.translation)
  }
}

//
//  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
{
  /// 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) -> bool where
    S : approx::RelativeEq <Epsilon=S>
  {
    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::relative_ne!(linear_map.determinant(), S::zero(),
      max_relative=max_relative, 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::default()))
    } else {
      None
    }
  }
  /// Checks if map is invertible with approximate equality
  pub fn new_approx (linear_map : M) -> Option <Self> where
    S : approx::RelativeEq <Epsilon=S>
  {
    if Self::is_invertible_approx (linear_map) {
      Some (LinearIso (linear_map, PhantomData::default()))
    } 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::default())
  }
}
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::default())
  }
}
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
{
  /// 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> MultiplicativeMonoid for LinearAuto <S, V> where
  V : Module <S>,
  V::LinearEndo : MaybeSerDes + PartialEq,
  S : Ring
{ }
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 Affinity
//
impl <S, A, B, M> Affinity <S, A, B, M> where
  M : LinearMap <S, A::Vector, B::Vector>,
  A : AffineSpace <S>,
  B : AffineSpace <S>,
  S : Field + std::fmt::Display,
  B::Vector : std::fmt::Display
{
  pub fn new (
    linear_iso  : LinearIso <S, A::Vector, B::Vector, M>,
    translation : B::Vector
  ) -> Self {
    Affinity { linear_iso, translation, _phantom: PhantomData::default() }
  }
}

//
//  impl Projectivity
//
impl <S, V, W, P, Q, M> Projectivity <S, V, W, P, Q, M> where
  M : LinearMap <S, P, Q>,
  P : ProjectiveSpace <S, V>,
  Q : ProjectiveSpace <S, W>,
  V : VectorSpace <S> + std::fmt::Display,
  W : VectorSpace <S> + std::fmt::Display,
  S : Field + std::fmt::Display
{
  pub fn new (linear_iso : LinearIso <S, P, Q, M>) -> Self {
    Projectivity (linear_iso, PhantomData::default())
  }
  pub fn transform (self, vector : P) -> Q {
    **self * vector
  }
}
impl <S, V, W, P, Q, M> std::ops::Deref for Projectivity <S, V, W, P, Q, M> where
  M : LinearMap <S, P, Q>,
  P : ProjectiveSpace <S, V>,
  Q : ProjectiveSpace <S, W>,
  V : VectorSpace <S> + std::fmt::Display,
  W : VectorSpace <S> + std::fmt::Display,
  S : Field + std::fmt::Display
{
  type Target = LinearIso <S, P, Q, M>;
  fn deref (&self) -> &Self::Target {
    &self.0
  }
}

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, PartialEq, PartialOrd)]
    pub struct $angle <S> (pub S);
    impl_numcast_primitive!($angle);
    impl <S : Real> AdditiveGroup  for $angle <S> { }
    impl <S : Real> AdditiveMonoid for $angle <S> { }
    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 : Ring> std::ops::Div for $angle <S> {
      type Output = S;
      fn div (self, other : Self) -> Self::Output {
        self.0 / other.0
      }
    }
    impl <S : Ring> 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 : Ring> 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()
      }
    }
  }
}

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

    impl_numcast!($angle);
    impl <S : Real> $angle <S> {
      pub fn angle (&self) -> Rad <S> {
        self.0
      }
    }
    impl <S : Real> AdditiveGroup  for $angle <S> { }
    impl <S : Real> AdditiveMonoid for $angle <S> { }
    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,
    [$($component:ident),+],
    [$(($axis_method:ident, $unit_method:ident)),+], [$($projective:ident)?],
    $matrix:ident, $affine:ident, $euclidean: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 = "affine space"]
    #[derive(Debug)]
    pub struct $affine <S> (PhantomData <S>);
    impl <S : Field> AffineSpace <S> for $affine <S> {
      type Point  = $point  <S>;
      type Vector = $vector <S>;
    }

    #[doc = $dimension]
    #[doc = "euclidean space"]
    #[derive(Debug)]
    pub struct $euclidean <S> (PhantomData <S>);
    impl <S : Real> EuclideanSpace <S> for $euclidean <S> { }
    impl <S : Real> AffineSpace <S> for $euclidean <S> {
      type Point  = $point  <S>;
      type Vector = $vector <S>;
    }

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

    //
    //  impl PointN
    //
    impl_numcast!($point);
    impl <S> $point <S> {
      pub const fn new ($($component : S),+) -> Self {
        $point ($vector::new ($($component),+))
      }
    }
    impl <S> Point <$vector <S>> for $point <S> where S : Ring {
      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> 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
    //
    impl <S : Field> InnerProductSpace <S> for $vector <S> { }
    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 : Ring> Module <S> for $vector <S> {
      type LinearEndo = $matrix <S>;
    }
    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 <S : Ring> MultiplicativeMonoid for $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 : Ring> AdditiveGroup  for $vector <S> { }
    impl <S : Ring> AdditiveMonoid for $vector <S> { }
    impl <S : AdditiveGroup> Group for $vector <S> {
      fn identity() -> Self {
        Self::zero()
      }
      fn operation (a : Self, b : Self) -> Self {
        a + b
      }
    }
    impl <S> Point <$vector <S>> for $vector <S> where S : Ring {
      fn to_vector (self) -> $vector <S> {
        self
      }
      fn from_vector (vector : $vector <S>) -> Self {
        vector
      }
    }
    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::Float + 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 <S> MultiplicativeGroup for LinearAuto <S, $vector <S>> where
      S : Ring + num::Float + MaybeSerDes
    { }
    // projective completion
    $(
      impl <S> ProjectiveSpace <S, $vector <S>> for $projective <S> where
        S : Field + std::fmt::Display
      {
        fn homography <A> (
          affinity : Affinity <S, A, A, <A::Vector as Module <S>>::LinearEndo>
        ) -> Projectivity <
          S, $vector <S>, $vector <S>, Self, Self, Self::LinearEndo
        > where
          A : AffineSpace <S, Vector=$vector <S>>
        {
          Projectivity::new (
            LinearIso::<S, Self, Self, Self::LinearEndo>::new (
              (*affinity.linear_iso).into()
            ).unwrap()
          )
        }
        fn homogeneous <A> (point_or_vector : Either <A::Point, A::Vector>)
          -> Self
        where
          A : AffineSpace <S, Vector=$vector <S>>,
          $vector <S> : GroupAction <A::Point>
        {
          match point_or_vector {
            Either::Left  (point)  => (point.to_vector(), S::one()).into(),
            Either::Right (vector) => (vector, S::zero()).into(),
          }
        }
      }
    )?
    //
    //  impl MatrixN
    //
    impl <S, V, W> num::Inv for LinearIso <S, V, W, $matrix <S>> where
      V : Module <S>,
      W : Module <S>,
      S : Ring + num::Float
    {
      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::Float
    {
      type Output = Self;
      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::Float
    {
      fn div_assign (&mut self, rhs : Self) {
        use num::Inv;
        *self *= rhs.inv()
      }
    }
    //
    //  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)
      }
    }
  }
}

#[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,
  [x, y],
  [(axis_x, unit_x), (axis_y, unit_y)], [Vector3], Matrix2, Affine2, Euclidean2,
  2, "2D");
impl_dimension!(Point3, Vector3, NonZero3, Unit3,
  [x, y, z],
  [(axis_x, unit_x), (axis_y, unit_y), (axis_z, unit_z)],
  [Vector4], Matrix3, Affine3, Euclidean3, 3, "3D");
impl_dimension!(Point4, Vector4, NonZero4, Unit4,
  [x, y, z, w],
  [(axis_x, unit_x), (axis_y, unit_y), (axis_z, unit_z), (axis_w, unit_w)],
  [], Matrix4, Affine4, Euclidean4, 4, "4D");

impl_numcast_primitive!(Positive);
impl_numcast_primitive!(NonNegative);
impl_numcast_primitive!(NonZero);
impl_numcast_primitive!(Normalized);
impl_numcast_primitive!(NormalSigned);

//
//  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, Vector3 <f32>, Vector3 <f32>, Matrix3 <f32>>::default(),
      AffineMap::new (Matrix3::identity(), Vector3::zero()));
    assert_eq!(
      Affinity::<f32, Vector3 <f32>, Vector3 <f32>, Matrix3 <f32>>::default(),
      Affinity::new (
        LinearIso::new (Matrix3::identity()).unwrap(), Vector3::zero()));
    assert_eq!(
      Projectivity::<
        f32, Vector3 <f32>, Vector3 <f32>, Vector4 <f32>, Vector4 <f32>,
        Matrix4 <f32>
      >::default().0,
      LinearIso::new (Matrix4::identity()).unwrap());
  }

  #[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::new (
        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,
      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,
      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,
      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,
      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,
      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,
      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::new (
      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, quat.into());
    let rot  = Rotation3::from_angle_x (Rad (-0.01));
    let quat = rot.versor().0;
    approx::assert_relative_eq!(*rot, quat.into());

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

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

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

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

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