From a7cf81951d0cb51ea40d9362a75385204596df42 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Tue, 30 Jul 2019 11:11:49 -0400
Subject: [PATCH] eja: compute the unit element of the algebra ourselves.

This brings us one step closer to eliminating the dependency
on the FiniteDimensionalAlgebra class.
---
 mjo/eja/eja_algebra.py | 69 +++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 68 insertions(+), 1 deletion(-)

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index fb840ed..362b03a 100644
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -232,7 +232,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         # the first equivalent thing instead so that we can pass in
         # standard coordinates.
         x = J(W(R.gens()))
-        l1 = [matrix.column(W.coordinates((x**k).vector())) for k in range(r)]
+
+        # Handle the zeroth power separately, because computing
+        # the unit element in J is mathematically suspect.
+        x0 = W.coordinates(self.one().vector())
+        l1  = [ matrix.column(x0) ]
+        l1 += [ matrix.column(W.coordinates((x**k).vector()))
+                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.coordinates((x**r).vector())
@@ -378,6 +384,67 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             return self._natural_basis
 
 
+    @cached_method
+    def one(self):
+        """
+        Return the unit element of this algebra.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
+            ....:                                  random_eja)
+
+        EXAMPLES::
+
+            sage: J = RealCartesianProductEJA(5)
+            sage: J.one()
+            e0 + e1 + e2 + e3 + e4
+
+        TESTS::
+
+        The identity element acts like the identity::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: x = J.random_element()
+            sage: J.one()*x == x and x*J.one() == x
+            True
+
+        The matrix of the unit element's operator is the identity::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: actual = J.one().operator().matrix()
+            sage: expected = matrix.identity(J.base_ring(), J.dimension())
+            sage: actual == expected
+            True
+
+        """
+        # We can brute-force compute the matrices of the operators
+        # that correspond to the basis elements of this algebra.
+        # If some linear combination of those basis elements is the
+        # algebra identity, then the same linear combination of
+        # their matrices has to be the identity matrix.
+        #
+        # Of course, matrices aren't vectors in sage, so we have to
+        # appeal to the "long vectors" isometry.
+        oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]
+
+        # Now we use basis linear algebra to find the coefficients,
+        # of the matrices-as-vectors-linear-combination, which should
+        # work for the original algebra basis too.
+        A = matrix.column(self.base_ring(), oper_vecs)
+
+        # We used the isometry on the left-hand side already, but we
+        # still need to do it for the right-hand side. Recall that we
+        # wanted something that summed to the identity matrix.
+        b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) )
+
+        # Now if there's an identity element in the algebra, this should work.
+        coeffs = A.solve_right(b)
+        return self.linear_combination(zip(coeffs,self.gens()))
+
+
     def rank(self):
         """
         Return the rank of this EJA.
-- 
2.43.2