eja: refactor the charpoly implementation... it's magically faster?

diff --git a/mjo/eja/euclidean_jordan_algebra.py b/mjo/eja/euclidean_jordan_algebra.py
index 4fb5f9b4a22b9e0184ca6d2e4640717321c1efc2..44ec225b35e4468dd34661983e186e2c429652dc 100644 (file)
--- a/mjo/eja/euclidean_jordan_algebra.py
+++ b/mjo/eja/euclidean_jordan_algebra.py
@@ -88,6 +88,55 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
 
 
+    @cached_method
+    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.
+        """
+        (A_of_x, x) = self._charpoly_matrix()
+        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()
+
+        # 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)
+
+
+    @cached_method
+    def _charpoly_matrix(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.
+        """
+        r = self.rank()
+        n = self.dimension()
+
+        # Construct a new algebra over a multivariate polynomial ring...
+        names = ['X' + str(i) for i in range(1,n+1)]
+        R = PolynomialRing(self.base_ring(), names)
+        J = FiniteDimensionalEuclideanJordanAlgebra(R,
+                                                    self._multiplication_table,
+                                                    rank=r)
+
+        idmat = identity_matrix(J.base_ring(), n)
+
+        x = J(vector(R, R.gens()))
+        l1 = [column_matrix((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)
+
+
     @cached_method
     def characteristic_polynomial(self):
         """
@@ -107,47 +156,19 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         r = self.rank()
         n = self.dimension()
 
-        # Now switch to the polynomial rings.
-        names = ['X' + str(i) for i in range(1,n+1)]
-        R = PolynomialRing(self.base_ring(), names)
-        J = FiniteDimensionalEuclideanJordanAlgebra(R,
-                                                    self._multiplication_table,
-                                                    rank=r)
-
-        def e(k):
-            # The coordinates of e_k with respect to the basis B.
-            # But, the e_k are elements of B...
-            return identity_matrix(J.base_ring(), n).column(k-1).column()
-
-        # A matrix implementation 1
-        x = J(vector(R, R.gens()))
-        l1 = [column_matrix((x**k).vector()) for k in range(r)]
-        l2 = [e(k) for k in range(r+1, n+1)]
-        A_of_x = block_matrix(1, n, (l1 + l2))
-        xr = (x**r).vector()
-        a = []
-        denominator = A_of_x.det() # This is constant
-        for i in range(n):
-            A_cols = A_of_x.columns()
-            A_cols[i] = xr
-            numerator = column_matrix(A_of_x.base_ring(), A_cols).det()
-            ai = numerator/denominator
-            a.append(ai)
+        # The list of coefficient polynomials a_1, a_2, ..., a_n.
+        a = [ self._charpoly_coeff(i) for i in range(n) ]
 
         # We go to a bit of trouble here to reorder the
         # indeterminates, so that it's easier to evaluate the
         # characteristic polynomial at x's coordinates and get back
         # something in terms of t, which is what we want.
+        R = a[0].parent()
         S = PolynomialRing(self.base_ring(),'t')
         t = S.gen(0)
         S = PolynomialRing(S, R.variable_names())
         t = S(t)
 
-        # We're relying on the theory here to ensure that each entry
-        # a[i] is indeed back in R, and the added negative signs are
-        # to make the whole expression sum to zero.
-        a = [R(-ai) for ai in a] # corresponds to powerx x^0 through x^(r-1)
-
         # Note: all entries past the rth should be zero. The
         # coefficient of the highest power (x^r) is 1, but it doesn't
         # appear in the solution vector which contains coefficients