Merge pull request #193 from eric-wieser/use-normalInv

Replace manual use of `(~B_c*(1/(B_c*~B_c)[0])` with a call to normalInv

clifford/_multivector.py

@@ -669,37 +669,55 @@ def normal(self) -> 'MultiVector':
         :math:`\frac{M}{|M|}` up to a sign
         """
 
-        return self / np.sqrt(abs(self.mag2()))
+        return self / abs(self)
 
     def leftLaInv(self) -> 'MultiVector':
         """Return left-inverse using a computational linear algebra method
         proposed by Christian Perwass.
         """
         return self._newMV(self.layout.inv_func(self.value))
 
-    def normalInv(self) -> 'MultiVector':
-        r"""The inverse of itself if :math:`M \tilde M = |M|^2`.
-
-        ..math::
+    def _pick_inv(self, fallback):
+        """Internal helper to choose an appropriate inverse method.
 
-            M^{-1} = \tilde M / (M \tilde M)
+        Parameters
+        ----------
+        fallback : bool, optional
+            If `None`, perform no checks on whether normal inv is appropriate.
+            If `True`, fallback to a linalg approach if necessary.
+            If `False`, raise an error if normal inv is not appropriate.
         """
-
         Madjoint = ~self
         MadjointM = (Madjoint * self)
+        if fallback is not None and not MadjointM.isScalar():
+            if fallback:
+                return self.leftLaInv()
+            else:
+                raise ValueError("no inverse exists for this multivector")
 
-        if MadjointM.isScalar() and abs(MadjointM[()]) > cf._eps:
-            # inverse exists
-            return Madjoint / MadjointM[()]
-        else:
+        MadjointM_scalar = MadjointM[()]
+        if fallback is not None and not abs(MadjointM_scalar) > cf._eps:
             raise ValueError("no inverse exists for this multivector")
 
+        return Madjoint / MadjointM_scalar
+
+    def normalInv(self, check=True) -> 'MultiVector':
+        r"""The inverse of itself if :math:`M \tilde M = |M|^2`.
+
+        Parameters
+        ----------
+        check : bool
+            When true, the default, validate that it is appropriate to use this
+            method of inversion.
+
+        ..math::
+
+            M^{-1} = \tilde M / (M \tilde M)
+        """
+        return self._pick_inv(fallback=False if check else None)
+
     def inv(self) -> 'MultiVector':
-        if (self*~self).isScalar():
-            it = self.normalInv()
-        else:
-            it = self.leftLaInv()
-        return it
+        return self._pick_inv(fallback=True)
 
     leftInv = leftLaInv
     rightInv = leftLaInv
@@ -803,9 +821,9 @@ def factorise(self) -> Tuple[List['MultiVector'], numbers.Number]:
         B_c = self/scale
         for ind in B_max_factors[1:]:
             ei = self.layout.blades_list[ind]
-            fi = (ei.lc(B_c)*(~B_c*(1/(B_c*~B_c)[0]))).normal()
+            fi = (ei.lc(B_c)*B_c.normalInv(check=False)).normal()
             factors.append(fi)
-            B_c = B_c * ~fi * (1 / (fi * ~fi)[0])
+            B_c = B_c * fi.normalInv(check=False)
         factors.append(B_c.normal())
         factors.reverse()
         return factors, scale

clifford/test/test_clifford.py

@@ -293,6 +293,25 @@ def test_indexing_blade_tuple(self):
         # three swaps does
         assert mv[1, 2, 3] == -mv[3, 2, 1]
 
+    def test_normalInv(self):
+        layout, blades = Cl(3)
+        e1 = layout.blades['e1']
+        e2 = layout.blades['e2']
+        e3 = layout.blades['e3']
+        assert (2*e1).normalInv() == (0.5*e1)
+
+        with pytest.raises(ValueError):
+            (0*e1).normalInv()  # divide by 0
+
+        with pytest.raises(ValueError):
+            (1 + e1 + e2).normalInv()  # mixed even and odd grades
+
+        # produces garbage, but doesn't crash
+        (1 + e1 + e2).normalInv(check=False)
+
+        # check that not requiring normalInv works fine
+        assert (1 + e1 + e2).inv() == -1 + e1 + e2
+
 
 class TestBasicConformal41:
     def test_metric(self):