eja: speed up a few slow examples.

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

diff --git a/mjo/eja/eja_element.py b/mjo/eja/eja_element.py
index 4d3b071923b2c275ecb559c8d60d183696c729ac..d99a7d873ab69e511186b7f0bc3ba3127a45c85a 100644 (file)
--- a/mjo/eja/eja_element.py
+++ b/mjo/eja/eja_element.py
@@ -1,16 +1,7 @@
-# -*- coding: utf-8 -*-
-
-from itertools import izip
-

 from sage.matrix.constructor import matrix
 from sage.modules.free_module import VectorSpace
 from sage.modules.with_basis.indexed_element import IndexedFreeModuleElement
 
 from sage.matrix.constructor import matrix
 from sage.modules.free_module import VectorSpace
 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')
-lazy_import('mjo.eja.eja_element_subalgebra',
-            'FiniteDimensionalEuclideanJordanElementSubalgebra')

 from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
 from mjo.eja.eja_utils import _mat2vec
 
 from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
 from mjo.eja.eja_utils import _mat2vec
 


@@ -98,7 +89,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::


@@ -106,7 +97,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
 


@@ -115,7 +106,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


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

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

 
         EXAMPLES:
 
 
         EXAMPLES:
 


@@ -147,14 +138,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
 


@@ -164,7 +155,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


@@ -172,7 +163,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()


@@ -185,7 +176,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())
 
 


@@ -236,9 +227,9 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         Ditto for the quaternions::
 
 
         Ditto for the quaternions::
 

-            sage: J = QuaternionHermitianEJA(3)
+            sage: J = QuaternionHermitianEJA(2)

             sage: J.one().inner_product(J.one())
             sage: J.one().inner_product(J.one())

-            3
+            2

 
         TESTS:
 
 
         TESTS:
 


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

+            ....:                                  TrivialEJA,
+            ....:                                  RealSymmetricEJA,
+            ....:                                  ComplexHermitianEJA,

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


@@ -366,6 +360,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


@@ -385,14 +390,51 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: (x*y).det() == x.det()*y.det()
             True
 
             sage: (x*y).det() == x.det()*y.det()
             True
 

+        The determinant in matrix algebras is just the usual determinant::
+
+            sage: set_random_seed()
+            sage: X = matrix.random(QQ,3)
+            sage: X = X + X.T
+            sage: J1 = RealSymmetricEJA(3)
+            sage: J2 = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
+            sage: expected = X.det()
+            sage: actual1 = J1(X).det()
+            sage: actual2 = J2(X).det()
+            sage: actual1 == expected
+            True
+            sage: actual2 == expected
+            True
+
+        ::
+
+            sage: set_random_seed()
+            sage: J1 = ComplexHermitianEJA(2)
+            sage: J2 = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: X = matrix.random(GaussianIntegers(), 2)
+            sage: X = X + X.H
+            sage: expected = AA(X.det())
+            sage: actual1 = J1(J1.real_embed(X)).det()
+            sage: actual2 = J2(J2.real_embed(X)).det()
+            sage: expected == actual1
+            True
+            sage: expected == actual2
+            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())
 
 


@@ -422,11 +464,11 @@ 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
 


@@ -493,6 +535,14 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         if not self.is_invertible():
             raise ValueError("element is not invertible")
 
         if not self.is_invertible():
             raise ValueError("element is not invertible")
 

+        if self.parent()._charpoly_coefficients.is_in_cache():
+            # We can invert using our charpoly if it will be fast to
+            # compute. If the coefficients are cached, our rank had
+            # better be too!
+            r = self.parent().rank()
+            a = self.characteristic_polynomial().coefficients(sparse=False)
+            return (-1)**(r+1)*sum(a[i+1]*self**i for i in range(r))/self.det()
+

         return (~self.quadratic_representation())(self)
 
 
         return (~self.quadratic_representation())(self)
 
 


@@ -506,9 +556,14 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         zero, but we need the characteristic polynomial for the
         determinant. The minimal polynomial is a lot easier to get,
         so we use Corollary 2 in Chapter V of Koecher to check
         zero, but we need the characteristic polynomial for the
         determinant. The minimal polynomial is a lot easier to get,
         so we use Corollary 2 in Chapter V of Koecher to check

-        whether or not the paren't algebra's zero element is a root
+        whether or not the parent algebra's zero element is a root

         of this element's minimal polynomial.
 
         of this element's minimal polynomial.
 

+        That is... unless the coefficients of our algebra's
+        "characteristic polynomial of" function are already cached!
+        In that case, we just use the determinant (which will be fast
+        as a result).
+

         Beware that we can't use the superclass method, because it
         relies on the algebra being associative.
 
         Beware that we can't use the superclass method, because it
         relies on the algebra being associative.
 


@@ -539,6 +594,11 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             else:
                 return False
 
             else:
                 return False
 

+        if self.parent()._charpoly_coefficients.is_in_cache():
+            # The determinant will be quicker than computing the minimal
+            # polynomial from scratch, most likely.
+            return (not self.det().is_zero())
+

         # In fact, we only need to know if the constant term is non-zero,
         # so we can pass in the field's zero element instead.
         zero = self.base_ring().zero()
         # In fact, we only need to know if the constant term is non-zero,
         # so we can pass in the field's zero element instead.
         zero = self.base_ring().zero()


@@ -546,10 +606,15 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         return not (p(zero) == zero)
 
 
         return not (p(zero) == zero)
 
 

-    def is_minimal_idempotent(self):
+    def is_primitive_idempotent(self):

         """
         """

-        Return whether or not this element is a minimal idempotent.
+        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
 
         An element of a Euclidean Jordan algebra is a minimal idempotent
         if it :meth:`is_idempotent` and if its Peirce subalgebra


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

+            ....:                                  TrivialEJA,

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


@@ -571,22 +637,22 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         The spectral decomposition of a non-regular element should always
         contain at least one non-minimal idempotent::
 
         The spectral decomposition of a non-regular element should always
         contain at least one non-minimal idempotent::
 

-            sage: J = RealSymmetricEJA(3, AA)
+            sage: J = RealSymmetricEJA(3)

             sage: x = sum(J.gens())
             sage: x.is_regular()
             False
             sage: x = sum(J.gens())
             sage: x.is_regular()
             False

-            sage: [ c.is_minimal_idempotent()
+            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::
 
             ....:   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, AA)
+            sage: J = JordanSpinEJA(4)

             sage: x = sum( i*J.gens()[i] for i in range(len(J.gens())) )
             sage: x.is_regular()
             True
             sage: x = sum( i*J.gens()[i] for i in range(len(J.gens())) )
             sage: x.is_regular()
             True

-            sage: [ c.is_minimal_idempotent()
+            sage: [ c.is_primitive_idempotent()

             ....:   for (l,c) in x.spectral_decomposition() ]
             [True, True]
 
             ....:   for (l,c) in x.spectral_decomposition() ]
             [True, True]
 


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

-            sage: J.rank() == 1 or not J.one().is_minimal_idempotent()
+            sage: J.rank() == 1 or not J.one().is_primitive_idempotent()

             True
 
         A non-idempotent cannot be a minimal idempotent::
             True
 
         A non-idempotent cannot be a minimal idempotent::


@@ -604,7 +670,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: set_random_seed()
             sage: J = JordanSpinEJA(4)
             sage: x = J.random_element()
             sage: set_random_seed()
             sage: J = JordanSpinEJA(4)
             sage: x = J.random_element()

-            sage: (not x.is_idempotent()) and x.is_minimal_idempotent()
+            sage: (not x.is_idempotent()) and x.is_primitive_idempotent()

             False
 
         Proposition 2.7.19 in Baes says that an element is a minimal
             False
 
         Proposition 2.7.19 in Baes says that an element is a minimal


@@ -615,14 +681,36 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: J = JordanSpinEJA(4)
             sage: x = J.random_element()
             sage: expected = (x.is_idempotent() and x.trace() == 1)
             sage: J = JordanSpinEJA(4)
             sage: x = J.random_element()
             sage: expected = (x.is_idempotent() and x.trace() == 1)

-            sage: actual = x.is_minimal_idempotent()
+            sage: actual = x.is_primitive_idempotent()

             sage: actual == expected
             True
 
             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 not self.is_idempotent():
             return False
 

+        if self.is_zero():
+            return False
+

         (_,_,J1) = self.parent().peirce_decomposition(self)
         return (J1.dimension() == 1)
 
         (_,_,J1) = self.parent().peirce_decomposition(self)
         return (J1.dimension() == 1)
 


@@ -752,8 +840,9 @@ 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:


@@ -831,13 +920,18 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         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: set_random_seed()

-            sage: J = random_eja(nontrivial=True)
-            sage: J.one().minimal_polynomial()
+            sage: J = random_eja()
+            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


@@ -855,7 +949,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()


@@ -881,11 +975,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,normalize_basis=False)
-            sage: X = random_matrix(QQ,n)
+            sage: J1 = RealSymmetricEJA(n)
+            sage: J2 = RealSymmetricEJA(n,orthonormalize=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)


@@ -905,20 +999,21 @@ 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()
 
 
 

-    def natural_representation(self):
+    def to_matrix(self):

         """
         """

-        Return a more-natural representation of this element.
+        Return an (often more natural) representation of this element as a
+        matrix.

 
 

-        Every finite-dimensional Euclidean Jordan Algebra is a
-        direct sum of five simple algebras, four of which comprise
-        Hermitian matrices. This method returns the original
-        "natural" representation of this element as a Hermitian
-        matrix, if it has one. If not, you get the usual representation.
+        Every finite-dimensional Euclidean Jordan Algebra is a direct
+        sum of five simple algebras, four of which comprise Hermitian
+        matrices. This method returns a "natural" matrix
+        representation of this element as either a Hermitian matrix or
+        column vector.

 
         SETUP::
 
 
         SETUP::
 


@@ -930,7 +1025,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: J = ComplexHermitianEJA(3)
             sage: J.one()
             e0 + e3 + e8
             sage: J = ComplexHermitianEJA(3)
             sage: J.one()
             e0 + e3 + e8

-            sage: J.one().natural_representation()
+            sage: J.one().to_matrix()

             [1 0 0 0 0 0]
             [0 1 0 0 0 0]
             [0 0 1 0 0 0]
             [1 0 0 0 0 0]
             [0 1 0 0 0 0]
             [0 0 1 0 0 0]


@@ -940,27 +1035,27 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         ::
 
 
         ::
 

-            sage: J = QuaternionHermitianEJA(3)
+            sage: J = QuaternionHermitianEJA(2)

             sage: J.one()
             sage: J.one()

-            e0 + e5 + e14
-            sage: J.one().natural_representation()
-            [1 0 0 0 0 0 0 0 0 0 0 0]
-            [0 1 0 0 0 0 0 0 0 0 0 0]
-            [0 0 1 0 0 0 0 0 0 0 0 0]
-            [0 0 0 1 0 0 0 0 0 0 0 0]
-            [0 0 0 0 1 0 0 0 0 0 0 0]
-            [0 0 0 0 0 1 0 0 0 0 0 0]
-            [0 0 0 0 0 0 1 0 0 0 0 0]
-            [0 0 0 0 0 0 0 1 0 0 0 0]
-            [0 0 0 0 0 0 0 0 1 0 0 0]
-            [0 0 0 0 0 0 0 0 0 1 0 0]
-            [0 0 0 0 0 0 0 0 0 0 1 0]
-            [0 0 0 0 0 0 0 0 0 0 0 1]
+            e0 + e5
+            sage: J.one().to_matrix()
+            [1 0 0 0 0 0 0 0]
+            [0 1 0 0 0 0 0 0]
+            [0 0 1 0 0 0 0 0]
+            [0 0 0 1 0 0 0 0]
+            [0 0 0 0 1 0 0 0]
+            [0 0 0 0 0 1 0 0]
+            [0 0 0 0 0 0 1 0]
+            [0 0 0 0 0 0 0 1]

 
         """
 
         """

-        B = self.parent().natural_basis()
-        W = self.parent().natural_basis_space()
-        return W.linear_combination(izip(B,self.to_vector()))
+        B = self.parent().matrix_basis()
+        W = self.parent().matrix_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()) )

 
 
     def norm(self):
 
 
     def norm(self):


@@ -970,14 +1065,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?

 
         ::
 
 
         ::
 


@@ -1036,16 +1131,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
 


@@ -1166,7 +1263,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         The spectral decomposition of the identity is ``1`` times itself,
         and the spectral decomposition of zero is ``0`` times the identity::
 
         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,AA)
+            sage: J = RealSymmetricEJA(3)

             sage: J.one()
             e0 + e2 + e5
             sage: J.one().spectral_decomposition()
             sage: J.one()
             e0 + e2 + e5
             sage: J.one().spectral_decomposition()


@@ -1176,7 +1273,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         TESTS::
 
 
         TESTS::
 

-            sage: J = RealSymmetricEJA(4,AA)
+            sage: J = RealSymmetricEJA(4)

             sage: x = sum(J.gens())
             sage: sd = x.spectral_decomposition()
             sage: l0 = sd[0][0]
             sage: x = sum(J.gens())
             sage: sd = x.spectral_decomposition()
             sage: l0 = sd[0][0]


@@ -1194,8 +1291,18 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: l0*c0 + l1*c1 == x
             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)]
+

         """
         """

-        P = self.parent()

         A = self.subalgebra_generated_by(orthonormalize_basis=True)
         result = []
         for (evalue, proj) in A(self).operator().spectral_decomposition():
         A = self.subalgebra_generated_by(orthonormalize_basis=True)
         result = []
         for (evalue, proj) in A(self).operator().spectral_decomposition():


@@ -1249,6 +1356,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             True
 
         """
             True
 
         """

+        from mjo.eja.eja_element_subalgebra import FiniteDimensionalEuclideanJordanElementSubalgebra

         return FiniteDimensionalEuclideanJordanElementSubalgebra(self, orthonormalize_basis)
 
 
         return FiniteDimensionalEuclideanJordanElementSubalgebra(self, orthonormalize_basis)
 
 


@@ -1264,12 +1372,13 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         TESTS:
 
         Ensure that we can find an idempotent in a non-trivial algebra
         TESTS:
 
         Ensure that we can find an idempotent in a non-trivial algebra

-        where there are non-nilpotent elements::
+        where there are non-nilpotent elements, or that we get the dumb
+        solution in the trivial algebra::

 
             sage: set_random_seed()
 
             sage: set_random_seed()

-            sage: J = random_eja(nontrivial=True)
+            sage: J = random_eja()

             sage: x = J.random_element()
             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
             ....:     x = J.random_element()
             sage: c = x.subalgebra_idempotent()
             sage: c^2 == c


@@ -1289,7 +1398,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


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

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

             ....:                                  TrivialEJA,
             ....:                                  random_eja)
 
             ....:                                  TrivialEJA,
             ....:                                  random_eja)
 


@@ -1342,7 +1451,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         ::
 
 
         ::
 

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

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


@@ -1364,7 +1473,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             # the trace is an empty sum.
             return P.base_ring().zero()
 
             # the trace is an empty sum.
             return P.base_ring().zero()
 

-        p = P._charpoly_coeff(r-1)
+        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,


@@ -1420,21 +1529,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()