eja: fix the fast matrix _charpoly_coefficients() method.

[sage.d.git] / mjo / eja / eja_algebra.py

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index 2b769ac447bce5b149782bbf96ada23631f9b2c9..4091d03fd6a40351e7d4c2789847d1db6fa0e807 100644 (file)
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -235,66 +235,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
     def product_on_basis(self, i, j):
         return self._multiplication_table[i][j]
 
-    def _a_regular_element(self):
-        """
-        Guess a regular element. Needed to compute the basis for our
-        characteristic polynomial coefficients.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import random_eja
-
-        TESTS:
-
-        Ensure that this hacky method succeeds for every algebra that we
-        know how to construct::
-
-            sage: set_random_seed()
-            sage: J = random_eja()
-            sage: J._a_regular_element().is_regular()
-            True
-
-        """
-        gs = self.gens()
-        z = self.sum( (i+1)*gs[i] for i in range(len(gs)) )
-        if not z.is_regular():
-            raise ValueError("don't know a regular element")
-        return z
-
-
-    @cached_method
-    def _charpoly_basis_space(self):
-        """
-        Return the vector space spanned by the basis used in our
-        characteristic polynomial coefficients. This is used not only to
-        compute those coefficients, but also any time we need to
-        evaluate the coefficients (like when we compute the trace or
-        determinant).
-        """
-        z = self._a_regular_element()
-        # Don't use the parent vector space directly here in case this
-        # happens to be a subalgebra. In that case, we would be e.g.
-        # two-dimensional but span_of_basis() would expect three
-        # coordinates.
-        V = VectorSpace(self.base_ring(), self.vector_space().dimension())
-        basis = [ (z**k).to_vector() for k in range(self.rank()) ]
-        V1 = V.span_of_basis( basis )
-        b =  (V1.basis() + V1.complement().basis())
-        return V.span_of_basis(b)
-
-
-    def _charpoly_coeff(self, i):
-        """
-        Return the coefficient polynomial "a_{i}" of this algebra's
-        general characteristic polynomial.
-
-        Having this be a separate cached method lets us compute and
-        store the trace/determinant (a_{r-1} and a_{0} respectively)
-        separate from the entire characteristic polynomial.
-        """
-        return self._charpoly_coefficients()[i]
-
-
     @cached_method
     def characteristic_polynomial(self):
         """
@@ -748,6 +688,14 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
                         for k in range(n) )
 
         L_x = matrix(F, n, n, L_x_i_j)
+
+        r = None
+        if self.rank.is_in_cache():
+            r = self.rank()
+            # There's no need to pad the system with redundant
+            # columns if we *know* they'll be redundant.
+            n = r
+
         # 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()
@@ -756,7 +704,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         AE = A.extended_echelon_form()
         E = AE[:,n:]
         A_rref = AE[:,:n]
-        r = A_rref.rank()
+        if r is None:
+            r = A_rref.rank()
         b = x_powers[r]
 
         # The theory says that only the first "r" coefficients are
@@ -1044,7 +993,24 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
             J = MatrixEuclideanJordanAlgebra(QQ,
                                              basis,
                                              normalize_basis=False)
-            return J._charpoly_coefficients()
+            a = J._charpoly_coefficients()
+
+            # Unfortunately, changing the basis does change the
+            # coefficients of the characteristic polynomial, but since
+            # these are really the coefficients of the "characteristic
+            # polynomial of" function, everything is still nice and
+            # unevaluated. It's therefore "obvious" how scaling the
+            # basis affects the coordinate variables X1, X2, et
+            # cetera. Scaling the first basis vector up by "n" adds a
+            # factor of 1/n into every "X1" term, for example. So here
+            # we simply undo the basis_normalizer scaling that we
+            # performed earlier.
+            #
+            # TODO: make this access safe.
+            XS = a[0].variables()
+            subs_dict = { XS[i]: self._basis_normalizers[i]*XS[i]
+                          for i in range(len(XS)) }
+            return tuple( a_i.subs(subs_dict) for a_i in a )
 
 
     @staticmethod