From 103fede29ce3213dab4a088ab8a7839470a9e341 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Sun, 10 Nov 2019 13:39:19 -0500
Subject: [PATCH] eja: fix the Peirce decomposition and add docs/tests for it.

---
 mjo/eja/eja_algebra.py | 103 +++++++++++++++++++++++++++++++++++------
 1 file changed, 90 insertions(+), 13 deletions(-)

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index c33352c..71cdd6f 100644
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -641,6 +641,82 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         In the future, this can be extended to a complete system of
         orthogonal idempotents.
+
+        INPUT:
+
+          - ``c`` -- an idempotent of this algebra.
+
+        OUTPUT:
+
+        A triple (J0, J5, J1) containing two subalgebras and one subspace
+        of this algebra,
+
+          - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
+            corresponding to the eigenvalue zero.
+
+          - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
+            corresponding to the eigenvalue one-half.
+
+          - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
+            corresponding to the eigenvalue one.
+
+        These are the only possible eigenspaces for that operator, and this
+        algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
+        orthogonal, and are subalgebras of this algebra with the appropriate
+        restrictions.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
+
+        EXAMPLES:
+
+        The canonical example comes from the symmetric matrices, which
+        decompose into diagonal and off-diagonal parts::
+
+            sage: J = RealSymmetricEJA(3)
+            sage: C = matrix(QQ, [ [1,0,0],
+            ....:                  [0,1,0],
+            ....:                  [0,0,0] ])
+            sage: c = J(C)
+            sage: J0,J5,J1 = J.peirce_decomposition(c)
+            sage: J0
+            Euclidean Jordan algebra of dimension 1...
+            sage: J5
+            Vector space of degree 6 and dimension 2...
+            sage: J1
+            Euclidean Jordan algebra of dimension 3...
+
+        TESTS:
+
+        Every algebra decomposes trivially with respect to its identity
+        element::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: J0,J5,J1 = J.peirce_decomposition(J.one())
+            sage: J0.dimension() == 0 and J5.dimension() == 0
+            True
+            sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
+            True
+
+        The identity elements in the two subalgebras are the
+        projections onto their respective subspaces of the
+        superalgebra's identity element::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: x = J.random_element()
+            sage: if not J.is_trivial():
+            ....:     while x.is_nilpotent():
+            ....:         x = J.random_element()
+            sage: c = x.subalgebra_idempotent()
+            sage: J0,J5,J1 = J.peirce_decomposition(c)
+            sage: J1(c) == J1.one()
+            True
+            sage: J0(J.one() - c) == J0.one()
+            True
+
         """
         if not c.is_idempotent():
             raise ValueError("element is not idempotent: %s" % c)
@@ -651,23 +727,24 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         # eigenspace corresponding to lambda=1 if you take the
         # decomposition relative to the identity element).
         trivial = FiniteDimensionalEuclideanJordanSubalgebra(self, ())
-        J0 = trivial # eigenvalue zero
-        J2 = trivial # eigenvalue one-half
-        J1 = trivial # eigenvalue one
+        J0 = trivial                          # eigenvalue zero
+        J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
+        J1 = trivial                          # eigenvalue one
 
         for (eigval, eigspace) in c.operator().matrix().left_eigenspaces():
-            gens = tuple( self.from_vector(b) for b in eigspace.basis() )
-            subalg = FiniteDimensionalEuclideanJordanSubalgebra(self, gens)
-            if eigval == 0:
-                J0 = subalg
-            elif eigval == ~(self.base_ring()(2)):
-                J2 = subalg
-            elif eigval == 1:
-                J1 = subalg
+            if eigval == ~(self.base_ring()(2)):
+                J5 = eigspace
             else:
-                raise ValueError("unexpected eigenvalue: %s" % eigval)
+                gens = tuple( self.from_vector(b) for b in eigspace.basis() )
+                subalg = FiniteDimensionalEuclideanJordanSubalgebra(self, gens)
+                if eigval == 0:
+                    J0 = subalg
+                elif eigval == 1:
+                    J1 = subalg
+                else:
+                    raise ValueError("unexpected eigenvalue: %s" % eigval)
 
-        return (J0, J2, J1)
+        return (J0, J5, J1)
 
 
     def random_elements(self, count):
-- 
2.43.2