]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_element.py
eja: factor out MatrixEJA initialization.
[sage.d.git] / mjo / eja / eja_element.py
index 1b5e3149e8e804f410403c066cc7e483eaa42a61..c98d9a2b64651562928a6489ef88e03fdd2b0dc9 100644 (file)
@@ -1,14 +1,12 @@
-from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_element import FiniteDimensionalAlgebraElement
 from sage.matrix.constructor import matrix
 from sage.matrix.constructor import matrix
+from sage.misc.cachefunc import cached_method
 from sage.modules.free_module import VectorSpace
 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')
-from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
-from mjo.eja.eja_utils import _mat2vec
+from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
+from mjo.eja.eja_utils import _mat2vec, _scale
 
 
-class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraElement):
+class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
     """
     An element of a Euclidean Jordan algebra.
     """
     """
     An element of a Euclidean Jordan algebra.
     """
@@ -23,75 +21,14 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
                       dir(self.__class__) )
 
 
                       dir(self.__class__) )
 
 
-    def __init__(self, A, elt=None):
-        """
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
-            ....:                                  random_eja)
-
-        EXAMPLES:
-
-        The identity in `S^n` is converted to the identity in the EJA::
-
-            sage: J = RealSymmetricEJA(3)
-            sage: I = matrix.identity(QQ,3)
-            sage: J(I) == J.one()
-            True
-
-        This skew-symmetric matrix can't be represented in the EJA::
-
-            sage: J = RealSymmetricEJA(3)
-            sage: A = matrix(QQ,3, lambda i,j: i-j)
-            sage: J(A)
-            Traceback (most recent call last):
-            ...
-            ArithmeticError: vector is not in free module
-
-        TESTS:
 
 
-        Ensure that we can convert any element of the parent's
-        underlying vector space back into an algebra element whose
-        vector representation is what we started with::
-
-            sage: set_random_seed()
-            sage: J = random_eja()
-            sage: v = J.vector_space().random_element()
-            sage: J(v).vector() == v
-            True
-
-        """
-        # Goal: if we're given a matrix, and if it lives in our
-        # parent algebra's "natural ambient space," convert it
-        # into an algebra element.
-        #
-        # The catch is, we make a recursive call after converting
-        # the given matrix into a vector that lives in the algebra.
-        # This we need to try the parent class initializer first,
-        # to avoid recursing forever if we're given something that
-        # already fits into the algebra, but also happens to live
-        # in the parent's "natural ambient space" (this happens with
-        # vectors in R^n).
-        try:
-            FiniteDimensionalAlgebraElement.__init__(self, A, elt)
-        except ValueError:
-            natural_basis = A.natural_basis()
-            if elt in natural_basis[0].matrix_space():
-                # Thanks for nothing! Matrix spaces aren't vector
-                # spaces in Sage, so we have to figure out its
-                # natural-basis coordinates ourselves.
-                V = VectorSpace(elt.base_ring(), elt.nrows()**2)
-                W = V.span( _mat2vec(s) for s in natural_basis )
-                coords =  W.coordinates(_mat2vec(elt))
-                FiniteDimensionalAlgebraElement.__init__(self, A, coords)
 
     def __pow__(self, n):
         """
         Return ``self`` raised to the power ``n``.
 
         Jordan algebras are always power-associative; see for
 
     def __pow__(self, n):
         """
         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,
@@ -137,7 +74,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
         elif n == 1:
             return self
         else:
         elif n == 1:
             return self
         else:
-            return (self.operator()**(n-1))(self)
+            return (self**(n-1))*self
 
 
     def apply_univariate_polynomial(self, p):
 
 
     def apply_univariate_polynomial(self, p):
@@ -153,7 +90,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
 
         SETUP::
 
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
             ....:                                  random_eja)
 
         EXAMPLES::
             ....:                                  random_eja)
 
         EXAMPLES::
@@ -161,7 +98,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
             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
 
@@ -170,7 +107,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
         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
@@ -194,7 +131,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
 
         SETUP::
 
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
+            sage: from mjo.eja.eja_algebra import HadamardEJA
 
         EXAMPLES:
 
 
         EXAMPLES:
 
@@ -202,14 +139,14 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
         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
 
@@ -219,14 +156,29 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
         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
 
+        The characteristic polynomials of the zero and unit elements
+        should be what we think they are in a subalgebra, too::
+
+            sage: J = HadamardEJA(3)
+            sage: p1 = J.one().characteristic_polynomial()
+            sage: q1 = J.zero().characteristic_polynomial()
+            sage: b0,b1,b2 = J.gens()
+            sage: A = (b0 + 2*b1 + 3*b2).subalgebra_generated_by() # dim 3
+            sage: p2 = A.one().characteristic_polynomial()
+            sage: q2 = A.zero().characteristic_polynomial()
+            sage: p1 == p2
+            True
+            sage: q1 == q2
+            True
+
         """
         """
-        p = self.parent().characteristic_polynomial()
-        return p(*self.vector())
+        p = self.parent().characteristic_polynomial_of()
+        return p(*self.to_vector())
 
 
     def inner_product(self, other):
 
 
     def inner_product(self, other):
@@ -253,7 +205,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
             sage: y = vector(QQ,[4,5,6])
             sage: x.inner_product(y)
             32
             sage: y = vector(QQ,[4,5,6])
             sage: x.inner_product(y)
             32
-            sage: J(x).inner_product(J(y))
+            sage: J.from_vector(x).inner_product(J.from_vector(y))
             32
 
         The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
             32
 
         The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
@@ -276,9 +228,9 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
 
         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:
 
@@ -287,9 +239,8 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
 
             sage: set_random_seed()
             sage: J = random_eja()
 
             sage: set_random_seed()
             sage: J = random_eja()
-            sage: x = J.random_element()
-            sage: y = J.random_element()
-            sage: x.inner_product(y) in RR
+            sage: x,y = J.random_elements(2)
+            sage: x.inner_product(y) in RLF
             True
 
         """
             True
 
         """
@@ -324,9 +275,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
         Test Lemma 1 from Chapter III of Koecher::
 
             sage: set_random_seed()
         Test Lemma 1 from Chapter III of Koecher::
 
             sage: set_random_seed()
-            sage: J = random_eja()
-            sage: u = J.random_element()
-            sage: v = J.random_element()
+            sage: u,v = random_eja().random_elements(2)
             sage: lhs = u.operator_commutes_with(u*v)
             sage: rhs = v.operator_commutes_with(u^2)
             sage: lhs == rhs
             sage: lhs = u.operator_commutes_with(u*v)
             sage: rhs = v.operator_commutes_with(u^2)
             sage: lhs == rhs
@@ -336,9 +285,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
         Chapter III, or from Baes (2.3)::
 
             sage: set_random_seed()
         Chapter III, or from Baes (2.3)::
 
             sage: set_random_seed()
-            sage: J = random_eja()
-            sage: x = J.random_element()
-            sage: y = J.random_element()
+            sage: x,y = random_eja().random_elements(2)
             sage: Lx = x.operator()
             sage: Ly = y.operator()
             sage: Lxx = (x*x).operator()
             sage: Lx = x.operator()
             sage: Ly = y.operator()
             sage: Lxx = (x*x).operator()
@@ -350,10 +297,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
         Baes (2.4)::
 
             sage: set_random_seed()
         Baes (2.4)::
 
             sage: set_random_seed()
-            sage: J = random_eja()
-            sage: x = J.random_element()
-            sage: y = J.random_element()
-            sage: z = J.random_element()
+            sage: x,y,z = random_eja().random_elements(3)
             sage: Lx = x.operator()
             sage: Ly = y.operator()
             sage: Lz = z.operator()
             sage: Lx = x.operator()
             sage: Ly = y.operator()
             sage: Lz = z.operator()
@@ -367,10 +311,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
         Baes (2.5)::
 
             sage: set_random_seed()
         Baes (2.5)::
 
             sage: set_random_seed()
-            sage: J = random_eja()
-            sage: u = J.random_element()
-            sage: y = J.random_element()
-            sage: z = J.random_element()
+            sage: u,y,z = random_eja().random_elements(3)
             sage: Lu = u.operator()
             sage: Ly = y.operator()
             sage: Lz = z.operator()
             sage: Lu = u.operator()
             sage: Ly = y.operator()
             sage: Lz = z.operator()
@@ -399,12 +340,14 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            ....:                                  TrivialEJA,
+            ....:                                  RealSymmetricEJA,
+            ....:                                  ComplexHermitianEJA,
             ....:                                  random_eja)
 
         EXAMPLES::
 
             sage: J = JordanSpinEJA(2)
             ....:                                  random_eja)
 
         EXAMPLES::
 
             sage: J = JordanSpinEJA(2)
-            sage: e0,e1 = J.gens()
             sage: x = sum( J.gens() )
             sage: x.det()
             0
             sage: x = sum( J.gens() )
             sage: x.det()
             0
@@ -412,11 +355,21 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
         ::
 
             sage: J = JordanSpinEJA(3)
         ::
 
             sage: J = JordanSpinEJA(3)
-            sage: e0,e1,e2 = J.gens()
             sage: x = sum( J.gens() )
             sage: x.det()
             -1
 
             sage: x = sum( J.gens() )
             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
@@ -427,29 +380,67 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
             sage: x.is_invertible() == (x.det() != 0)
             True
 
             sage: x.is_invertible() == (x.det() != 0)
             True
 
+        Ensure that the determinant is multiplicative on an associative
+        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
+
+        The determinant in real matrix algebras is 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
+
         """
         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.
-        return ((-1)**r)*p(*self.vector())
+
+        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())
 
 
 
 
+    @cached_method
     def inverse(self):
         """
         Return the Jordan-multiplicative inverse of this element.
 
         ALGORITHM:
 
     def inverse(self):
         """
         Return the Jordan-multiplicative inverse of this element.
 
         ALGORITHM:
 
-        We appeal to the quadratic representation as in Koecher's
-        Theorem 12 in Chapter III, Section 5.
+        In general we appeal to the quadratic representation as in
+        Koecher's Theorem 12 in Chapter III, Section 5. But if the
+        parent algebra's "characteristic polynomial of" coefficients
+        happen to be cached, then we use Proposition II.2.4 in Faraut
+        and Korányi which gives a formula for the inverse based on the
+        characteristic polynomial and the Cayley-Hamilton theorem for
+        Euclidean Jordan algebras::
 
         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:
@@ -458,20 +449,26 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
         Example 11.11::
 
             sage: set_random_seed()
         Example 11.11::
 
             sage: set_random_seed()
-            sage: n = ZZ.random_element(1,10)
-            sage: J = JordanSpinEJA(n)
+            sage: J = JordanSpinEJA.random_instance()
             sage: x = J.random_element()
             sage: while not x.is_invertible():
             ....:     x = J.random_element()
             sage: x = J.random_element()
             sage: while not x.is_invertible():
             ....:     x = J.random_element()
-            sage: x_vec = x.vector()
-            sage: x0 = x_vec[0]
+            sage: x_vec = x.to_vector()
+            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: x.inverse() == J(x_inverse)
+            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
 
             True
 
+        Trying to invert a non-invertible element throws an error:
+
+            sage: JordanSpinEJA(3).zero().inverse()
+            Traceback (most recent call last):
+            ...
+            ZeroDivisionError: element is not invertible
+
         TESTS:
 
         The identity element is its own inverse::
         TESTS:
 
         The identity element is its own inverse::
@@ -497,36 +494,63 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
             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
 
 
+        Check that the fast (cached) and slow algorithms give the same
+        answer::
+
+            sage: set_random_seed()                              # long time
+            sage: J = random_eja(field=QQ, orthonormalize=False) # long time
+            sage: x = J.random_element()                         # long time
+            sage: while not x.is_invertible():                   # long time
+            ....:     x = J.random_element()                     # long time
+            sage: slow = x.inverse()                             # long time
+            sage: _ = J._charpoly_coefficients()                 # long time
+            sage: fast = x.inverse()                             # long time
+            sage: slow == fast                                   # long time
+            True
         """
         """
-        if not self.is_invertible():
-            raise ValueError("element is not invertible")
+        not_invertible_msg = "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!
+            if self.det().is_zero():
+                raise ZeroDivisionError(not_invertible_msg)
+            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)
+        try:
+            inv = (~self.quadratic_representation())(self)
+            self.is_invertible.set_cache(True)
+            return inv
+        except ZeroDivisionError:
+            self.is_invertible.set_cache(False)
+            raise ZeroDivisionError(not_invertible_msg)
 
 
 
 
+    @cached_method
     def is_invertible(self):
         """
         Return whether or not this element is invertible.
 
         ALGORITHM:
 
     def is_invertible(self):
         """
         Return whether or not this element is invertible.
 
         ALGORITHM:
 
-        The usual way to do this is to check if the determinant is
-        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
-        of this element's minimal polynomial.
-
-        Beware that we can't use the superclass method, because it
-        relies on the algebra being associative.
+        If computing my determinant will be fast, we do so and compare
+        with zero (Proposition II.2.4 in Faraut and
+        Koranyi). Otherwise, Proposition II.3.2 in Faraut and Koranyi
+        reduces the problem to the invertibility of my quadratic
+        representation.
 
         SETUP::
 
 
         SETUP::
 
@@ -541,17 +565,152 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
             sage: J.one().is_invertible()
             True
 
             sage: J.one().is_invertible()
             True
 
-        The zero element is never invertible::
+        The zero element is never invertible in a non-trivial algebra::
 
             sage: set_random_seed()
             sage: J = random_eja()
 
             sage: set_random_seed()
             sage: J = random_eja()
-            sage: J.zero().is_invertible()
+            sage: (not J.is_trivial()) and J.zero().is_invertible()
             False
 
             False
 
+        Test that the fast (cached) and slow algorithms give the same
+        answer::
+
+            sage: set_random_seed()                              # long time
+            sage: J = random_eja(field=QQ, orthonormalize=False) # long time
+            sage: x = J.random_element()                         # long time
+            sage: slow = x.is_invertible()                       # long time
+            sage: _ = J._charpoly_coefficients()                 # long time
+            sage: fast = x.is_invertible()                       # long time
+            sage: slow == fast                                   # long time
+            True
         """
         """
-        zero = self.parent().zero()
-        p = self.minimal_polynomial()
-        return not (p(zero) == zero)
+        if self.is_zero():
+            if self.parent().is_trivial():
+                return True
+            else:
+                return False
+
+        if self.parent()._charpoly_coefficients.is_in_cache():
+            # The determinant will be quicker than inverting the
+            # quadratic representation, most likely.
+            return (not self.det().is_zero())
+
+        # The easiest way to determine if I'm invertible is to try.
+        try:
+            inv = (~self.quadratic_representation())(self)
+            self.inverse.set_cache(inv)
+            return True
+        except ZeroDivisionError:
+            return False
+
+
+    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.monomial(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):
 
 
     def is_nilpotent(self):
@@ -581,10 +740,11 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
 
         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::
@@ -616,7 +776,9 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
             sage: J = JordanSpinEJA(5)
             sage: J.one().is_regular()
             False
             sage: J = JordanSpinEJA(5)
             sage: J.one().is_regular()
             False
-            sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
+            sage: b0, b1, b2, b3, b4 = J.gens()
+            sage: b0 == J.one()
+            True
             sage: for x in J.gens():
             ....:     (J.one() + x).is_regular()
             False
             sage: for x in J.gens():
             ....:     (J.one() + x).is_regular()
             False
@@ -627,19 +789,20 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
 
         TESTS:
 
 
         TESTS:
 
-        The zero element should never be regular::
+        The zero element should never be regular, unless the parent
+        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.zero().is_regular()
-            False
+            sage: J.dimension() <= 1 or not J.zero().is_regular()
+            True
 
         The unit element isn't regular unless the algebra happens to
         consist of only its scalar multiples::
 
             sage: set_random_seed()
             sage: J = random_eja()
 
         The unit element isn't regular unless the algebra happens to
         consist of only its scalar multiples::
 
             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
 
         """
@@ -653,10 +816,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
 
         ALGORITHM:
 
 
         ALGORITHM:
 
-        For now, we skip the messy minimal polynomial computation
-        and instead return the dimension of the vector space spanned
-        by the powers of this element. The latter is a bit more
-        straightforward to compute.
+        .........
 
         SETUP::
 
 
         SETUP::
 
@@ -668,30 +828,33 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
             sage: J = JordanSpinEJA(4)
             sage: J.one().degree()
             1
             sage: J = JordanSpinEJA(4)
             sage: J.one().degree()
             1
-            sage: e0,e1,e2,e3 = J.gens()
-            sage: (e0 - e1).degree()
+            sage: b0,b1,b2,b3 = J.gens()
+            sage: (b0 - b1).degree()
             2
 
         In the spin factor algebra (of rank two), all elements that
         aren't multiples of the identity are regular::
 
             sage: set_random_seed()
             2
 
         In the spin factor algebra (of rank two), all elements that
         aren't multiples of the identity are regular::
 
             sage: set_random_seed()
-            sage: n = ZZ.random_element(1,10)
-            sage: J = JordanSpinEJA(n)
+            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::
 
@@ -701,7 +864,59 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
             True
 
         """
             True
 
         """
-        return self.span_of_powers().dimension()
+        n = self.parent().dimension()
+
+        if n == 0:
+            # The minimal polynomial is an empty product, i.e. the
+            # constant polynomial "1" having degree zero.
+            return 0
+        elif self.is_zero():
+            # The minimal polynomial of zero in a nontrivial algebra
+            # is "t", and is of degree one.
+            return 1
+        elif n == 1:
+            # If this is a nonzero element of a nontrivial algebra, it
+            # has degree at least one. It follows that, in an algebra
+            # of dimension one, the degree must be actually one.
+            return 1
+
+        # BEWARE: The subalgebra_generated_by() method uses the result
+        # of this method to construct a basis for the subalgebra. That
+        # means, in particular, that we cannot implement this method
+        # as ``self.subalgebra_generated_by().dimension()``.
+
+        # Algorithm: keep appending (vector representations of) powers
+        # self as rows to a matrix and echelonizing it. When its rank
+        # stops increasing, we've reached a redundancy.
+
+        # Given the special cases above, we can assume that "self" is
+        # nonzero, the algebra is nontrivial, and that its dimension
+        # is at least two.
+        M = matrix([(self.parent().one()).to_vector()])
+        old_rank = 1
+
+        # Specifying the row-reduction algorithm can e.g. help over
+        # AA because it avoids the RecursionError that gets thrown
+        # when we have to look too hard for a root.
+        #
+        # Beware: QQ supports an entirely different set of "algorithm"
+        # keywords than do AA and RR.
+        algo = None
+        from sage.rings.all import QQ
+        if self.parent().base_ring() is not QQ:
+            algo = "scaled_partial_pivoting"
+
+        for d in range(1,n):
+            M = matrix(M.rows() + [(self**d).to_vector()])
+            M.echelonize(algo)
+            new_rank = M.rank()
+            if new_rank == old_rank:
+                return new_rank
+            else:
+                old_rank = new_rank
+
+        return n
+
 
 
     def left_matrix(self):
 
 
     def left_matrix(self):
@@ -729,18 +944,38 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
         SETUP::
 
             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
@@ -754,16 +989,18 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
         The minimal polynomial and the characteristic polynomial coincide
         and are known (see Alizadeh, Example 11.11) for all elements of
         the spin factor algebra that aren't scalar multiples of the
         The minimal polynomial and the characteristic polynomial coincide
         and are known (see Alizadeh, Example 11.11) for all elements of
         the spin factor algebra that aren't scalar multiples of the
-        identity::
+        identity. We require the dimension of the algebra to be at least
+        two here so that said elements actually exist::
 
             sage: set_random_seed()
 
             sage: set_random_seed()
-            sage: n = ZZ.random_element(2,10)
+            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: while y == y.coefficient(0)*J.one():
             ....:     y = J.random_element()
             sage: J = JordanSpinEJA(n)
             sage: y = J.random_element()
             sage: while y == y.coefficient(0)*J.one():
             ....:     y = J.random_element()
-            sage: y0 = y.vector()[0]
-            sage: y_bar = y.vector()[1:]
+            sage: y0 = y.to_vector()[0]
+            sage: y_bar = y.to_vector()[1:]
             sage: actual = y.minimal_polynomial()
             sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
             sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
             sage: actual = y.minimal_polynomial()
             sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
             sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
@@ -778,68 +1015,150 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
             sage: x.apply_univariate_polynomial(p)
             0
 
             sage: x.apply_univariate_polynomial(p)
             0
 
+        The minimal polynomial is invariant under a change of basis,
+        and in particular, a re-scaling of the basis::
+
+            sage: set_random_seed()
+            sage: n_max = RealSymmetricEJA._max_random_instance_size()
+            sage: n = ZZ.random_element(1, n_max)
+            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: x1.minimal_polynomial() == x2.minimal_polynomial()
+            True
+
         """
         """
-        V = self.span_of_powers()
-        assoc_subalg = self.subalgebra_generated_by()
-        # Mis-design warning: the basis used for span_of_powers()
-        # and subalgebra_generated_by() must be the same, and in
-        # the same order!
-        elt = assoc_subalg(V.coordinates(self.vector()))
-        return elt.operator().minimal_polynomial()
+        if self.is_zero():
+            # Pretty sure we know what the minimal polynomial of
+            # the zero operator is going to be. This ensures
+            # consistency of e.g. the polynomial variable returned
+            # in the "normal" case without us having to think about it.
+            return self.operator().minimal_polynomial()
+
+        # If we don't orthonormalize the subalgebra's basis, then the
+        # first two monomials in the subalgebra will be self^0 and
+        # self^1... assuming that self^1 is not a scalar multiple of
+        # self^0 (the unit element). We special case these to avoid
+        # having to solve a system to coerce self into the subalgebra.
+        A = self.subalgebra_generated_by(orthonormalize=False)
+
+        if A.dimension() == 1:
+            # Does a solve to find the scalar multiple alpha such that
+            # alpha*unit = self. We have to do this because the basis
+            # for the subalgebra will be [ self^0 ], and not [ self^1 ]!
+            unit = self.parent().one()
+            alpha = self.to_vector() / unit.to_vector()
+            return (unit.operator()*alpha).minimal_polynomial()
+        else:
+            # If the dimension of the subalgebra is >= 2, then we just
+            # use the second basis element.
+            return A.monomial(1).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::
 
             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
 
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
-            ....:                                  QuaternionHermitianEJA)
+            ....:                                  HadamardEJA,
+            ....:                                  QuaternionHermitianEJA,
+            ....:                                  RealSymmetricEJA)
 
         EXAMPLES::
 
             sage: J = ComplexHermitianEJA(3)
             sage: J.one()
 
         EXAMPLES::
 
             sage: J = ComplexHermitianEJA(3)
             sage: J.one()
-            e0 + e5 + e8
-            sage: J.one().natural_representation()
-            [1 0 0 0 0 0]
-            [0 1 0 0 0 0]
-            [0 0 1 0 0 0]
-            [0 0 0 1 0 0]
-            [0 0 0 0 1 0]
-            [0 0 0 0 0 1]
+            b0 + b3 + b8
+            sage: J.one().to_matrix()
+            +---+---+---+
+            | 1 | 0 | 0 |
+            +---+---+---+
+            | 0 | 1 | 0 |
+            +---+---+---+
+            | 0 | 0 | 1 |
+            +---+---+---+
 
         ::
 
 
         ::
 
-            sage: J = QuaternionHermitianEJA(3)
+            sage: J = QuaternionHermitianEJA(2)
             sage: J.one()
             sage: J.one()
-            e0 + e9 + 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]
+            b0 + b5
+            sage: J.one().to_matrix()
+            +---+---+
+            | 1 | 0 |
+            +---+---+
+            | 0 | 1 |
+            +---+---+
+
+        This also works in Cartesian product algebras::
+
+            sage: J1 = HadamardEJA(1)
+            sage: J2 = RealSymmetricEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: x = sum(J.gens())
+            sage: x.to_matrix()[0]
+            [1]
+            sage: x.to_matrix()[1]
+            [                  1 0.7071067811865475?]
+            [0.7071067811865475?                   1]
+
+        """
+        B = self.parent().matrix_basis()
+        W = self.parent().matrix_space()
+
+        if hasattr(W, 'cartesian_factors'):
+            # Aaaaand linear combinations don't work in Cartesian
+            # product spaces, even though they provide a method with
+            # that name. This is hidden behind an "if" because the
+            # _scale() function is slow.
+            pairs = zip(B, self.to_vector())
+            return W.sum( _scale(b, alpha) for (b,alpha) in pairs )
+        else:
+            # 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):
         """
         """
-        B = self.parent().natural_basis()
-        W = B[0].matrix_space()
-        return W.linear_combination(zip(self.vector(), B))
+        The norm of this element with respect to :meth:`inner_product`.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            ....:                                  HadamardEJA)
+
+        EXAMPLES::
+
+            sage: J = HadamardEJA(2)
+            sage: x = sum(J.gens())
+            sage: x.norm()
+            1.414213562373095?
+
+        ::
+
+            sage: J = JordanSpinEJA(4)
+            sage: x = sum(J.gens())
+            sage: x.norm()
+            2
+
+        """
+        return self.inner_product(self).sqrt()
 
 
     def operator(self):
 
 
     def operator(self):
@@ -855,8 +1174,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
 
             sage: set_random_seed()
             sage: J = random_eja()
 
             sage: set_random_seed()
             sage: J = random_eja()
-            sage: x = J.random_element()
-            sage: y = J.random_element()
+            sage: x,y = J.random_elements(2)
             sage: x.operator()(y) == x*y
             True
             sage: y.operator()(x) == x*y
             sage: x.operator()(y) == x*y
             True
             sage: y.operator()(x) == x*y
@@ -864,11 +1182,9 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
 
         """
         P = self.parent()
 
         """
         P = self.parent()
-        fda_elt = FiniteDimensionalAlgebraElement(P, self)
-        return FiniteDimensionalEuclideanJordanAlgebraOperator(
-                 P,
-                 P,
-                 fda_elt.matrix().transpose() )
+        left_mult_by_self = lambda y: self*y
+        L = P.module_morphism(function=left_mult_by_self, codomain=P)
+        return FiniteDimensionalEJAOperator(P, P, L.matrix() )
 
 
     def quadratic_representation(self, other=None):
 
 
     def quadratic_representation(self, other=None):
@@ -886,19 +1202,20 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
         Alizadeh's Example 11.12::
 
             sage: set_random_seed()
         Alizadeh's Example 11.12::
 
             sage: set_random_seed()
-            sage: n = ZZ.random_element(1,10)
-            sage: J = JordanSpinEJA(n)
-            sage: x = J.random_element()
-            sage: x_vec = x.vector()
-            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: 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: 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
 
@@ -906,8 +1223,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
 
             sage: set_random_seed()
             sage: J = random_eja()
 
             sage: set_random_seed()
             sage: J = random_eja()
-            sage: x = J.random_element()
-            sage: y = J.random_element()
+            sage: x,y = J.random_elements(2)
             sage: Lx = x.operator()
             sage: Lxx = (x*x).operator()
             sage: Qx = x.quadratic_representation()
             sage: Lx = x.operator()
             sage: Lxx = (x*x).operator()
             sage: Qx = x.quadratic_representation()
@@ -924,7 +1240,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
 
         Property 2 (multiply on the right for :trac:`28272`):
 
 
         Property 2 (multiply on the right for :trac:`28272`):
 
-            sage: alpha = QQ.random_element()
+            sage: alpha = J.base_ring().random_element()
             sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
             True
 
             sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
             True
 
@@ -952,10 +1268,10 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
             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:
@@ -991,35 +1307,119 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
         return ( L*M + M*L - (self*other).operator() )
 
 
         return ( L*M + M*L - (self*other).operator() )
 
 
-    def span_of_powers(self):
-        """
-        Return the vector space spanned by successive powers of
-        this element.
+
+    def spectral_decomposition(self):
         """
         """
-        # The dimension of the subalgebra can't be greater than
-        # the big algebra, so just put everything into a list
-        # and let span() get rid of the excess.
-        #
-        # We do the extra ambient_vector_space() in case we're messing
-        # with polynomials and the direct parent is a module.
-        V = self.parent().vector_space()
-        return V.span( (self**d).vector() for d in xrange(V.dimension()) )
+        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()
+            b0 + b2 + b5
+            sage: J.one().spectral_decomposition()
+            [(1, b0 + b2 + b5)]
+            sage: J.zero().spectral_decomposition()
+            [(0, b0 + b2 + b5)]
+
+        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: (b0, b1, b2, b3, b4, b5, b6, b7, b8, b9) = J.gens()
+            sage: A = 2*b5 - 2*b8
+            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, c2), (1, c0)]
 
 
-    def subalgebra_generated_by(self):
+        """
+        A = self.subalgebra_generated_by(orthonormalize=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, **kwargs):
         """
         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::
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import random_eja
+            sage: from mjo.eja.eja_algebra import (random_eja,
+            ....:                                  HadamardEJA,
+            ....:                                  RealSymmetricEJA)
 
 
-        TESTS::
+        EXAMPLES:
+
+        We can create subalgebras of Cartesian product EJAs that are not
+        themselves Cartesian product EJAs (they're just "regular" EJAs)::
+
+            sage: J1 = HadamardEJA(3)
+            sage: J2 = RealSymmetricEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: J.one().subalgebra_generated_by()
+            Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
+
+        TESTS:
+
+        This subalgebra, being composed of only powers, is associative::
 
             sage: set_random_seed()
 
             sage: set_random_seed()
-            sage: x = random_eja().random_element()
-            sage: x.subalgebra_generated_by().is_associative()
+            sage: x0 = random_eja().random_element()
+            sage: A = x0.subalgebra_generated_by()
+            sage: x,y,z = A.random_elements(3)
+            sage: (x*y)*z == x*(y*z)
             True
 
         Squaring in the subalgebra should work the same as in
             True
 
         Squaring in the subalgebra should work the same as in
@@ -1027,54 +1427,29 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
 
             sage: set_random_seed()
             sage: x = random_eja().random_element()
 
             sage: set_random_seed()
             sage: x = random_eja().random_element()
-            sage: u = x.subalgebra_generated_by().random_element()
-            sage: u.operator()(u) == u^2
+            sage: A = x.subalgebra_generated_by()
+            sage: A(x^2) == A(x)*A(x)
+            True
+
+        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: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
             True
 
         """
             True
 
         """
-        # First get the subspace spanned by the powers of myself...
-        V = self.span_of_powers()
-        F = self.base_ring()
-
-        # Now figure out the entries of the right-multiplication
-        # matrix for the successive basis elements b0, b1,... of
-        # that subspace.
-        mats = []
-        for b_right in V.basis():
-            eja_b_right = self.parent()(b_right)
-            b_right_rows = []
-            # The first row of the right-multiplication matrix by
-            # b1 is what we get if we apply that matrix to b1. The
-            # second row of the right multiplication matrix by b1
-            # is what we get when we apply that matrix to b2...
-            #
-            # IMPORTANT: this assumes that all vectors are COLUMN
-            # vectors, unlike our superclass (which uses row vectors).
-            for b_left in V.basis():
-                eja_b_left = self.parent()(b_left)
-                # Multiply in the original EJA, but then get the
-                # coordinates from the subalgebra in terms of its
-                # basis.
-                this_row = V.coordinates((eja_b_left*eja_b_right).vector())
-                b_right_rows.append(this_row)
-            b_right_matrix = matrix(F, b_right_rows)
-            mats.append(b_right_matrix)
-
-        # It's an algebra of polynomials in one element, and EJAs
-        # are power-associative.
-        #
-        # TODO: choose generator names intelligently.
-        #
-        # The rank is the highest possible degree of a minimal polynomial,
-        # and is bounded above by the dimension. We know in this case that
-        # there's an element whose minimal polynomial has the same degree
-        # as the space's dimension, so that must be its rank too.
-        return FiniteDimensionalEuclideanJordanAlgebra(
-                 F,
-                 mats,
-                 V.dimension(),
-                 assume_associative=True,
-                 names='f')
+        powers = tuple( self**k for k in range(self.degree()) )
+        A = self.parent().subalgebra(powers,
+                                     associative=True,
+                                     check_field=False,
+                                     check_axioms=False,
+                                     **kwargs)
+        A.one.set_cache(A(self.parent().one()))
+        return A
 
 
     def subalgebra_idempotent(self):
 
 
     def subalgebra_idempotent(self):
@@ -1086,33 +1461,36 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
 
             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!")
 
-        V = self.span_of_powers()
         J = self.subalgebra_generated_by()
         J = self.subalgebra_generated_by()
-        # Mis-design warning: the basis used for span_of_powers()
-        # and subalgebra_generated_by() must be the same, and in
-        # the same order!
-        u = J(V.coordinates(self.vector()))
+        u = J(self)
 
         # The image of the matrix of left-u^m-multiplication
         # will be minimal for some natural number s...
         s = 0
 
         # The image of the matrix of left-u^m-multiplication
         # will be minimal for some natural number s...
         s = 0
-        minimal_dim = V.dimension()
-        for i in xrange(1, V.dimension()):
+        minimal_dim = J.dimension()
+        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
@@ -1131,29 +1509,33 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
         # Our FiniteDimensionalAlgebraElement superclass uses rows.
         u_next = u**(s+1)
         A = u_next.operator().matrix()
         # Our FiniteDimensionalAlgebraElement superclass uses rows.
         u_next = u**(s+1)
         A = u_next.operator().matrix()
-        c_coordinates = A.solve_right(u_next.vector())
+        c = J.from_vector(A.solve_right(u_next.to_vector()))
 
 
-        # Now c_coordinates is the idempotent we want, but it's in
-        # the coordinate system of the subalgebra.
-        #
-        # We need the basis for J, but as elements of the parent algebra.
-        #
-        basis = [self.parent(v) for v in V.basis()]
-        return self.parent().linear_combination(zip(c_coordinates, basis))
+        # Now c is the idempotent we want, but it still lives in the subalgebra.
+        return c.superalgebra_element()
 
 
     def trace(self):
         """
         Return my trace, the sum of my eigenvalues.
 
 
 
     def trace(self):
         """
         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()
@@ -1161,7 +1543,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
 
         ::
 
 
         ::
 
-            sage: J = RealCartesianProductEJA(5)
+            sage: J = HadamardEJA(5)
             sage: J.one().trace()
             5
 
             sage: J.one().trace()
             5
 
@@ -1171,18 +1553,33 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
 
             sage: set_random_seed()
             sage: J = random_eja()
 
             sage: set_random_seed()
             sage: J = random_eja()
-            sage: J.random_element().trace() in J.base_ring()
+            sage: J.random_element().trace() in RLF
+            True
+
+        The trace is linear::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: x,y = J.random_elements(2)
+            sage: alpha = J.base_ring().random_element()
+            sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace()
             True
 
         """
         P = self.parent()
         r = P.rank()
             True
 
         """
         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,
         # we want the negative of THAT for the trace.
         # 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,
         # we want the negative of THAT for the trace.
-        return -p(*self.vector())
+        return -p(*self.to_vector())
 
 
     def trace_inner_product(self, other):
 
 
     def trace_inner_product(self, other):
@@ -1195,22 +1592,16 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
 
         TESTS:
 
 
         TESTS:
 
-        The trace inner product is commutative::
+        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()
-            sage: x = J.random_element(); y = J.random_element()
+            sage: x,y,z = J.random_elements(3)
+            sage: # commutative
             sage: x.trace_inner_product(y) == y.trace_inner_product(x)
             True
             sage: x.trace_inner_product(y) == y.trace_inner_product(x)
             True
-
-        The trace inner product is bilinear::
-
-            sage: set_random_seed()
-            sage: J = random_eja()
-            sage: x = J.random_element()
-            sage: y = J.random_element()
-            sage: z = J.random_element()
-            sage: a = QQ.random_element();
+            sage: # bilinear
+            sage: a = J.base_ring().random_element();
             sage: actual = (a*(x+z)).trace_inner_product(y)
             sage: expected = ( a*x.trace_inner_product(y) +
             ....:              a*z.trace_inner_product(y) )
             sage: actual = (a*(x+z)).trace_inner_product(y)
             sage: expected = ( a*x.trace_inner_product(y) +
             ....:              a*z.trace_inner_product(y) )
@@ -1221,15 +1612,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
             ....:              a*x.trace_inner_product(z) )
             sage: actual == expected
             True
             ....:              a*x.trace_inner_product(z) )
             sage: actual == expected
             True
-
-        The trace inner product satisfies the compatibility
-        condition in the definition of a Euclidean Jordan algebra::
-
-            sage: set_random_seed()
-            sage: J = random_eja()
-            sage: x = J.random_element()
-            sage: y = J.random_element()
-            sage: z = J.random_element()
+            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
 
@@ -1238,3 +1621,30 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(FiniteDimensionalAlgebraEle
             raise TypeError("'other' must live in the same algebra")
 
         return (self*other).trace()
             raise TypeError("'other' must live in the same algebra")
 
         return (self*other).trace()
+
+
+    def trace_norm(self):
+        """
+        The norm of this element with respect to :meth:`trace_inner_product`.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            ....:                                  HadamardEJA)
+
+        EXAMPLES::
+
+            sage: J = HadamardEJA(2)
+            sage: x = sum(J.gens())
+            sage: x.trace_norm()
+            1.414213562373095?
+
+        ::
+
+            sage: J = JordanSpinEJA(4)
+            sage: x = sum(J.gens())
+            sage: x.trace_norm()
+            2.828427124746190?
+
+        """
+        return self.trace_inner_product(self).sqrt()