From ef4a0674cdc17613ce65103fc679506e7c7008de Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Sat, 10 Oct 2020 21:47:50 -0400
Subject: [PATCH] mjo/eja/eja_algebra.py: add WIP orthogonal_idempotents()
 method.

---
 mjo/eja/eja_algebra.py | 59 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 59 insertions(+)

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index 8bee729..4b8d466 100644
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -773,6 +773,65 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return (J0, J5, J1)
 
 
+    def orthogonal_idempotents(self):
+        r"""
+        Generate a set of `r` orthogonal idempotents for this algebra,
+        where `r` is its rank.
+
+        This implementation is based on the so-called "central
+        orthogonal idempotents" implemented for (semisimple) centers
+        of SageMath ``FiniteDimensionalAlgebrasWithBasis``. Since all
+        Euclidean Jordan algebas are commutative (and thus equal to
+        their own centers) and semisimple, the method should work more
+        or less as implemented, if it ever worked in the first place.
+        (I don't know the justification for the original implementation.
+        yet).
+
+        How it works: we loop through the algebras generators, looking
+        for their eigenspaces. If there's more than one eigenspace,
+        and if they result in more than one subalgebra, then we split
+        those subalgebras recursively until we get to subalgebras of
+        dimension one (whose idempotent is the unit element). Why does
+        some generator have to produce at least two subalgebras? I
+        dunno. But it seems to work.
+
+        Beware that Koecher defines the "center" of a Jordan algebra to
+        be something else, because the usual definition is stupid in a
+        (necessarily commutative) Jordan algebra.
+        """
+        if self.dimension() == 1:
+            return [self.one()]
+
+        for g in self.gens():
+            eigenpairs = g.operator().matrix().right_eigenspaces()
+            if len(eigenpairs) >= 2:
+                subalgebras = []
+                for eigval, eigspace in eigenpairs:
+                    # Make sub-EJAs from the matrix eigenspaces...
+                    sb = tuple( self.from_vector(b) for b in eigspace.basis() )
+                    try:
+                        # This will fail if e.g. the eigenspace basis
+                        # contains two elements and their product
+                        # isn't a linear combination of the two of
+                        # them (i.e. the generated EJA isn't actually
+                        # two dimensional).
+                        s = FiniteDimensionalEuclideanJordanSubalgebra(self, sb)
+                        subalgebras.append(s)
+                    except:
+                        pass
+                if len(subalgebras) >= 2:
+                    # apply this method recursively.
+                    return tuple( c.superalgebra_element()
+                                  for subalgebra in subalgebras
+                                  for c in subalgebra.orthogonal_idempotents() )
+
+        # If we got here, the algebra didn't decompose, at least not when we looked at
+        # the eigenspaces corresponding only to basis elements of the algebra. The
+        # implementation I stole says that this should work because of Schur's Lemma,
+        # so I personally blame Schur's Lemma if it does not.
+        raise Exception("Schur's Lemma didn't work!")
+
+
     def random_elements(self, count):
         """
         Return ``count`` random elements as a tuple.
-- 
2.43.2