From 119af99174a49d71d1b11fdfd1b0de451fee733f Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Mon, 22 Jul 2019 01:55:51 -0400
Subject: [PATCH] eja: begin to stub out a characteristic_polynomial() for
 EJAs.

---
 mjo/eja/euclidean_jordan_algebra.py | 48 +++++++++++++++++++++++++++++
 1 file changed, 48 insertions(+)

diff --git a/mjo/eja/euclidean_jordan_algebra.py b/mjo/eja/euclidean_jordan_algebra.py
index 2b4a886..29248fc 100644
--- a/mjo/eja/euclidean_jordan_algebra.py
+++ b/mjo/eja/euclidean_jordan_algebra.py
@@ -87,6 +87,54 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         return fmt.format(self.degree(), self.base_ring())
 
 
+    def characteristic_polynomial(self):
+        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()
+        c = 1
+        for g in J.gens():
+            x0 += c*g
+            c +=1
+        if not x0.is_regular():
+            raise ValueError("don't know a regular element")
+
+        # 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()
+        W = V.span_of_basis(B)
+
+        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(W.coordinates((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 = W.coordinates((x**r).vector())
+        a = []
+        for i in range(n):
+            A_cols = A.columns()
+            A_cols[i] = xr
+            numerator = column_matrix(A.base_ring(), A_cols).det()
+            denominator = A.det()
+            ai = numerator/denominator
+            a.append(ai)
+
+        # Note: all entries past the rth should be zero.
+        return a
+
+
     def inner_product(self, x, y):
         """
         The inner product associated with this Euclidean Jordan algebra.
-- 
2.43.2