eja: start experimenting with single-element Peirce decomposition.

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index 12207b7c5a8e897738ab21a73361883cae03626f..c33352c85d92c11022804ac1d2e6ce57aad7177f 100644 (file)
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -13,6 +13,7 @@ from sage.combinat.free_module import CombinatorialFreeModule
 from sage.matrix.constructor import matrix
 from sage.matrix.matrix_space import MatrixSpace
 from sage.misc.cachefunc import cached_method
+from sage.misc.lazy_import import lazy_import
 from sage.misc.prandom import choice
 from sage.misc.table import table
 from sage.modules.free_module import FreeModule, VectorSpace
@@ -20,6 +21,8 @@ from sage.rings.all import (ZZ, QQ, RR, RLF, CLF,
                             PolynomialRing,
                             QuadraticField)
 from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement
+lazy_import('mjo.eja.eja_subalgebra',
+            'FiniteDimensionalEuclideanJordanSubalgebra')
 from mjo.eja.eja_utils import _mat2vec
 
 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
@@ -631,6 +634,42 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return self.linear_combination(zip(self.gens(), coeffs))
 
 
+    def peirce_decomposition(self, c):
+        """
+        The Peirce decomposition of this algebra relative to the
+        idempotent ``c``.
+
+        In the future, this can be extended to a complete system of
+        orthogonal idempotents.
+        """
+        if not c.is_idempotent():
+            raise ValueError("element is not idempotent: %s" % c)
+
+        # Default these to what they should be if they turn out to be
+        # trivial, because eigenspaces_left() won't return eigenvalues
+        # corresponding to trivial spaces (e.g. it returns only the
+        # 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
+
+        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
+            else:
+                raise ValueError("unexpected eigenvalue: %s" % eigval)
+
+        return (J0, J2, J1)
+
+
     def random_elements(self, count):
         """
         Return ``count`` random elements as a tuple.