From: Michael Orlitzky <michael@orlitzky.com>
Date: Wed, 4 Nov 2020 13:01:11 +0000 (-0500)
Subject: eja: replace rank computation with length of charpoly coefficients.
X-Git-Url: http://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=417406c6a509462be431961a3dd1f13829b7cd96;p=sage.d.git

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 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):