From 637903a30dd3bade46075a96661bae7db060faac Mon Sep 17 00:00:00 2001 From: ickk Date: Sat, 4 May 2024 15:13:01 +1000 Subject: [PATCH] impl `Outermorphism` for `Matrix4x4` --- ega/src/lib.rs | 39 +++ ega/src/operators/add.rs | 1 - ega/src/operators/dot.rs | 1 - ega/src/operators/geometric_product.rs | 1 - ega/src/operators/join.rs | 1 - ega/src/operators/meet.rs | 1 - ega/src/operators/mod.rs | 32 --- ega/src/operators/scalar_product.rs | 1 - ega/src/operators/sub.rs | 1 - ega/src/outermorphisms/matrix_4x4.rs | 352 +++++++++++++++++++++++++ ega/src/outermorphisms/mod.rs | 7 + 11 files changed, 398 insertions(+), 39 deletions(-) create mode 100644 ega/src/outermorphisms/matrix_4x4.rs create mode 100644 ega/src/outermorphisms/mod.rs diff --git a/ega/src/lib.rs b/ega/src/lib.rs index 5c91a72..ec6efcb 100644 --- a/ega/src/lib.rs +++ b/ega/src/lib.rs @@ -5,11 +5,43 @@ use ::libm::Libm; mod operators; mod optional_features; +mod outermorphisms; mod values; pub use operators::*; +pub use outermorphisms::*; pub use values::*; +/// return `Empty` +#[inline] +fn return_empty(_: Lhs, _: Rhs) -> Empty { + Empty +} + +/// return the left-hand-side +#[inline] +fn return_lhs(lhs: Lhs, _: Rhs) -> Lhs { + lhs +} + +/// return the right-hand-side +#[inline] +fn return_rhs(_: Lhs, rhs: Rhs) -> Rhs { + rhs +} + +/// negate and return the right-hand-side +#[inline] +fn return_neg_rhs>(_: Lhs, rhs: Rhs) -> Rhs { + -rhs +} + +/// return zero +#[inline] +fn return_zero(_: Lhs, _: Rhs) -> Output { + Output::zero() +} + #[rustfmt::skip] #[cfg(any(test, doctest))] mod test_values { @@ -83,4 +115,11 @@ mod test_values { pub const PSEUDOSCALAR_A: Pseudoscalar = Pseudoscalar { e0123: 397. }; pub const PSEUDOSCALAR_B: Pseudoscalar = Pseudoscalar { e0123: 401. }; pub const PSEUDOSCALAR_C: Pseudoscalar = Pseudoscalar { e0123: -409. }; + + pub const MATRIX_4X4_A: Matrix4x4 = Matrix4x4 { + m00: 1.0, m01: 0.5, m02: 0.5, m03: 0.5, + m10: 0.5, m11: 0.9, m12: -0.5, m13: -1.5, + m20: 0.5, m21: -1.0, m22: 0.5, m23: -2.0, + m30: 0.5, m31: -0.9, m32: 2.0, m33: -1.0, + }; } diff --git a/ega/src/operators/add.rs b/ega/src/operators/add.rs index 03d26df..c38fec3 100644 --- a/ega/src/operators/add.rs +++ b/ega/src/operators/add.rs @@ -1,4 +1,3 @@ -use super::{return_lhs, return_rhs}; use crate::*; pub use ::core::ops::Add; diff --git a/ega/src/operators/dot.rs b/ega/src/operators/dot.rs index d613bb7..de54528 100644 --- a/ega/src/operators/dot.rs +++ b/ega/src/operators/dot.rs @@ -1,4 +1,3 @@ -use super::return_empty; use crate::*; /// The inner product diff --git a/ega/src/operators/geometric_product.rs b/ega/src/operators/geometric_product.rs index ce52503..bf9b4e5 100644 --- a/ega/src/operators/geometric_product.rs +++ b/ega/src/operators/geometric_product.rs @@ -1,4 +1,3 @@ -use super::return_empty; use crate::*; /// The geometric product diff --git a/ega/src/operators/join.rs b/ega/src/operators/join.rs index 3c62ad9..d93cf6e 100644 --- a/ega/src/operators/join.rs +++ b/ega/src/operators/join.rs @@ -1,4 +1,3 @@ -use super::return_empty; use crate::*; /// The regressive product diff --git a/ega/src/operators/meet.rs b/ega/src/operators/meet.rs index e41706d..7356371 100644 --- a/ega/src/operators/meet.rs +++ b/ega/src/operators/meet.rs @@ -1,4 +1,3 @@ -use super::return_empty; use crate::*; /// The outer product diff --git a/ega/src/operators/mod.rs b/ega/src/operators/mod.rs index a396a56..5e3809b 100644 --- a/ega/src/operators/mod.rs +++ b/ega/src/operators/mod.rs @@ -38,8 +38,6 @@ pub use reverse::Reverse; pub use scalar_product::ScalarProduct; pub use sub::Sub; -use crate::{values::Empty, Zero}; - /// Grade Involution pub trait Involution { type Output; @@ -55,33 +53,3 @@ pub trait Exponent { /// Exponentiation fn exp(self, rhs: Rhs) -> Self::Output; } - -/// return `Empty` -#[inline] -fn return_empty(_: Lhs, _: Rhs) -> Empty { - Empty -} - -/// return the left-hand-side -#[inline] -fn return_lhs(lhs: Lhs, _: Rhs) -> Lhs { - lhs -} - -/// return the right-hand-side -#[inline] -fn return_rhs(_: Lhs, rhs: Rhs) -> Rhs { - rhs -} - -/// negate and return the right-hand-side -#[inline] -fn return_neg_rhs>(_: Lhs, rhs: Rhs) -> Rhs { - -rhs -} - -/// return zero -#[inline] -fn return_zero(_: Lhs, _: Rhs) -> Output { - Output::zero() -} diff --git a/ega/src/operators/scalar_product.rs b/ega/src/operators/scalar_product.rs index 76825df..86fe832 100644 --- a/ega/src/operators/scalar_product.rs +++ b/ega/src/operators/scalar_product.rs @@ -1,4 +1,3 @@ -use super::return_zero; use crate::*; /// Scalar Product diff --git a/ega/src/operators/sub.rs b/ega/src/operators/sub.rs index 9e15ecc..ecd9fe5 100644 --- a/ega/src/operators/sub.rs +++ b/ega/src/operators/sub.rs @@ -1,4 +1,3 @@ -use super::{return_lhs, return_neg_rhs}; use crate::*; pub use core::ops::Sub; diff --git a/ega/src/outermorphisms/matrix_4x4.rs b/ega/src/outermorphisms/matrix_4x4.rs new file mode 100644 index 0000000..564ce50 --- /dev/null +++ b/ega/src/outermorphisms/matrix_4x4.rs @@ -0,0 +1,352 @@ +use crate::*; + +/// A matrix whose column vectors have basis [e0, e1, e2, e3] +/// +/// Note: This is different to matrices from other libraries, which often have +/// a basis ordered more like [e1, e2, e3, e0] (homogeneous coordinates) +#[rustfmt::skip] +#[derive(Debug, Copy, Clone)] +pub struct Matrix4x4 { + pub m00: f32, pub m01: f32, pub m02: f32, pub m03: f32, + pub m10: f32, pub m11: f32, pub m12: f32, pub m13: f32, + pub m20: f32, pub m21: f32, pub m22: f32, pub m23: f32, + pub m30: f32, pub m31: f32, pub m32: f32, pub m33: f32, +} + +impl Matrix4x4 { + #[rustfmt::skip] + pub const ZERO: Self = Self { + m00: 0., m01: 0., m02: 0., m03: 0., + m10: 0., m11: 0., m12: 0., m13: 0., + m20: 0., m21: 0., m22: 0., m23: 0., + m30: 0., m31: 0., m32: 0., m33: 0., + }; + + #[rustfmt::skip] + pub const IDENTITY: Self = Self { + m00: 1., m01: 0., m02: 0., m03: 0., + m10: 0., m11: 1., m12: 0., m13: 0., + m20: 0., m21: 0., m22: 1., m23: 0., + m30: 0., m31: 0., m32: 0., m33: 1., + }; +} + +impl Zero for Matrix4x4 { + #[inline] + fn zero() -> Self { + Matrix4x4::ZERO + } +} + +macro_rules! impl_outermorphism { + ($morph_fn:ident: $morph:ty, $arg:ty => $ret:ty) => { + impl Outermorphism<$arg> for $morph { + fn morph(self, arg: $arg) -> $ret { + $morph_fn(self, arg) + } + } + }; +} + +impl_outermorphism! { return_empty: Matrix4x4, Empty => Empty } +impl_outermorphism! { matrix4x4_morph_scalar: Matrix4x4, Scalar => Scalar } +impl_outermorphism! { matrix4x4_morph_vector: Matrix4x4, Vector => Vector } +impl_outermorphism! { matrix4x4_morph_bivector: Matrix4x4, Bivector => Bivector } +impl_outermorphism! { matrix4x4_morph_trivector: Matrix4x4, Trivector => Trivector } +impl_outermorphism! { matrix4x4_morph_pseudo_scalar: Matrix4x4, Pseudoscalar => Pseudoscalar } +impl_outermorphism! { matrix4x4_morph_multivector: Matrix4x4, Multivector => Multivector } + +#[inline] +fn matrix4x4_morph_scalar(_: Matrix4x4, s: Scalar) -> Scalar { + s +} + +#[rustfmt::skip] +#[inline] +fn matrix4x4_morph_vector( + Matrix4x4 { + m00, m01, m02, m03, + m10, m11, m12, m13, + m20, m21, m22, m23, + m30, m31, m32, m33, + }: Matrix4x4, + Vector { + e0: v0, + e1: v1, + e2: v2, + e3: v3 + }: Vector +) -> Vector { + + let e0 = m00*v0 + m01*v1 + m02*v2 + m03*v3; + let e1 = m10*v0 + m11*v1 + m12*v2 + m13*v3; + let e2 = m20*v0 + m21*v1 + m22*v2 + m23*v3; + let e3 = m30*v0 + m31*v1 + m32*v2 + m33*v3; + + Vector { e0, e1, e2, e3 } +} + +#[rustfmt::skip] +#[inline] +fn matrix4x4_morph_bivector( + Matrix4x4 { + m00, m01, m02, m03, + m10, m11, m12, m13, + m20, m21, m22, m23, + m30, m31, m32, m33, + }: Matrix4x4, + Bivector { + e23: b23, e31: b31, e12: b12, + e01: b01, e02: b02, e03: b03, + }: Bivector +) -> Bivector { + + let w1 = m11*b01 + m12*b02 + m13*b03; + let x1 = m10*b01 - m12*b12 + m13*b31; + let y1 = m10*b02 - m13*b23 + m11*b12; + let z1 = m10*b03 - m11*b31 + m12*b23; + + let w2 = m21*b01 + m22*b02 + m23*b03; + let x2 = m20*b01 - m22*b12 + m23*b31; + let y2 = m20*b02 - m23*b23 + m21*b12; + let z2 = m20*b03 - m21*b31 + m22*b23; + + let w3 = m31*b01 + m32*b02 + m33*b03; + let x3 = m30*b01 - m32*b12 + m33*b31; + let y3 = m30*b02 - m33*b23 + m31*b12; + let z3 = m30*b03 - m31*b31 + m32*b23; + + let e23 = m20*w3 - m21*x3 - m22*y3 - m23*z3; + let e03 = m00*w3 - m01*x3 - m02*y3 - m03*z3; + + let e31 = m30*w1 - m31*x1 - m32*y1 - m33*z1; + let e01 = m00*w1 - m01*x1 - m02*y1 - m03*z1; + + let e12 = m10*w2 - m11*x2 - m12*y2 - m13*z2; + let e02 = m00*w2 - m01*x2 - m02*y2 - m03*z2; + + Bivector { e23, e31, e12, e01, e02, e03 } +} + +#[rustfmt::skip] +#[inline] +fn matrix4x4_morph_trivector( + Matrix4x4 { + m00, m01, m02, m03, + m10, m11, m12, m13, + m20, m21, m22, m23, + m30, m31, m32, m33 + }: Matrix4x4, + Trivector { + e123: t123, + e032: t032, + e013: t013, + e021: t021 + }: Trivector +) -> Trivector { + + let a = m22*m33 - m23*m32; + let b = m23*m31 - m21*m33; + let c = m21*m32 - m22*m31; + let d = m22*m30 - m20*m32; + let e = m23*m30 - m20*m33; + let f = m21*m30 - m20*m31; + + let g = m02*m13 - m03*m12; + let h = m03*m11 - m01*m13; + let i = m01*m12 - m02*m11; + let j = m02*m10 - m00*m12; + let k = m03*m10 - m00*m13; + let l = m01*m10 - m00*m11; + + let e123 = + t123*(m12*b + m13*c + m11*a) + - t032*(m10*a + m12*e - m13*d) + - t013*(m10*b + m13*f - m11*e) + - t021*(m10*c + m11*d - m12*f) + ; + let e032 = -( + t123*(m02*b + m03*c + m01*a) + - t032*(m00*a + m02*e - m03*d) + - t013*(m00*b + m03*f - m01*e) + - t021*(m00*c + m01*d - m02*f) + ); + + let e013 = + t123*(m32*h + m33*i + m31*g) + - t032*(m30*g + m32*k - m33*j) + - t013*(m30*h + m33*l - m31*k) + - t021*(m30*i + m31*j - m32*l) + ; + let e021 = -( + t123*(m22*h + m23*i + m21*g) + - t032*(m20*g + m22*k - m23*j) + - t013*(m20*h + m23*l - m21*k) + - t021*(m20*i + m21*j - m22*l) + ); + + Trivector { + e123, e032, e013, e021 + } +} + +#[rustfmt::skip] +#[inline] +fn matrix4x4_morph_pseudo_scalar( + Matrix4x4 { + m00, m01, m02, m03, + m10, m11, m12, m13, + m20, m21, m22, m23, + m30, m31, m32, m33, + }: Matrix4x4, + Pseudoscalar { + e0123: p0123 + }: Pseudoscalar, +) -> Pseudoscalar { + + let e0123 = p0123*( + (m00*m11 - m01*m10)*(m22*m33 - m23*m32) + + (m00*m12 - m02*m10)*(m23*m31 - m21*m33) + + (m00*m13 - m03*m10)*(m21*m32 - m22*m31) + + (m01*m12 - m02*m11)*(m20*m33 - m23*m30) + + (m02*m13 - m03*m12)*(m20*m31 - m21*m30) + + (m03*m11 - m01*m13)*(m20*m32 - m22*m30) + ); + + Pseudoscalar { e0123 } +} + +#[rustfmt::skip] +#[inline] +fn matrix4x4_morph_multivector( + m: Matrix4x4, + Multivector { + e0, e1, e2, e3, + s, e23, e31, e12, + e01, e02, e03, e0123, + e123, e032, e013, e021, + }: Multivector, +) -> Multivector { + let Vector { e0, e1, e2, e3 } = + matrix4x4_morph_vector(m, Vector { e0, e1, e2, e3 }); + let Scalar { s } = matrix4x4_morph_scalar(m, Scalar { s }); + let Bivector { e23, e31, e12, e01, e02, e03 } = matrix4x4_morph_bivector( + m, + Bivector { e23, e31, e12, e01, e02, e03 }, + ); + let Pseudoscalar { e0123 } = + matrix4x4_morph_pseudo_scalar(m, Pseudoscalar { e0123 }); + let Trivector { e123, e032, e013, e021 } = matrix4x4_morph_trivector( + m, + Trivector { e123, e032, e013, e021 } + ); + + Multivector { + e0, e1, e2, e3, + s, e23, e31, e12, + e01, e02, e03, e0123, + e123, e032, e013, e021, + } +} + +#[cfg(any(test, doctest))] +mod tests { + use super::*; + use crate::test_values::*; + use ::approx::assert_ulps_eq; + + #[test] + fn identity_morph_scalar() { + let result = Matrix4x4::IDENTITY.morph(SCALAR_A); + let expected = SCALAR_A; + assert_eq!(dbg!(result), dbg!(expected)); + } + + #[test] + fn identity_morph_vector() { + let result = Matrix4x4::IDENTITY.morph(VECTOR_A); + let expected = VECTOR_A; + assert_eq!(dbg!(result), dbg!(expected)); + } + + #[test] + fn identity_morph_bivector() { + let result = Matrix4x4::IDENTITY.morph(BIVECTOR_A); + let expected = BIVECTOR_A; + assert_eq!(dbg!(result), dbg!(expected)); + } + + #[test] + fn identity_morph_trivector() { + let result = Matrix4x4::IDENTITY.morph(TRIVECTOR_A); + let expected = TRIVECTOR_A; + assert_ulps_eq!(dbg!(result), dbg!(expected), max_ulps = 1); + } + + #[test] + fn identity_morph_pseudoscalar() { + let result = Matrix4x4::IDENTITY.morph(PSEUDOSCALAR_A); + let expected = PSEUDOSCALAR_A; + assert_eq!(dbg!(result), dbg!(expected)); + } + + #[test] + fn identity_morph_multivector() { + let result = Matrix4x4::IDENTITY.morph(MULTIVECTOR_A); + let expected = MULTIVECTOR_A; + assert_eq!(dbg!(result), dbg!(expected)); + } + + #[test] + fn morph_scalar() { + let result = MATRIX_4X4_A.morph(SCALAR_A); + let expected = SCALAR_A; + assert_eq!(dbg!(result), dbg!(expected)); + } + + #[test] + #[rustfmt::skip] + fn morph_vector() { + let result = MATRIX_4X4_A.morph(VECTOR_A); + let expected = Vector { + e0: 394.5, e1: -115.2, e2: -334., e3: 93.2 + }; + assert_ulps_eq!(dbg!(result), dbg!(expected), max_ulps = 2); + } + + #[test] + fn morph_bivector() { + let a = + Meet::meet(MATRIX_4X4_A.morph(VECTOR_A), MATRIX_4X4_A.morph(VECTOR_B)); + let b = MATRIX_4X4_A.morph(Meet::meet(VECTOR_A, VECTOR_B)); + assert_ulps_eq!(dbg!(a), dbg!(b), max_ulps = 90); + } + + #[test] + fn morph_trivector() { + let a = + Meet::meet(MATRIX_4X4_A.morph(VECTOR_A), MATRIX_4X4_A.morph(BIVECTOR_A)); + let b = MATRIX_4X4_A.morph(Meet::meet(VECTOR_A, BIVECTOR_A)); + assert_ulps_eq!(dbg!(a), dbg!(b), max_ulps = 2); + } + + #[test] + fn morph_pseudoscalar() { + let a = Meet::meet( + MATRIX_4X4_A.morph(VECTOR_A), + MATRIX_4X4_A.morph(TRIVECTOR_A), + ); + let b = MATRIX_4X4_A.morph(Meet::meet(VECTOR_A, TRIVECTOR_A)); + assert_ulps_eq!(dbg!(a), dbg!(b), max_ulps = 1); + } + + #[test] + fn morph_multivector() { + let a = Meet::meet( + MATRIX_4X4_A.morph(MULTIVECTOR_A), + MATRIX_4X4_A.morph(MULTIVECTOR_B), + ); + let b = MATRIX_4X4_A.morph(Meet::meet(MULTIVECTOR_A, MULTIVECTOR_B)); + assert_ulps_eq!(dbg!(a), dbg!(b), max_ulps = 10); + } +} diff --git a/ega/src/outermorphisms/mod.rs b/ega/src/outermorphisms/mod.rs new file mode 100644 index 0000000..0737bf9 --- /dev/null +++ b/ega/src/outermorphisms/mod.rs @@ -0,0 +1,7 @@ +mod matrix_4x4; +pub use matrix_4x4::Matrix4x4; + +/// A linear mapping extended to multivectors +pub trait Outermorphism { + fn morph(self, arg: Arg) -> Arg; +}