impl Inverse for Bivector

ega/src/operators/geometric_product.rs

@@ -523,17 +523,17 @@ fn scalar_mul_vector(
 #[rustfmt::skip]
 #[inline]
 fn scalar_mul_bivector(
- lhs: Scalar,
- mut rhs: Bivector,
+ Scalar { s }: Scalar,
+ Bivector { e23, e31, e12, e01, e02, e03 }: Bivector,
 ) -> Bivector {
- rhs.e01 *= lhs.s;
- rhs.e02 *= lhs.s;
- rhs.e03 *= lhs.s;
- rhs.e12 *= lhs.s;
- rhs.e31 *= lhs.s;
- rhs.e23 *= lhs.s;
-
- rhs
+ let e23 = s*e23;
+ let e31 = s*e31;
+ let e12 = s*e12;
+ let e01 = s*e01;
+ let e02 = s*e02;
+ let e03 = s*e03;
+
+ Bivector { e23, e31, e12, e01, e02, e03 }
 }
 
 #[rustfmt::skip]

ega/src/operators/inverse.rs

@@ -1,51 +1,103 @@
 use crate::*;
 
 pub trait Inverse {
- fn inverse(self) -> Self;
+ type Output;
+
+ fn inverse(self) -> Self::Output;
 }
 
 impl Inverse for Scalar {
+ type Output = Self;
+
  #[inline]
  fn inverse(self) -> Self {
  Scalar { s: 1.0 / self.s }
  }
 }
 
 impl Inverse for Vector {
+ type Output = Self;
+
  #[inline]
  fn inverse(self) -> Self {
  simple_inverse(self)
  }
 }
 
+impl Inverse for Bivector {
+ type Output = Self;
+
+ #[inline]
+ fn inverse(self) -> Self {
+ bivector_inverse(self)
+ }
+}
+
 impl Inverse for Trivector {
+ type Output = Self;
+
  #[inline]
  fn inverse(self) -> Self {
  simple_inverse(self)
  }
 }
 
-// this is only valid for some kinds of elements
+impl Inverse for Pseudoscalar {
+ type Output = Option<Multivector>;
+
+ /// The inverse for a Pseudoscalar does not exist, as there is no pair of
+ /// pseudoscalar & multivector, P & M, that can satisfy the geometric product
+ /// PM = 1
+ #[inline]
+ fn inverse(self) -> Option<Multivector> {
+ None
+ }
+}
+
+// note: this is only valid for some kinds of "simple" elements
+// note: this is valid for bivectors only when 0 = e01 = e02 = e03
 #[inline]
 fn simple_inverse<T: Copy + Reverse + NormSquared + Mul<f32, Output = T>>(
  value: T,
 ) -> T {
  value.reverse() * (1.0 / value.norm_squared().s)
 }
 
+// works for general bivectors
+#[rustfmt::skip]
+#[inline]
+fn bivector_inverse(
+ Bivector { e23: be23, e31: be31, e12: be12, e01: be01, e02: be02, e03: be03 }: Bivector
+) -> Bivector {
+ let x = -(be23*be23 + be31*be31 + be12*be12);
+ let y = 2.*(be03*be12 + be02*be31 + be01*be23);
+ let q = 1./(x*x);
+
+ let oe23 = q*be23*x;
+ let oe31 = q*be31*x;
+ let oe12 = q*be12*x;
+ let oe01 = q*(be01*x + be23*y);
+ let oe02 = q*(be02*x + be31*y);
+ let oe03 = q*(be03*x + be12*y);
+
+ Bivector {
+ e23: oe23, e31: oe31, e12: oe12, e01: oe01, e02: oe02, e03: oe03
+ }
+}
+
 #[cfg(any(test, doctest))]
 mod tests {
  use super::*;
  use crate::test_values::*;
- use ::approx::assert_relative_eq;
+ use ::approx::assert_ulps_eq;
 
  #[test]
  fn inverse_scalar_a() {
  let inverse = SCALAR_A.inverse();
  let product = SCALAR_A * inverse;
 
  let expected = Scalar::UNIT;
- assert_relative_eq!(dbg!(expected), dbg!(product));
+ assert_ulps_eq!(dbg!(expected), dbg!(product), max_ulps = 0);
  }
 
  #[test]
@@ -54,7 +106,7 @@ mod tests {
  let product = SCALAR_B * inverse;
 
  let expected = Scalar::UNIT;
- assert_relative_eq!(dbg!(expected), dbg!(product));
+ assert_ulps_eq!(dbg!(expected), dbg!(product), max_ulps = 0);
  }
 
  #[test]
@@ -63,7 +115,7 @@ mod tests {
  let product = SCALAR_C * inverse;
 
  let expected = Scalar::UNIT;
- assert_relative_eq!(dbg!(expected), dbg!(product));
+ assert_ulps_eq!(dbg!(expected), dbg!(product), max_ulps = 0);
  }
 
  #[test]
@@ -72,7 +124,7 @@ mod tests {
  let product = VECTOR_A * inverse;
 
  let expected = Multivector::from(Scalar::UNIT);
- assert_relative_eq!(dbg!(expected), dbg!(product));
+ assert_ulps_eq!(dbg!(expected), dbg!(product), max_ulps = 0);
  }
 
  #[test]
@@ -81,7 +133,7 @@ mod tests {
  let product = VECTOR_B * inverse;
 
  let expected = Multivector::from(Scalar::UNIT);
- assert_relative_eq!(dbg!(expected), dbg!(product));
+ assert_ulps_eq!(dbg!(expected), dbg!(product), max_ulps = 0);
  }
 
  #[test]
@@ -90,7 +142,34 @@ mod tests {
  let product = VECTOR_C * inverse;
 
  let expected = Multivector::from(Scalar::UNIT);
- assert_relative_eq!(dbg!(expected), dbg!(product));
+ assert_ulps_eq!(dbg!(expected), dbg!(product), max_ulps = 0);
+ }
+
+ #[test]
+ fn inverse_bivector_a() {
+ let inverse = BIVECTOR_A.inverse();
+ let product = BIVECTOR_A * inverse;
+
+ let expected = Multivector::from(Scalar::UNIT);
+ assert_ulps_eq!(dbg!(expected), dbg!(product), max_ulps = 0);
+ }
+
+ #[test]
+ fn inverse_bivector_b() {
+ let inverse = BIVECTOR_B.inverse();
+ let product = BIVECTOR_B * inverse;
+
+ let expected = Multivector::from(Scalar::UNIT);
+ assert_ulps_eq!(dbg!(expected), dbg!(product), max_ulps = 0);
+ }
+
+ #[test]
+ fn inverse_bivector_c() {
+ let inverse = BIVECTOR_C.inverse();
+ let product = BIVECTOR_C * inverse;
+
+ let expected = Multivector::from(Scalar::UNIT);
+ assert_ulps_eq!(dbg!(expected), dbg!(product), max_ulps = 0);
  }
 
  #[test]
@@ -99,7 +178,7 @@ mod tests {
  let product = TRIVECTOR_A * inverse;
 
  let expected = Multivector::from(Scalar::UNIT);
- assert_relative_eq!(dbg!(expected), dbg!(product));
+ assert_ulps_eq!(dbg!(expected), dbg!(product), max_ulps = 0);
  }
 
  #[test]
@@ -108,7 +187,7 @@ mod tests {
  let product = TRIVECTOR_B * inverse;
 
  let expected = Multivector::from(Scalar::UNIT);
- assert_relative_eq!(dbg!(expected), dbg!(product));
+ assert_ulps_eq!(dbg!(expected), dbg!(product), max_ulps = 0);
  }
 
  #[test]
@@ -117,6 +196,6 @@ mod tests {
  let product = TRIVECTOR_C * inverse;
 
  let expected = Multivector::from(Scalar::UNIT);
- assert_relative_eq!(dbg!(expected), dbg!(product));
+ assert_ulps_eq!(dbg!(expected), dbg!(product), max_ulps = 0);
  }
 }

ega/src/operators/norm.rs

@@ -70,7 +70,7 @@ impl NormSquared for Vector {
 impl NormSquared for Bivector {
  #[inline]
  fn norm_squared(self) -> Scalar {
- let s = self.e12 * self.e12 + self.e31 * self.e31 + self.e23 * self.e23;
+ let s = self.e23 * self.e23 + self.e31 * self.e31 + self.e12 * self.e12;
 
  Scalar { s }
  }