From 9fca3c9328ffbdbdb853a889362f75f4c38a6d7d Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Mon, 22 Jul 2019 16:50:28 -0400
Subject: [PATCH] eja: get the characteristic_polynomial() for EJAs working.

---
 mjo/eja/euclidean_jordan_algebra.py | 80 +++++++++++++++++++++++------
 1 file changed, 65 insertions(+), 15 deletions(-)

diff --git a/mjo/eja/euclidean_jordan_algebra.py b/mjo/eja/euclidean_jordan_algebra.py
index da1145d..b5cef0a 100644
--- a/mjo/eja/euclidean_jordan_algebra.py
+++ b/mjo/eja/euclidean_jordan_algebra.py
@@ -69,6 +69,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             True
 
         """
+        self._charpoly = None # for caching
         self._rank = rank
         self._natural_basis = natural_basis
         self._multiplication_table = mult_table
@@ -87,28 +88,52 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         return fmt.format(self.degree(), self.base_ring())
 
 
+
     def characteristic_polynomial(self):
+        """
+        EXAMPLES:
+
+        The characteristic polynomial in the spin algebra is given in
+        Alizadeh, Example 11.11::
+
+            sage: J = JordanSpinEJA(3)
+            sage: p = J.characteristic_polynomial(); p
+            X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
+            sage: xvec = J.one().vector()
+            sage: p(*xvec)
+            t^2 - 2*t + 1
+
+        """
+        if self._charpoly is not None:
+            return self._charpoly
+
         r = self.rank()
         n = self.dimension()
 
-        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)
-        x0 = J.zero()
+        # First, compute the basis B...
+        x0 = self.zero()
         c = 1
-        for g in J.gens():
+        for g in self.gens():
             x0 += c*g
             c +=1
         if not x0.is_regular():
             raise ValueError("don't know a regular element")
 
+        V = x0.vector().parent().ambient_vector_space()
+        V1 = V.span_of_basis( (x0**k).vector() for k in range(self.rank()) )
+        B =  (V1.basis() + V1.complement().basis())
+
+        # 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)
+        B = [ b.change_ring(R.fraction_field()) for b in B ]
         # Get the vector space (as opposed to module) so that
         # span_of_basis() works.
-        V = x0.vector().parent().ambient_vector_space()
-        V1 = V.span_of_basis( (x0**k).vector() for k in range(r) )
-        B = V1.basis() + V1.complement().basis()
+        V = J.zero().vector().parent().ambient_vector_space()
         W = V.span_of_basis(B)
 
         def e(k):
@@ -123,16 +148,41 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         A_of_x = block_matrix(1, n, (l1 + l2))
         xr = W.coordinates((x**r).vector())
         a = []
+        denominator = A_of_x.det() # This is constant
         for i in range(n):
-            A_cols = A.columns()
+            A_cols = A_of_x.columns()
             A_cols[i] = xr
-            numerator = column_matrix(A.base_ring(), A_cols).det()
-            denominator = A.det()
+            numerator = column_matrix(A_of_x.base_ring(), A_cols).det()
             ai = numerator/denominator
             a.append(ai)
 
-        # Note: all entries past the rth should be zero.
-        return a
+        # 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.
+        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
+        # for the other powers (to make them sum to x^r).
+        if (r < n):
+            a[r] = 1 # corresponds to x^r
+        else:
+            # When the rank is equal to the dimension, trying to
+            # assign a[r] goes out-of-bounds.
+            a.append(1) # corresponds to x^r
+
+        self._charpoly = sum( a[k]*(t**k) for k in range(len(a)) )
+        return self._charpoly
 
 
     def inner_product(self, x, y):
-- 
2.44.2