eja: replace rank computation with length of charpoly coefficients.

If we cache the solution to the charpoly system instead, then the
result is useful in more than one place (for the charpoly itself,
in particular).

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index 01deb1b65140675ffbc81c9c673c21e6d4a2e77f..d1c0fa2eafd5d3d5e17ccab175e16ee44c3e9d82 100644 (file)
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -826,8 +826,45 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return cls(n, field, **kwargs)
 
     @cached_method
-    def rank(self):
+    def _charpoly_coefficients(self):
+        r"""
+        The `r` polynomial coefficients of the "characteristic polynomial
+        of" function.
         """
+        n = self.dimension()
+        var_names = [ "X" + str(z) for z in range(1,n+1) ]
+        R = PolynomialRing(self.base_ring(), var_names)
+        vars = R.gens()
+        F = R.fraction_field()
+
+        def L_x_i_j(i,j):
+            # From a result in my book, these are the entries of the
+            # basis representation of L_x.
+            return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
+                        for k in range(n) )
+
+        L_x = matrix(F, n, n, L_x_i_j)
+        # Compute an extra power in case the rank is equal to
+        # the dimension (otherwise, we would stop at x^(r-1)).
+        x_powers = [ (L_x**k)*self.one().to_vector()
+                     for k in range(n+1) ]
+        A = matrix.column(F, x_powers[:n])
+        AE = A.extended_echelon_form()
+        E = AE[:,n:]
+        A_rref = AE[:,:n]
+        r = A_rref.rank()
+        b = x_powers[r]
+
+        # The theory says that only the first "r" coefficients are
+        # nonzero, and they actually live in the original polynomial
+        # ring and not the fraction field. We negate them because
+        # in the actual characteristic polynomial, they get moved
+        # to the other side where x^r lives.
+        return -A_rref.solve_right(E*b).change_ring(R)[:r]
+
+    @cached_method
+    def rank(self):
+        r"""
         Return the rank of this EJA.
 
         ALGORITHM:
@@ -920,43 +957,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             2
 
         """
-        n = self.dimension()
-        if n == 0:
-            return 0
-        elif n == 1:
-            return 1
-
-        var_names = [ "X" + str(z) for z in range(1,n+1) ]
-        R = PolynomialRing(self.base_ring(), var_names)
-        vars = R.gens()
-
-        def L_x_i_j(i,j):
-            # From a result in my book, these are the entries of the
-            # basis representation of L_x.
-            return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
-                        for k in range(n) )
-
-        L_x = matrix(R, n, n, L_x_i_j)
-        x_powers = [ vars[k]*(L_x**k)*self.one().to_vector()
-                     for k in range(n) ]
-
-        # Can assume n >= 2
-        M = matrix([x_powers[0]])
-        old_rank = 1
-
-        for d in range(1,n):
-            M = matrix(M.rows() + [x_powers[d]])
-            M.echelonize()
-            # TODO: we've basically solved the system here.
-            # We should save the echelonized matrix somehow
-            # so that it can be reused in the charpoly method.
-            new_rank = M.rank()
-            if new_rank == old_rank:
-                return new_rank
-            else:
-                old_rank = new_rank
-
-        return n
+        return len(self._charpoly_coefficients())
 
 
     def vector_space(self):