eja: add a random comment.

[sage.d.git] / mjo / eja / eja_element.py

diff --git a/mjo/eja/eja_element.py b/mjo/eja/eja_element.py
index 7c861834723344ea6403f9f5da289af8aa299ae7..a5880c4c4c3fc083726fbf2adcef6ee071e7cd63 100644 (file)
--- a/mjo/eja/eja_element.py
+++ b/mjo/eja/eja_element.py
@@ -5,7 +5,7 @@ from sage.modules.with_basis.indexed_element import IndexedFreeModuleElement
 # TODO: make this unnecessary somehow.
 from sage.misc.lazy_import import lazy_import
 lazy_import('mjo.eja.eja_algebra', 'FiniteDimensionalEuclideanJordanAlgebra')
 # TODO: make this unnecessary somehow.
 from sage.misc.lazy_import import lazy_import
 lazy_import('mjo.eja.eja_algebra', 'FiniteDimensionalEuclideanJordanAlgebra')

-lazy_import('mjo.eja.eja_subalgebra',
+lazy_import('mjo.eja.eja_element_subalgebra',

             'FiniteDimensionalEuclideanJordanElementSubalgebra')
 from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
 from mjo.eja.eja_utils import _mat2vec
             'FiniteDimensionalEuclideanJordanElementSubalgebra')
 from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
 from mjo.eja.eja_utils import _mat2vec


@@ -32,7 +32,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         Return ``self`` raised to the power ``n``.
 
         Jordan algebras are always power-associative; see for
         Return ``self`` raised to the power ``n``.
 
         Jordan algebras are always power-associative; see for

-        example Faraut and Koranyi, Proposition II.1.2 (ii).
+        example Faraut and Korányi, Proposition II.1.2 (ii).

 
         We have to override this because our superclass uses row
         vectors instead of column vectors! We, on the other hand,
 
         We have to override this because our superclass uses row
         vectors instead of column vectors! We, on the other hand,


@@ -94,7 +94,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         SETUP::
 
 
         SETUP::
 

-            sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,

             ....:                                  random_eja)
 
         EXAMPLES::
             ....:                                  random_eja)
 
         EXAMPLES::


@@ -102,7 +102,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: R = PolynomialRing(QQ, 't')
             sage: t = R.gen(0)
             sage: p = t^4 - t^3 + 5*t - 2
             sage: R = PolynomialRing(QQ, 't')
             sage: t = R.gen(0)
             sage: p = t^4 - t^3 + 5*t - 2

-            sage: J = RealCartesianProductEJA(5)
+            sage: J = HadamardEJA(5)

             sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
             True
 
             sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
             True
 


@@ -111,7 +111,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         We should always get back an element of the algebra::
 
             sage: set_random_seed()
         We should always get back an element of the algebra::
 
             sage: set_random_seed()

-            sage: p = PolynomialRing(QQ, 't').random_element()
+            sage: p = PolynomialRing(AA, 't').random_element()

             sage: J = random_eja()
             sage: x = J.random_element()
             sage: x.apply_univariate_polynomial(p) in J
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: x.apply_univariate_polynomial(p) in J


@@ -135,7 +135,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         SETUP::
 
 
         SETUP::
 

-            sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
+            sage: from mjo.eja.eja_algebra import HadamardEJA

 
         EXAMPLES:
 
 
         EXAMPLES:
 


@@ -143,14 +143,14 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         the identity element is `(t-1)` from which it follows that
         the characteristic polynomial should be `(t-1)^3`::
 
         the identity element is `(t-1)` from which it follows that
         the characteristic polynomial should be `(t-1)^3`::
 

-            sage: J = RealCartesianProductEJA(3)
+            sage: J = HadamardEJA(3)

             sage: J.one().characteristic_polynomial()
             t^3 - 3*t^2 + 3*t - 1
 
         Likewise, the characteristic of the zero element in the
         rank-three algebra `R^{n}` should be `t^{3}`::
 
             sage: J.one().characteristic_polynomial()
             t^3 - 3*t^2 + 3*t - 1
 
         Likewise, the characteristic of the zero element in the
         rank-three algebra `R^{n}` should be `t^{3}`::
 

-            sage: J = RealCartesianProductEJA(3)
+            sage: J = HadamardEJA(3)

             sage: J.zero().characteristic_polynomial()
             t^3
 
             sage: J.zero().characteristic_polynomial()
             t^3
 


@@ -160,7 +160,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         to zero on that element::
 
             sage: set_random_seed()
         to zero on that element::
 
             sage: set_random_seed()

-            sage: x = RealCartesianProductEJA(3).random_element()
+            sage: x = HadamardEJA(3).random_element()

             sage: p = x.characteristic_polynomial()
             sage: x.apply_univariate_polynomial(p)
             0
             sage: p = x.characteristic_polynomial()
             sage: x.apply_univariate_polynomial(p)
             0


@@ -168,7 +168,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         The characteristic polynomials of the zero and unit elements
         should be what we think they are in a subalgebra, too::
 
         The characteristic polynomials of the zero and unit elements
         should be what we think they are in a subalgebra, too::
 

-            sage: J = RealCartesianProductEJA(3)
+            sage: J = HadamardEJA(3)

             sage: p1 = J.one().characteristic_polynomial()
             sage: q1 = J.zero().characteristic_polynomial()
             sage: e0,e1,e2 = J.gens()
             sage: p1 = J.one().characteristic_polynomial()
             sage: q1 = J.zero().characteristic_polynomial()
             sage: e0,e1,e2 = J.gens()


@@ -181,7 +181,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             True
 
         """
             True
 
         """

-        p = self.parent().characteristic_polynomial()
+        p = self.parent().characteristic_polynomial_of()

         return p(*self.to_vector())
 
 
         return p(*self.to_vector())
 
 


@@ -344,6 +344,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,

+            ....:                                  TrivialEJA,

             ....:                                  random_eja)
 
         EXAMPLES::
             ....:                                  random_eja)
 
         EXAMPLES::


@@ -362,6 +363,17 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: x.det()
             -1
 
             sage: x.det()
             -1
 

+        The determinant of the sole element in the rank-zero trivial
+        algebra is ``1``, by three paths of reasoning. First, its
+        characteristic polynomial is a constant ``1``, so the constant
+        term in that polynomial is ``1``. Second, the characteristic
+        polynomial evaluated at zero is again ``1``. And finally, the
+        (empty) product of its eigenvalues is likewise just unity::
+
+            sage: J = TrivialEJA()
+            sage: J.zero().det()
+            1
+

         TESTS:
 
         An element is invertible if and only if its determinant is
         TESTS:
 
         An element is invertible if and only if its determinant is


@@ -373,22 +385,28 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             True
 
         Ensure that the determinant is multiplicative on an associative
             True
 
         Ensure that the determinant is multiplicative on an associative

-        subalgebra as in Faraut and Koranyi's Proposition II.2.2::
+        subalgebra as in Faraut and Korányi's Proposition II.2.2::

 
             sage: set_random_seed()
             sage: J = random_eja().random_element().subalgebra_generated_by()
             sage: x,y = J.random_elements(2)
             sage: (x*y).det() == x.det()*y.det()
             True
 
             sage: set_random_seed()
             sage: J = random_eja().random_element().subalgebra_generated_by()
             sage: x,y = J.random_elements(2)
             sage: (x*y).det() == x.det()*y.det()
             True

-

         """
         P = self.parent()
         r = P.rank()
         """
         P = self.parent()
         r = P.rank()

-        p = P._charpoly_coeff(0)
-        # The _charpoly_coeff function already adds the factor of
-        # -1 to ensure that _charpoly_coeff(0) is really what
-        # appears in front of t^{0} in the charpoly. However,
-        # we want (-1)^r times THAT for the determinant.
+
+        if r == 0:
+            # Special case, since we don't get the a0=1
+            # coefficient when the rank of the algebra
+            # is zero.
+            return P.base_ring().one()
+
+        p = P._charpoly_coefficients()[0]
+        # The _charpoly_coeff function already adds the factor of -1
+        # to ensure that _charpoly_coefficients()[0] is really what
+        # appears in front of t^{0} in the charpoly. However, we want
+        # (-1)^r times THAT for the determinant.

         return ((-1)**r)*p(*self.to_vector())
 
 
         return ((-1)**r)*p(*self.to_vector())
 
 


@@ -403,7 +421,8 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         SETUP::
 
 
         SETUP::
 

-            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
+            ....:                                  JordanSpinEJA,

             ....:                                  random_eja)
 
         EXAMPLES:
             ....:                                  random_eja)
 
         EXAMPLES:


@@ -417,14 +436,21 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: while not x.is_invertible():
             ....:     x = J.random_element()
             sage: x_vec = x.to_vector()
             sage: while not x.is_invertible():
             ....:     x = J.random_element()
             sage: x_vec = x.to_vector()

-            sage: x0 = x_vec[0]
+            sage: x0 = x_vec[:1]

             sage: x_bar = x_vec[1:]
             sage: x_bar = x_vec[1:]

-            sage: coeff = ~(x0^2 - x_bar.inner_product(x_bar))
-            sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
-            sage: x_inverse = coeff*inv_vec
+            sage: coeff = x0.inner_product(x0) - x_bar.inner_product(x_bar)
+            sage: x_inverse = x_vec.parent()(x0.list() + (-x_bar).list())
+            sage: if not coeff.is_zero(): x_inverse = x_inverse/coeff

             sage: x.inverse() == J.from_vector(x_inverse)
             True
 
             sage: x.inverse() == J.from_vector(x_inverse)
             True
 

+        Trying to invert a non-invertible element throws an error:
+
+            sage: JordanSpinEJA(3).zero().inverse()
+            Traceback (most recent call last):
+            ...
+            ValueError: element is not invertible
+

         TESTS:
 
         The identity element is its own inverse::
         TESTS:
 
         The identity element is its own inverse::


@@ -450,13 +476,32 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
             True
 
             sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
             True
 

-        The zero element is never invertible::
+        Proposition II.2.3 in Faraut and Korányi says that the inverse
+        of an element is the inverse of its left-multiplication operator
+        applied to the algebra's identity, when that inverse exists::

 
             sage: set_random_seed()
 
             sage: set_random_seed()

-            sage: J = random_eja().zero().inverse()
-            Traceback (most recent call last):
-            ...
-            ValueError: element is not invertible
+            sage: J = random_eja()
+            sage: x = J.random_element()
+            sage: (not x.operator().is_invertible()) or (
+            ....:    x.operator().inverse()(J.one()) == x.inverse() )
+            True
+
+        Proposition II.2.4 in Faraut and Korányi gives a formula for
+        the inverse based on the characteristic polynomial and the
+        Cayley-Hamilton theorem for Euclidean Jordan algebras::
+
+            sage: set_random_seed()
+            sage: J = ComplexHermitianEJA(3)
+            sage: x = J.random_element()
+            sage: while not x.is_invertible():
+            ....:     x = J.random_element()
+            sage: r = J.rank()
+            sage: a = x.characteristic_polynomial().coefficients(sparse=False)
+            sage: expected  = (-1)^(r+1)/x.det()
+            sage: expected *= sum( a[i+1]*x^i for i in range(r) )
+            sage: x.inverse() == expected
+            True

 
         """
         if not self.is_invertible():
 
         """
         if not self.is_invertible():


@@ -515,6 +560,115 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         return not (p(zero) == zero)
 
 
         return not (p(zero) == zero)
 
 

+    def is_primitive_idempotent(self):
+        """
+        Return whether or not this element is a primitive (or minimal)
+        idempotent.
+
+        A primitive idempotent is a non-zero idempotent that is not
+        the sum of two other non-zero idempotents. Remark 2.7.15 in
+        Baes shows that this is what he refers to as a "minimal
+        idempotent."
+
+        An element of a Euclidean Jordan algebra is a minimal idempotent
+        if it :meth:`is_idempotent` and if its Peirce subalgebra
+        corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes,
+        Proposition 2.7.17).
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            ....:                                  RealSymmetricEJA,
+            ....:                                  TrivialEJA,
+            ....:                                  random_eja)
+
+        WARNING::
+
+        This method is sloooooow.
+
+        EXAMPLES:
+
+        The spectral decomposition of a non-regular element should always
+        contain at least one non-minimal idempotent::
+
+            sage: J = RealSymmetricEJA(3)
+            sage: x = sum(J.gens())
+            sage: x.is_regular()
+            False
+            sage: [ c.is_primitive_idempotent()
+            ....:   for (l,c) in x.spectral_decomposition() ]
+            [False, True]
+
+        On the other hand, the spectral decomposition of a regular
+        element should always be in terms of minimal idempotents::
+
+            sage: J = JordanSpinEJA(4)
+            sage: x = sum( i*J.gens()[i] for i in range(len(J.gens())) )
+            sage: x.is_regular()
+            True
+            sage: [ c.is_primitive_idempotent()
+            ....:   for (l,c) in x.spectral_decomposition() ]
+            [True, True]
+
+        TESTS:
+
+        The identity element is minimal only in an EJA of rank one::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
+            True
+
+        A non-idempotent cannot be a minimal idempotent::
+
+            sage: set_random_seed()
+            sage: J = JordanSpinEJA(4)
+            sage: x = J.random_element()
+            sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
+            False
+
+        Proposition 2.7.19 in Baes says that an element is a minimal
+        idempotent if and only if it's idempotent with trace equal to
+        unity::
+
+            sage: set_random_seed()
+            sage: J = JordanSpinEJA(4)
+            sage: x = J.random_element()
+            sage: expected = (x.is_idempotent() and x.trace() == 1)
+            sage: actual = x.is_primitive_idempotent()
+            sage: actual == expected
+            True
+
+        Primitive idempotents must be non-zero::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: J.zero().is_idempotent()
+            True
+            sage: J.zero().is_primitive_idempotent()
+            False
+
+        As a consequence of the fact that primitive idempotents must
+        be non-zero, there are no primitive idempotents in a trivial
+        Euclidean Jordan algebra::
+
+            sage: J = TrivialEJA()
+            sage: J.one().is_idempotent()
+            True
+            sage: J.one().is_primitive_idempotent()
+            False
+
+        """
+        if not self.is_idempotent():
+            return False
+
+        if self.is_zero():
+            return False
+
+        (_,_,J1) = self.parent().peirce_decomposition(self)
+        return (J1.dimension() == 1)
+
+

     def is_nilpotent(self):
         """
         Return whether or not some power of this element is zero.
     def is_nilpotent(self):
         """
         Return whether or not some power of this element is zero.


@@ -542,10 +696,11 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         TESTS:
 
 
         TESTS:
 

-        The identity element is never nilpotent::
+        The identity element is never nilpotent, except in a trivial EJA::

 
             sage: set_random_seed()
 
             sage: set_random_seed()

-            sage: random_eja().one().is_nilpotent()
+            sage: J = random_eja()
+            sage: J.one().is_nilpotent() and not J.is_trivial()

             False
 
         The additive identity is always nilpotent::
             False
 
         The additive identity is always nilpotent::


@@ -589,11 +744,11 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         TESTS:
 
         The zero element should never be regular, unless the parent
         TESTS:
 
         The zero element should never be regular, unless the parent

-        algebra has dimension one::
+        algebra has dimension less than or equal to one::

 
             sage: set_random_seed()
             sage: J = random_eja()
 
             sage: set_random_seed()
             sage: J = random_eja()

-            sage: J.dimension() == 1 or not J.zero().is_regular()
+            sage: J.dimension() <= 1 or not J.zero().is_regular()

             True
 
         The unit element isn't regular unless the algebra happens to
             True
 
         The unit element isn't regular unless the algebra happens to


@@ -601,7 +756,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
             sage: set_random_seed()
             sage: J = random_eja()
 
             sage: set_random_seed()
             sage: J = random_eja()

-            sage: J.dimension() == 1 or not J.one().is_regular()
+            sage: J.dimension() <= 1 or not J.one().is_regular()

             True
 
         """
             True
 
         """


@@ -639,20 +794,24 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
             sage: set_random_seed()
             sage: J = JordanSpinEJA.random_instance()
 
             sage: set_random_seed()
             sage: J = JordanSpinEJA.random_instance()

+            sage: n = J.dimension()

             sage: x = J.random_element()
             sage: x = J.random_element()

-            sage: x == x.coefficient(0)*J.one() or x.degree() == 2
+            sage: x.degree() == min(n,2) or (x == x.coefficient(0)*J.one())

             True
 
         TESTS:
 
             True
 
         TESTS:
 

-        The zero and unit elements are both of degree one::
+        The zero and unit elements are both of degree one in nontrivial
+        algebras::

 
             sage: set_random_seed()
             sage: J = random_eja()
 
             sage: set_random_seed()
             sage: J = random_eja()

-            sage: J.zero().degree()
-            1
-            sage: J.one().degree()
-            1
+            sage: d = J.zero().degree()
+            sage: (J.is_trivial() and d == 0) or d == 1
+            True
+            sage: d = J.one().degree()
+            sage: (J.is_trivial() and d == 0) or d == 1
+            True

 
         Our implementation agrees with the definition::
 
 
         Our implementation agrees with the definition::
 


@@ -696,18 +855,37 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
             ....:                                  RealSymmetricEJA,
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
             ....:                                  RealSymmetricEJA,

+            ....:                                  TrivialEJA,

             ....:                                  random_eja)
 
             ....:                                  random_eja)
 

+        EXAMPLES:
+
+        Keeping in mind that the polynomial ``1`` evaluates the identity
+        element (also the zero element) of the trivial algebra, it is clear
+        that the polynomial ``1`` is the minimal polynomial of the only
+        element in a trivial algebra::
+
+            sage: J = TrivialEJA()
+            sage: J.one().minimal_polynomial()
+            1
+            sage: J.zero().minimal_polynomial()
+            1
+

         TESTS:
 
         The minimal polynomial of the identity and zero elements are
         TESTS:
 
         The minimal polynomial of the identity and zero elements are

-        always the same::
+        always the same, except in trivial algebras where the minimal
+        polynomial of the unit/zero element is ``1``::

 
             sage: set_random_seed()
             sage: J = random_eja()
 
             sage: set_random_seed()
             sage: J = random_eja()

-            sage: J.one().minimal_polynomial()
+            sage: mu = J.one().minimal_polynomial()
+            sage: t = mu.parent().gen()
+            sage: mu + int(J.is_trivial())*(t-2)

             t - 1
             t - 1

-            sage: J.zero().minimal_polynomial()
+            sage: mu = J.zero().minimal_polynomial()
+            sage: t = mu.parent().gen()
+            sage: mu + int(J.is_trivial())*(t-1)

             t
 
         The degree of an element is (by one definition) the degree
             t
 
         The degree of an element is (by one definition) the degree


@@ -725,7 +903,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         two here so that said elements actually exist::
 
             sage: set_random_seed()
         two here so that said elements actually exist::
 
             sage: set_random_seed()

-            sage: n_max = max(2, JordanSpinEJA._max_test_case_size())
+            sage: n_max = max(2, JordanSpinEJA._max_random_instance_size())

             sage: n = ZZ.random_element(2, n_max)
             sage: J = JordanSpinEJA(n)
             sage: y = J.random_element()
             sage: n = ZZ.random_element(2, n_max)
             sage: J = JordanSpinEJA(n)
             sage: y = J.random_element()


@@ -751,11 +929,11 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         and in particular, a re-scaling of the basis::
 
             sage: set_random_seed()
         and in particular, a re-scaling of the basis::
 
             sage: set_random_seed()

-            sage: n_max = RealSymmetricEJA._max_test_case_size()
+            sage: n_max = RealSymmetricEJA._max_random_instance_size()

             sage: n = ZZ.random_element(1, n_max)
             sage: n = ZZ.random_element(1, n_max)

-            sage: J1 = RealSymmetricEJA(n,QQ)
-            sage: J2 = RealSymmetricEJA(n,QQ,False)
-            sage: X = random_matrix(QQ,n)
+            sage: J1 = RealSymmetricEJA(n)
+            sage: J2 = RealSymmetricEJA(n,normalize_basis=False)
+            sage: X = random_matrix(AA,n)

             sage: X = X*X.transpose()
             sage: x1 = J1(X)
             sage: x2 = J2(X)
             sage: X = X*X.transpose()
             sage: x1 = J1(X)
             sage: x2 = J2(X)


@@ -775,7 +953,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
                 # in the "normal" case without us having to think about it.
                 return self.operator().minimal_polynomial()
 
                 # in the "normal" case without us having to think about it.
                 return self.operator().minimal_polynomial()
 

-        A = self.subalgebra_generated_by()
+        A = self.subalgebra_generated_by(orthonormalize_basis=False)

         return A(self).operator().minimal_polynomial()
 
 
         return A(self).operator().minimal_polynomial()
 
 


@@ -830,6 +1008,10 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         """
         B = self.parent().natural_basis()
         W = self.parent().natural_basis_space()
         """
         B = self.parent().natural_basis()
         W = self.parent().natural_basis_space()

+
+        # This is just a manual "from_vector()", but of course
+        # matrix spaces aren't vector spaces in sage, so they
+        # don't have a from_vector() method.

         return W.linear_combination(zip(B,self.to_vector()))
 
 
         return W.linear_combination(zip(B,self.to_vector()))
 
 


@@ -840,14 +1022,14 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,

-            ....:                                  RealCartesianProductEJA)
+            ....:                                  HadamardEJA)

 
         EXAMPLES::
 
 
         EXAMPLES::
 

-            sage: J = RealCartesianProductEJA(2)
+            sage: J = HadamardEJA(2)

             sage: x = sum(J.gens())
             sage: x.norm()
             sage: x = sum(J.gens())
             sage: x.norm()

-            sqrt(2)
+            1.414213562373095?

 
         ::
 
 
         ::
 


@@ -906,16 +1088,18 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: set_random_seed()
             sage: x = JordanSpinEJA.random_instance().random_element()
             sage: x_vec = x.to_vector()
             sage: set_random_seed()
             sage: x = JordanSpinEJA.random_instance().random_element()
             sage: x_vec = x.to_vector()

+            sage: Q = matrix.identity(x.base_ring(), 0)

             sage: n = x_vec.degree()
             sage: n = x_vec.degree()

-            sage: x0 = x_vec[0]
-            sage: x_bar = x_vec[1:]
-            sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
-            sage: B = 2*x0*x_bar.row()
-            sage: C = 2*x0*x_bar.column()
-            sage: D = matrix.identity(QQ, n-1)
-            sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
-            sage: D = D + 2*x_bar.tensor_product(x_bar)
-            sage: Q = matrix.block(2,2,[A,B,C,D])
+            sage: if n > 0:
+            ....:     x0 = x_vec[0]
+            ....:     x_bar = x_vec[1:]
+            ....:     A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
+            ....:     B = 2*x0*x_bar.row()
+            ....:     C = 2*x0*x_bar.column()
+            ....:     D = matrix.identity(x.base_ring(), n-1)
+            ....:     D = (x0^2 - x_bar.inner_product(x_bar))*D
+            ....:     D = D + 2*x_bar.tensor_product(x_bar)
+            ....:     Q = matrix.block(2,2,[A,B,C,D])

             sage: Q == x.quadratic_representation().matrix()
             True
 
             sage: Q == x.quadratic_representation().matrix()
             True
 


@@ -968,10 +1152,10 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: not x.is_invertible() or (
             ....:   x.quadratic_representation(x.inverse())*Qx
             ....:   ==
             sage: not x.is_invertible() or (
             ....:   x.quadratic_representation(x.inverse())*Qx
             ....:   ==

-            ....:   2*x.operator()*Qex - Qx )
+            ....:   2*Lx*Qex - Qx )

             True
 
             True
 

-            sage: 2*x.operator()*Qex - Qx == Lxx
+            sage: 2*Lx*Qex - Qx == Lxx

             True
 
         Property 5:
             True
 
         Property 5:


@@ -1008,12 +1192,92 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
 
 
 
 
 

+    def spectral_decomposition(self):
+        """
+        Return the unique spectral decomposition of this element.
+
+        ALGORITHM:
+
+        Following Faraut and Korányi's Theorem III.1.1, we restrict this
+        element's left-multiplication-by operator to the subalgebra it
+        generates. We then compute the spectral decomposition of that
+        operator, and the spectral projectors we get back must be the
+        left-multiplication-by operators for the idempotents we
+        seek. Thus applying them to the identity element gives us those
+        idempotents.
+
+        Since the eigenvalues are required to be distinct, we take
+        the spectral decomposition of the zero element to be zero
+        times the identity element of the algebra (which is idempotent,
+        obviously).
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import RealSymmetricEJA
+
+        EXAMPLES:
+
+        The spectral decomposition of the identity is ``1`` times itself,
+        and the spectral decomposition of zero is ``0`` times the identity::
+
+            sage: J = RealSymmetricEJA(3)
+            sage: J.one()
+            e0 + e2 + e5
+            sage: J.one().spectral_decomposition()
+            [(1, e0 + e2 + e5)]
+            sage: J.zero().spectral_decomposition()
+            [(0, e0 + e2 + e5)]
+
+        TESTS::
+
+            sage: J = RealSymmetricEJA(4)
+            sage: x = sum(J.gens())
+            sage: sd = x.spectral_decomposition()
+            sage: l0 = sd[0][0]
+            sage: l1 = sd[1][0]
+            sage: c0 = sd[0][1]
+            sage: c1 = sd[1][1]
+            sage: c0.inner_product(c1) == 0
+            True
+            sage: c0.is_idempotent()
+            True
+            sage: c1.is_idempotent()
+            True
+            sage: c0 + c1 == J.one()
+            True
+            sage: l0*c0 + l1*c1 == x
+            True
+
+        The spectral decomposition should work in subalgebras, too::
+
+            sage: J = RealSymmetricEJA(4)
+            sage: (e0, e1, e2, e3, e4, e5, e6, e7, e8, e9) = J.gens()
+            sage: A = 2*e5 - 2*e8
+            sage: (lambda1, c1) = A.spectral_decomposition()[1]
+            sage: (J0, J5, J1) = J.peirce_decomposition(c1)
+            sage: (f0, f1, f2) = J1.gens()
+            sage: f0.spectral_decomposition()
+            [(0, f2), (1, f0)]
+
+        """
+        A = self.subalgebra_generated_by(orthonormalize_basis=True)
+        result = []
+        for (evalue, proj) in A(self).operator().spectral_decomposition():
+            result.append( (evalue, proj(A.one()).superalgebra_element()) )
+        return result

 
 

-    def subalgebra_generated_by(self):
+    def subalgebra_generated_by(self, orthonormalize_basis=False):

         """
         Return the associative subalgebra of the parent EJA generated
         by this element.
 
         """
         Return the associative subalgebra of the parent EJA generated
         by this element.
 

+        Since our parent algebra is unital, we want "subalgebra" to mean
+        "unital subalgebra" as well; thus the subalgebra that an element
+        generates will itself be a Euclidean Jordan algebra after
+        restricting the algebra operations appropriately. This is the
+        subalgebra that Faraut and Korányi work with in section II.2, for
+        example.
+

         SETUP::
 
             sage: from mjo.eja.eja_algebra import random_eja
         SETUP::
 
             sage: from mjo.eja.eja_algebra import random_eja


@@ -1038,17 +1302,18 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: A(x^2) == A(x)*A(x)
             True
 
             sage: A(x^2) == A(x)*A(x)
             True
 

-        The subalgebra generated by the zero element is trivial::
+        By definition, the subalgebra generated by the zero element is
+        the one-dimensional algebra generated by the identity
+        element... unless the original algebra was trivial, in which
+        case the subalgebra is trivial too::

 
             sage: set_random_seed()
             sage: A = random_eja().zero().subalgebra_generated_by()
 
             sage: set_random_seed()
             sage: A = random_eja().zero().subalgebra_generated_by()

-            sage: A
-            Euclidean Jordan algebra of dimension 0 over...
-            sage: A.one()
-            0
+            sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
+            True

 
         """
 
         """

-        return FiniteDimensionalEuclideanJordanElementSubalgebra(self)
+        return FiniteDimensionalEuclideanJordanElementSubalgebra(self, orthonormalize_basis)

 
 
     def subalgebra_idempotent(self):
 
 
     def subalgebra_idempotent(self):


@@ -1060,18 +1325,25 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
             sage: from mjo.eja.eja_algebra import random_eja
 
 
             sage: from mjo.eja.eja_algebra import random_eja
 

-        TESTS::
+        TESTS:
+
+        Ensure that we can find an idempotent in a non-trivial algebra
+        where there are non-nilpotent elements, or that we get the dumb
+        solution in the trivial algebra::

 
             sage: set_random_seed()
             sage: J = random_eja()
             sage: x = J.random_element()
 
             sage: set_random_seed()
             sage: J = random_eja()
             sage: x = J.random_element()

-            sage: while x.is_nilpotent():
+            sage: while x.is_nilpotent() and not J.is_trivial():

             ....:     x = J.random_element()
             sage: c = x.subalgebra_idempotent()
             sage: c^2 == c
             True
 
         """
             ....:     x = J.random_element()
             sage: c = x.subalgebra_idempotent()
             sage: c^2 == c
             True
 
         """

+        if self.parent().is_trivial():
+            return self
+

         if self.is_nilpotent():
             raise ValueError("this only works with non-nilpotent elements!")
 
         if self.is_nilpotent():
             raise ValueError("this only works with non-nilpotent elements!")
 


@@ -1082,7 +1354,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         # will be minimal for some natural number s...
         s = 0
         minimal_dim = J.dimension()
         # will be minimal for some natural number s...
         s = 0
         minimal_dim = J.dimension()

-        for i in xrange(1, minimal_dim):
+        for i in range(1, minimal_dim):

             this_dim = (u**i).operator().matrix().image().dimension()
             if this_dim < minimal_dim:
                 minimal_dim = this_dim
             this_dim = (u**i).operator().matrix().image().dimension()
             if this_dim < minimal_dim:
                 minimal_dim = this_dim


@@ -1111,14 +1383,23 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         """
         Return my trace, the sum of my eigenvalues.
 
         """
         Return my trace, the sum of my eigenvalues.
 

+        In a trivial algebra, however you want to look at it, the trace is
+        an empty sum for which we declare the result to be zero.
+

         SETUP::
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,

-            ....:                                  RealCartesianProductEJA,
+            ....:                                  HadamardEJA,
+            ....:                                  TrivialEJA,

             ....:                                  random_eja)
 
         EXAMPLES::
 
             ....:                                  random_eja)
 
         EXAMPLES::
 

+            sage: J = TrivialEJA()
+            sage: J.zero().trace()
+            0
+
+        ::

             sage: J = JordanSpinEJA(3)
             sage: x = sum(J.gens())
             sage: x.trace()
             sage: J = JordanSpinEJA(3)
             sage: x = sum(J.gens())
             sage: x.trace()


@@ -1126,7 +1407,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         ::
 
 
         ::
 

-            sage: J = RealCartesianProductEJA(5)
+            sage: J = HadamardEJA(5)

             sage: J.one().trace()
             5
 
             sage: J.one().trace()
             5
 


@@ -1142,7 +1423,13 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         """
         P = self.parent()
         r = P.rank()
         """
         P = self.parent()
         r = P.rank()

-        p = P._charpoly_coeff(r-1)
+
+        if r == 0:
+            # Special case for the trivial algebra where
+            # the trace is an empty sum.
+            return P.base_ring().zero()
+
+        p = P._charpoly_coefficients()[r-1]

         # The _charpoly_coeff function already adds the factor of
         # -1 to ensure that _charpoly_coeff(r-1) is really what
         # appears in front of t^{r-1} in the charpoly. However,
         # The _charpoly_coeff function already adds the factor of
         # -1 to ensure that _charpoly_coeff(r-1) is really what
         # appears in front of t^{r-1} in the charpoly. However,


@@ -1160,8 +1447,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         TESTS:
 
 
         TESTS:
 

-        The trace inner product is commutative, bilinear, and satisfies
-        the Jordan axiom:
+        The trace inner product is commutative, bilinear, and associative::

 
             sage: set_random_seed()
             sage: J = random_eja()
 
             sage: set_random_seed()
             sage: J = random_eja()


@@ -1181,7 +1467,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             ....:              a*x.trace_inner_product(z) )
             sage: actual == expected
             True
             ....:              a*x.trace_inner_product(z) )
             sage: actual == expected
             True

-            sage: # jordan axiom
+            sage: # associative

             sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
             True
 
             sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
             True
 


@@ -1199,21 +1485,21 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,

-            ....:                                  RealCartesianProductEJA)
+            ....:                                  HadamardEJA)

 
         EXAMPLES::
 
 
         EXAMPLES::
 

-            sage: J = RealCartesianProductEJA(2)
+            sage: J = HadamardEJA(2)

             sage: x = sum(J.gens())
             sage: x.trace_norm()
             sage: x = sum(J.gens())
             sage: x.trace_norm()

-            sqrt(2)
+            1.414213562373095?

 
         ::
 
             sage: J = JordanSpinEJA(4)
             sage: x = sum(J.gens())
             sage: x.trace_norm()
 
         ::
 
             sage: J = JordanSpinEJA(4)
             sage: x = sum(J.gens())
             sage: x.trace_norm()

-            2*sqrt(2)
+            2.828427124746190?

 
         """
         return self.trace_inner_product(self).sqrt()
 
         """
         return self.trace_inner_product(self).sqrt()