From: Michael Orlitzky <michael@orlitzky.com>
Date: Wed, 7 Aug 2019 00:09:09 +0000 (-0400)
Subject: eja: rewrite the hacky process used for characteristic polynomials.
X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=13a49fd59d57cf34b026460ace004ddbd60d4b68;p=sage.d.git

eja: rewrite the hacky process used for characteristic polynomials.
---

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index a9bb6c7..3855f0e 100644
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -259,16 +259,34 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         r = self.rank()
         n = self.dimension()
 
-        # Construct a new algebra over a multivariate polynomial ring...
+        # Turn my vector space into a module so that "vectors" can
+        # have multivatiate polynomial entries.
         names = tuple('X' + str(i) for i in range(1,n+1))
         R = PolynomialRing(self.base_ring(), names)
-        # Hack around the fact that our multiplication table is in terms of
-        # algebra elements but the constructor wants it in terms of vectors.
-        vmt = [ tuple(map(lambda x: x.to_vector(), ls))
-                for ls in self._multiplication_table ]
-        J = FiniteDimensionalEuclideanJordanAlgebra(R, tuple(vmt), r)
-
-        idmat = matrix.identity(J.base_ring(), n)
+        V = self.vector_space().change_ring(R)
+
+        # Now let x = (X1,X2,...,Xn) be the vector whose entries are
+        # indeterminates...
+        x = V(names)
+
+        # And figure out the "left multiplication by x" matrix in
+        # that setting.
+        lmbx_cols = []
+        monomial_matrices = [ self.monomial(i).operator().matrix()
+                              for i in range(n) ] # don't recompute these!
+        for k in range(n):
+            ek = self.monomial(k).to_vector()
+            lmbx_cols.append(
+              sum( x[i]*(monomial_matrices[i]*ek)
+                   for i in range(n) ) )
+        Lx = matrix.column(R, lmbx_cols)
+
+        # Now we can compute powers of x "symbolically"
+        x_powers = [self.one().to_vector(), x]
+        for d in range(2, r+1):
+            x_powers.append( Lx*(x_powers[-1]) )
+
+        idmat = matrix.identity(R, n)
 
         W = self._charpoly_basis_space()
         W = W.change_ring(R.fraction_field())
@@ -288,18 +306,10 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         # 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.from_vector(W(R.gens()))
-
-        # Handle the zeroth power separately, because computing
-        # the unit element in J is mathematically suspect.
-        x0 = W.coordinate_vector(self.one().to_vector())
-        l1  = [ x0.column() ]
-        l1 += [ W.coordinate_vector((x**k).to_vector()).column()
-                for k in range(1,r) ]
-        l2 = [idmat.column(k-1).column() for k in range(r+1, n+1)]
-        A_of_x = matrix.block(R, 1, n, (l1 + l2))
-        xr = W.coordinate_vector((x**r).to_vector())
-        return (A_of_x, x, xr, A_of_x.det())
+        x_powers = [ W.coordinate_vector(xp) for xp in x_powers ]
+        l2 = [idmat.column(k-1) for k in range(r+1, n+1)]
+        A_of_x = matrix.column(R, n, (x_powers[:r] + l2))
+        return (A_of_x, x, x_powers[r], A_of_x.det())
 
 
     @cached_method