eja: add back the charpoly basis trickery needed for the theory.

diff --git a/mjo/eja/euclidean_jordan_algebra.py b/mjo/eja/euclidean_jordan_algebra.py
index b31322f9b1a2ee23a4451e249efd8734455d4595..2fa4800b09f4f26ca784532e154cde1665ed7b84 100644 (file)
--- a/mjo/eja/euclidean_jordan_algebra.py
+++ b/mjo/eja/euclidean_jordan_algebra.py
@@ -87,6 +87,33 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         return fmt.format(self.degree(), self.base_ring())
 
 
+    def _a_regular_element(self):
+        """
+        Guess a regular element. Needed to compute the basis for our
+        characteristic polynomial coefficients.
+        """
+        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()
+        V = z.vector().parent().ambient_vector_space()
+        V1 = V.span_of_basis( (z**k).vector() for k in range(self.rank()) )
+        b =  (V1.basis() + V1.complement().basis())
+        return V.span_of_basis(b)
+
 
     @cached_method
     def _charpoly_coeff(self, i):
@@ -98,25 +125,38 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         store the trace/determinant (a_{r-1} and a_{0} respectively)
         separate from the entire characteristic polynomial.
         """
-        (A_of_x, x) = self._charpoly_matrix()
+        (A_of_x, x, xr, detA) = self._charpoly_matrix_system()
         R = A_of_x.base_ring()
-        A_cols = A_of_x.columns()
-        A_cols[i] = (x**self.rank()).vector()
-        numerator = column_matrix(A_of_x.base_ring(), A_cols).det()
-        denominator = A_of_x.det()
+        if i >= self.rank():
+            # Guaranteed by theory
+            return R.zero()
+
+        # Danger: the in-place modification is done for performance
+        # reasons (reconstructing a matrix with huge polynomial
+        # entries is slow), but I don't know how cached_method works,
+        # so it's highly possible that we're modifying some global
+        # list variable by reference, here. In other words, you
+        # probably shouldn't call this method twice on the same
+        # algebra, at the same time, in two threads
+        Ai_orig = A_of_x.column(i)
+        A_of_x.set_column(i,xr)
+        numerator = A_of_x.det()
+        A_of_x.set_column(i,Ai_orig)
 
         # We're relying on the theory here to ensure that each a_i is
         # indeed back in R, and the added negative signs are to make
         # the whole charpoly expression sum to zero.
-        return R(-numerator/denominator)
+        return R(-numerator/detA)
 
 
     @cached_method
-    def _charpoly_matrix(self):
+    def _charpoly_matrix_system(self):
         """
         Compute the matrix whose entries A_ij are polynomials in
-        X1,...,XN. This same matrix is used in more than one method and
-        it's not so fast to construct.
+        X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
+        corresponding to `x^r` and the determinent of the matrix A =
+        [A_ij]. In other words, all of the fixed (cachable) data needed
+        to compute the coefficients of the characteristic polynomial.
         """
         r = self.rank()
         n = self.dimension()
@@ -130,11 +170,30 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
         idmat = identity_matrix(J.base_ring(), n)
 
+        W = self._charpoly_basis_space()
+        W = W.change_ring(R.fraction_field())
+
+        # Starting with the standard coordinates x = (X1,X2,...,Xn)
+        # and then converting the entries to W-coordinates allows us
+        # to pass in the standard coordinates to the charpoly and get
+        # back the right answer. Specifically, with x = (X1,X2,...,Xn),
+        # we have
+        #
+        #   W.coordinates(x^2) eval'd at (standard z-coords)
+        #     =
+        #   W-coords of (z^2)
+        #     =
+        #   W-coords of (standard coords of x^2 eval'd at std-coords of z)
+        #
+        # We want the middle equivalent thing in our matrix, but use
+        # the first equivalent thing instead so that we can pass in
+        # standard coordinates.
         x = J(vector(R, R.gens()))
-        l1 = [column_matrix((x**k).vector()) for k in range(r)]
+        l1 = [column_matrix(W.coordinates((x**k).vector())) for k in range(r)]
         l2 = [idmat.column(k-1).column() for k in range(r+1, n+1)]
         A_of_x = block_matrix(R, 1, n, (l1 + l2))
-        return (A_of_x, x)
+        xr = W.coordinates((x**r).vector())
+        return (A_of_x, x, xr, A_of_x.det())
 
 
     @cached_method