From 417406c6a509462be431961a3dd1f13829b7cd96 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Wed, 4 Nov 2020 08:01:11 -0500
Subject: [PATCH] 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).
---
 mjo/eja/eja_algebra.py | 77 +++++++++++++++++++++---------------------
 1 file changed, 39 insertions(+), 38 deletions(-)

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index 01deb1b..d1c0fa2 100644
--- 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):
-- 
2.44.2