some sandwich products

ega/src/operators/mod.rs

@@ -11,12 +11,12 @@ mod meet;
 mod mul;
 mod neg;
 // mod norm;
+mod div;
 mod normalise;
 mod partial_eq;
 mod reverse;
 mod scalar_product;
 mod sub;
-mod div;
 
 pub use add::Add;
 pub use conjugate::Conjugate;

rga/src/operators/inverse.rs

@@ -86,10 +86,10 @@ impl Inverse for DualNumber {
  let DualNumber { s: a, e1234: b } = self;
 
  if a != 0. {
- let q = 1./(a*a);
+ let q = 1./a;
 
- let s = a*q;
- let e1234 = -b*q;
+ let s = q;
+ let e1234 = -b*q*q;
 
  Some(DualNumber { s, e1234 })
  } else {

rga/src/operators/mod.rs

@@ -49,6 +49,8 @@ mod neg;
 mod partial_eq;
 mod sub;
 
+mod sandwich_product;
+
 pub use antidot_product::AntidotProduct;
 pub use antiwedge_product::AntiwedgeProduct;
 pub use dot_product::DotProduct;
@@ -89,6 +91,8 @@ pub use grade_select::{
  GradeSelect0, GradeSelect1, GradeSelect2, GradeSelect3, GradeSelect4,
 };
 
+pub use sandwich_product::SandwichProduct;
+
 pub mod free_functions {
  use crate::{operators::*, values::*};
 

rga/src/operators/sandwich_product.rs

@@ -0,0 +1,194 @@
+use crate::values::*;
+
+pub trait SandwichProduct<Arg> {
+ type Output;
+
+ fn sandwich_product(self, arg: Arg) -> Self::Output;
+}
+
+macro_rules! impl_sandwich_product {
+ ($($operator:ty, $arg:ty => $output:ty: $product_fn:ident;)*) => {
+ $(impl SandwichProduct<$arg> for $operator {
+ type Output = $output;
+
+ #[inline]
+ fn sandwich_product(self, arg: $arg) -> Self::Output {
+ $product_fn(self, arg)
+ }
+ })*
+ };
+}
+
+impl_sandwich_product! {
+ // Bivector, Vector => Vector: bivector_sandwich_vector;
+ Bivector, Bivector => Bivector: bivector_sandwich_bivector;
+ // Bivector, Trivector => Trivector: bivector_sandwich_trivector;
+ Trivector, Vector => Vector: trivector_sandwich_vector;
+ Trivector, Bivector => Bivector: trivector_sandwich_bivector;
+ Trivector, Trivector => Trivector: trivector_sandwich_trivector;
+}
+
+#[rustfmt::skip]
+#[inline]
+fn bivector_sandwich_bivector(
+ Bivector {
+ e41: o41, e42: o42, e43: o43, e23: o23, e31: o31, e12: o12,
+ }: Bivector,
+ Bivector {
+ e41: a41, e42: a42, e43: a43, e23: a23, e31: a31, e12: a12,
+ }: Bivector,
+) -> Bivector {
+
+ let q = 1./(o12*o12 + o23*o23 + o31*o31);
+ let a = o23*o23 - o31*o31 - o12*o12;
+ let b = o31*o31 - o12*o12 - o23*o23;
+ let c = o12*o12 - o23*o23 - o31*o31;
+ let x = o41*o23;
+ let y = o42*o31;
+ let z = o43*o12;
+
+ let e41 = q*(a41*a
+ + 2.*(q*(x + y + z)*(-2.*o23*(o31*a31 + o12*a12) - a23*a)
+ + a23*(x - y - z)
+ + a31*(o41*o31 + o42*o23) + a12*(o41*o12 + o43*o23)
+ + a42*(o23*o31) + a43*(o12*o23)));
+ let e42 = q*(a42*b
+ + 2.*(q*(x + y + z)*(-2.*o31*(o12*a12 + o23*a23) - a31*b)
+ + a31*(y - z - x)
+ + a12*(o42*o12 + o43*o31) + a23*(o42*o23 + o41*o31)
+ + a43*(o31*o12) + a41*(o23*o31)));
+ let e43 = q*(a43*c
+ + 2.*(q*(x + y + z)*(-2.*o12*(o23*a23 + o31*a31) - a12*c)
+ + a12*(z - x - y)
+ + a23*(o43*o23 + o41*o12) + a31*(o43*o31 + o42*o12)
+ + a42*(o31*o12) + a41*(o12*o23)));
+ let e23 = q*(a23*(o23*o23 - o31*o31 - o12*o12)
+ + 2.*o23*(o31*a31 + o12*a12));
+ let e31 = q*(a31*(o31*o31 - o12*o12 - o23*o23)
+ + 2.*o31*(o12*a12 + o23*a23));
+ let e12 = q*(a12*(o12*o12 - o23*o23 - o31*o31)
+ + 2.*o12*(o23*a23 + o31*a31));
+
+ Bivector { e41, e42, e43, e23, e31, e12 }
+}
+
+#[rustfmt::skip]
+#[inline]
+fn trivector_sandwich_vector(
+ Trivector { e423: o423, e431: o431, e412: o412, e321: o321 }: Trivector,
+ Vector { e1: a1, e2: a2, e3: a3, e4: a4 }: Vector,
+) -> Vector {
+
+ let q = -2./o321;
+
+ let e1 = a1;
+ let e2 = a2;
+ let e3 = a3;
+ let e4 = q*(a1*o423 + a2*o431 + a3*o412) - a4;
+
+ Vector { e1, e2, e3, e4 }
+}
+
+#[rustfmt::skip]
+#[inline]
+fn trivector_sandwich_bivector(
+ Trivector { e423: o423, e431: o431, e412: o412, e321: o321 }: Trivector,
+ Bivector {
+ e41: a41, e42: a42, e43: a43, e23: a23, e31: a31, e12: a12,
+ }: Bivector,
+) -> Bivector {
+
+ let q = 2./o321;
+
+ let e41 = q*(o431*a12 - o412*a31) - a41;
+ let e42 = q*(o412*a23 - o423*a12) - a42;
+ let e43 = q*(o423*a31 - o431*a23) - a43;
+ let e23 = a23;
+ let e31 = a31;
+ let e12 = a12;
+
+ Bivector { e41, e42, e43, e23, e31, e12 }
+}
+
+#[rustfmt::skip]
+#[inline]
+fn trivector_sandwich_trivector(
+ Trivector { e423: o423, e431: o431, e412: o412, e321: o321 }: Trivector,
+ Trivector { e423: a423, e431: a431, e412: a412, e321: a321 }: Trivector,
+) -> Trivector {
+
+ let q = 2.*a321/o321;
+
+ let e423 = q*o423 - a423;
+ let e431 = q*o431 - a431;
+ let e412 = q*o412 - a412;
+ let e321 = a321;
+
+ Trivector { e423, e431, e412, e321 }
+}
+
+#[cfg(any(test, doctest))]
+mod tests {
+ use {
+ crate::{free_functions::*, test_values::*, *},
+ ::approx::assert_ulps_eq,
+ };
+
+ #[test]
+ fn bivector_sandwich_bivector() {
+ let bivector_a = grade_2(MULTIVECTOR_A);
+ let bivector_b = grade_2(MULTIVECTOR_B);
+
+ assert_ulps_eq!(
+ dbg!(geometric_product(
+ geometric_product(bivector_a, bivector_b),
+ inverse(bivector_a).unwrap()
+ )),
+ dbg!(bivector_a.sandwich_product(bivector_b)).into(),
+ epsilon = 0.00_000_001
+ );
+ }
+
+ #[test]
+ fn trivector_sandwich_vector() {
+ let trivector = grade_3(MULTIVECTOR_A);
+ let vector = grade_1(MULTIVECTOR_B);
+
+ assert_ulps_eq!(
+ geometric_product(
+ geometric_product(trivector, vector),
+ inverse(trivector).unwrap()
+ ),
+ trivector.sandwich_product(vector).into(),
+ epsilon = 0.00_000_000_000_1
+ );
+ }
+
+ #[test]
+ fn trivector_sandwich_bivector() {
+ let trivector = grade_3(MULTIVECTOR_A);
+ let bivector = grade_2(MULTIVECTOR_B);
+
+ assert_ulps_eq!(
+ geometric_product(
+ geometric_product(trivector, bivector),
+ inverse(trivector).unwrap()
+ ),
+ trivector.sandwich_product(bivector).into()
+ );
+ }
+
+ #[test]
+ fn trivector_sandwich_trivector() {
+ let trivector_a = grade_3(MULTIVECTOR_A);
+ let trivector_b = grade_3(MULTIVECTOR_B);
+
+ assert_ulps_eq!(
+ geometric_product(
+ geometric_product(trivector_a, trivector_b),
+ inverse(trivector_a).unwrap()
+ ),
+ trivector_a.sandwich_product(trivector_b).into()
+ );
+ }
+}