X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=d1c0fa2eafd5d3d5e17ccab175e16ee44c3e9d82;hb=417406c6a509462be431961a3dd1f13829b7cd96;hp=01deb1b65140675ffbc81c9c673c21e6d4a2e77f;hpb=8e4c84f17025f6c3f2d049902a6a01f954d484ee;p=sage.d.git

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