From d65df87de905c6ba0912b84ad971bee5f181a30a Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Wed, 28 Oct 2020 11:04:52 -0400
Subject: [PATCH] eja: add two more experimental rank computations.

---
 mjo/eja/eja_algebra.py | 66 +++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 65 insertions(+), 1 deletion(-)

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index fec8a39..34fc0c3 100644
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -794,7 +794,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return  tuple( self.random_element() for idx in range(count) )
 
 
-    def _rank_computation(self):
+    def _rank_computation1(self):
         r"""
         Compute the rank of this algebra using highly suspicious voodoo.
 
@@ -872,6 +872,70 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         # answer anyway; when d=3, we go up to x^2.
         return d-2
 
+    def _rank_computation2(self):
+        r"""
+        Instead of using the dimension of an ideal, find the rank of a
+        matrix containing polynomials.
+        """
+        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()
+
+        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()).row()
+                     for k in range(n) ]
+
+        from sage.matrix.constructor import block_matrix
+        M = block_matrix(n,1,x_powers)
+        return M.rank()
+
+    def _rank_computation3(self):
+        r"""
+        Similar to :meth:`_rank_computation2`, but it stops echelonizing
+        as soon as it hits the first zero row.
+        """
+        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
+        rows = [x_powers[0]]
+        M = matrix(rows)
+        old_rank = 1
+
+        for d in range(1,n):
+            rows = M.rows() + [x_powers[d]]
+            M = matrix(rows)
+            M.echelonize()
+            new_rank = M.rank()
+            if new_rank == old_rank:
+                return new_rank
+            else:
+                old_rank = new_rank
+
     def rank(self):
         """
         Return the rank of this EJA.
-- 
2.43.2