]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: remove a completed TODO, and doctest it.
[sage.d.git] / mjo / eja / eja_algebra.py
index 58d5ce03a495ddf2eb7a9c4cba0042145be0c037..3a832c1051694dc8e19d6c8dfb002ceda1203279 100644 (file)
@@ -20,7 +20,9 @@ from itertools import repeat
 
 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
 from sage.categories.magmatic_algebras import MagmaticAlgebras
-from sage.combinat.free_module import CombinatorialFreeModule
+from sage.categories.sets_cat import cartesian_product
+from sage.combinat.free_module import (CombinatorialFreeModule,
+                                       CombinatorialFreeModule_CartesianProduct)
 from sage.matrix.constructor import matrix
 from sage.matrix.matrix_space import MatrixSpace
 from sage.misc.cachefunc import cached_method
@@ -31,7 +33,7 @@ from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
                             QuadraticField)
 from mjo.eja.eja_element import FiniteDimensionalEJAElement
 from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
-from mjo.eja.eja_utils import _mat2vec
+from mjo.eja.eja_utils import _all2list, _mat2vec
 
 class FiniteDimensionalEJA(CombinatorialFreeModule):
     r"""
@@ -39,7 +41,15 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
     INPUT:
 
-      - basis -- a tuple of basis elements in their matrix form.
+      - basis -- a tuple of basis elements in "matrix form," which
+        must be the same form as the arguments to ``jordan_product``
+        and ``inner_product``. In reality, "matrix form" can be either
+        vectors, matrices, or a Cartesian product (ordered tuple)
+        of vectors or matrices. All of these would ideally be vector
+        spaces in sage with no special-casing needed; but in reality
+        we turn vectors into column-matrices and Cartesian products
+        `(a,b)` into column matrices `(a,b)^{T}` after converting
+        `a` and `b` themselves.
 
       - jordan_product -- function of two elements (in matrix form)
         that returns their jordan product in this algebra; this will
@@ -60,10 +70,23 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                  field=AA,
                  orthonormalize=True,
                  associative=False,
+                 cartesian_product=False,
                  check_field=True,
                  check_axioms=True,
                  prefix='e'):
 
+        # Keep track of whether or not the matrix basis consists of
+        # tuples, since we need special cases for them damned near
+        # everywhere.  This is INDEPENDENT of whether or not the
+        # algebra is a cartesian product, since a subalgebra of a
+        # cartesian product will have a basis of tuples, but will not
+        # in general itself be a cartesian product algebra.
+        self._matrix_basis_is_cartesian = False
+        n = len(basis)
+        if n > 0:
+            if hasattr(basis[0], 'cartesian_factors'):
+                self._matrix_basis_is_cartesian = True
+
         if check_field:
             if not field.is_subring(RR):
                 # Note: this does return true for the real algebraic
@@ -73,7 +96,18 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         # If the basis given to us wasn't over the field that it's
         # supposed to be over, fix that. Or, you know, crash.
-        basis = tuple( b.change_ring(field) for b in basis )
+        if not cartesian_product:
+            # The field for a cartesian product algebra comes from one
+            # of its factors and is the same for all factors, so
+            # there's no need to "reapply" it on product algebras.
+            if self._matrix_basis_is_cartesian:
+                # OK since if n == 0, the basis does not consist of tuples.
+                P = basis[0].parent()
+                basis = tuple( P(tuple(b_i.change_ring(field) for b_i in b))
+                               for b in basis )
+            else:
+                basis = tuple( b.change_ring(field) for b in basis )
+
 
         if check_axioms:
             # Check commutativity of the Jordan and inner-products.
@@ -96,15 +130,17 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         if associative:
             # Element subalgebras can take advantage of this.
             category = category.Associative()
+        if cartesian_product:
+            category = category.CartesianProducts()
 
         # Call the superclass constructor so that we can use its from_vector()
         # method to build our multiplication table.
-        n = len(basis)
-        super().__init__(field,
-                         range(n),
-                         prefix=prefix,
-                         category=category,
-                         bracket=False)
+        CombinatorialFreeModule.__init__(self,
+                                         field,
+                                         range(n),
+                                         prefix=prefix,
+                                         category=category,
+                                         bracket=False)
 
         # Now comes all of the hard work. We'll be constructing an
         # ambient vector space V that our (vectorized) basis lives in,
@@ -115,8 +151,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         degree = 0
         if n > 0:
-            # Works on both column and square matrices...
-            degree = len(basis[0].list())
+            degree = len(_all2list(basis[0]))
 
         # Build an ambient space that fits our matrix basis when
         # written out as "long vectors."
@@ -130,7 +165,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             # Save a copy of the un-orthonormalized basis for later.
             # Convert it to ambient V (vector) coordinates while we're
             # at it, because we'd have to do it later anyway.
-            deortho_vector_basis = tuple( V(b.list()) for b in basis )
+            deortho_vector_basis = tuple( V(_all2list(b)) for b in basis )
 
             from mjo.eja.eja_utils import gram_schmidt
             basis = tuple(gram_schmidt(basis, inner_product))
@@ -142,7 +177,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # Now create the vector space for the algebra, which will have
         # its own set of non-ambient coordinates (in terms of the
         # supplied basis).
-        vector_basis = tuple( V(b.list()) for b in basis )
+        vector_basis = tuple( V(_all2list(b)) for b in basis )
         W = V.span_of_basis( vector_basis, check=check_axioms)
 
         if orthonormalize:
@@ -174,7 +209,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                 # The jordan product returns a matrixy answer, so we
                 # have to convert it to the algebra coordinates.
                 elt = jordan_product(q_i, q_j)
-                elt = W.coordinate_vector(V(elt.list()))
+                elt = W.coordinate_vector(V(_all2list(elt)))
                 self._multiplication_table[i][j] = self.from_vector(elt)
 
                 if not orthonormalize:
@@ -279,22 +314,32 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: y = J.random_element()
             sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
             True
+
         """
         B = self._inner_product_matrix
         return (B*x.to_vector()).inner_product(y.to_vector())
 
 
-    def _is_commutative(self):
+    def is_associative(self):
         r"""
-        Whether or not this algebra's multiplication table is commutative.
+        Return whether or not this algebra's Jordan product is associative.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
+
+        EXAMPLES::
+
+            sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
+            sage: J.is_associative()
+            False
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
+            sage: A.is_associative()
+            True
 
-        This method should of course always return ``True``, unless
-        this algebra was constructed with ``check_axioms=False`` and
-        passed an invalid multiplication table.
         """
-        return all( self.product_on_basis(i,j) == self.product_on_basis(i,j)
-                    for i in range(self.dimension())
-                    for j in range(self.dimension()) )
+        return "Associative" in self.category().axioms()
 
     def _is_jordanian(self):
         r"""
@@ -303,13 +348,13 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         We only check one arrangement of `x` and `y`, so for a
         ``True`` result to be truly true, you should also check
-        :meth:`_is_commutative`. This method should of course always
+        :meth:`is_commutative`. This method should of course always
         return ``True``, unless this algebra was constructed with
         ``check_axioms=False`` and passed an invalid multiplication table.
         """
-        return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
+        return all( (self.gens()[i]**2)*(self.gens()[i]*self.gens()[j])
                     ==
-                    (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j))
+                    (self.gens()[i])*((self.gens()[i]**2)*self.gens()[j])
                     for i in range(self.dimension())
                     for j in range(self.dimension()) )
 
@@ -331,9 +376,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         for i in range(self.dimension()):
             for j in range(self.dimension()):
                 for k in range(self.dimension()):
-                    x = self.monomial(i)
-                    y = self.monomial(j)
-                    z = self.monomial(k)
+                    x = self.gens()[i]
+                    y = self.gens()[j]
+                    z = self.gens()[k]
                     diff = (x*y).inner_product(z) - x.inner_product(y*z)
 
                     if self.base_ring().is_exact():
@@ -377,6 +422,15 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             ...
             ValueError: not an element of this algebra
 
+        Tuples work as well, provided that the matrix basis for the
+        algebra consists of them::
+
+            sage: J1 = HadamardEJA(3)
+            sage: J2 = RealSymmetricEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
+            e(0, 1) + e(1, 2)
+
         TESTS:
 
         Ensure that we can convert any element of the two non-matrix
@@ -393,13 +447,23 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: J(x.to_vector().column()) == x
             True
 
+        We cannot coerce elements between algebras just because their
+        matrix representations are compatible::
+
+            sage: J1 = HadamardEJA(3)
+            sage: J2 = JordanSpinEJA(3)
+            sage: J2(J1.one())
+            Traceback (most recent call last):
+            ...
+            ValueError: not an element of this algebra
+            sage: J1(J2.zero())
+            Traceback (most recent call last):
+            ...
+            ValueError: not an element of this algebra
+
         """
         msg = "not an element of this algebra"
-        if elt == 0:
-            # The superclass implementation of random_element()
-            # needs to be able to coerce "0" into the algebra.
-            return self.zero()
-        elif elt in self.base_ring():
+        if elt in self.base_ring():
             # Ensure that no base ring -> algebra coercion is performed
             # by this method. There's some stupidity in sage that would
             # otherwise propagate to this method; for example, sage thinks
@@ -407,9 +471,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             raise ValueError(msg)
 
         try:
+            # Try to convert a vector into a column-matrix...
             elt = elt.column()
         except (AttributeError, TypeError):
-            # Try to convert a vector into a column-matrix
+            # and ignore failure, because we weren't really expecting
+            # a vector as an argument anyway.
             pass
 
         if elt not in self.matrix_space():
@@ -422,14 +488,20 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # closure whereas the base ring of the 3-by-3 identity matrix
         # could be QQ instead of QQbar.
         #
+        # And, we also have to handle Cartesian product bases (when
+        # the matrix basis consists of tuples) here. The "good news"
+        # is that we're already converting everything to long vectors,
+        # and that strategy works for tuples as well.
+        #
         # We pass check=False because the matrix basis is "guaranteed"
         # to be linearly independent... right? Ha ha.
-        V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols())
-        W = V.span_of_basis( (_mat2vec(s) for s in self.matrix_basis()),
+        elt = _all2list(elt)
+        V = VectorSpace(self.base_ring(), len(elt))
+        W = V.span_of_basis( (V(_all2list(s)) for s in self.matrix_basis()),
                              check=False)
 
         try:
-            coords =  W.coordinate_vector(_mat2vec(elt))
+            coords = W.coordinate_vector(V(elt))
         except ArithmeticError:  # vector is not in free module
             raise ValueError(msg)
 
@@ -656,8 +728,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         # And to each subsequent row, prepend an entry that belongs to
         # the left-side "header column."
-        M += [ [self.monomial(i)] + [ self.product_on_basis(i,j)
-                                      for j in range(n) ]
+        M += [ [self.gens()[i]] + [ self.product_on_basis(i,j)
+                                    for j in range(n) ]
                for i in range(n) ]
 
         return table(M, header_row=True, header_column=True, frame=True)
@@ -687,7 +759,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         Why implement this for non-matrix algebras? Avoiding special
         cases for the :class:`BilinearFormEJA` pays with simplicity in
         its own right. But mainly, we would like to be able to assume
-        that elements of a :class:`DirectSumEJA` can be displayed
+        that elements of a :class:`CartesianProductEJA` can be displayed
         nicely, without having to have special classes for direct sums
         one of whose components was a matrix algebra.
 
@@ -737,7 +809,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         if self.is_trivial():
             return MatrixSpace(self.base_ring(), 0)
         else:
-            return self._matrix_basis[0].matrix_space()
+            return self.matrix_basis()[0].parent()
 
 
     @cached_method
@@ -993,14 +1065,12 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         if not c.is_idempotent():
             raise ValueError("element is not idempotent: %s" % c)
 
-        from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
-
         # Default these to what they should be if they turn out to be
         # trivial, because eigenspaces_left() won't return eigenvalues
         # corresponding to trivial spaces (e.g. it returns only the
         # eigenspace corresponding to lambda=1 if you take the
         # decomposition relative to the identity element).
-        trivial = FiniteDimensionalEJASubalgebra(self, ())
+        trivial = self.subalgebra(())
         J0 = trivial                          # eigenvalue zero
         J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
         J1 = trivial                          # eigenvalue one
@@ -1010,9 +1080,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                 J5 = eigspace
             else:
                 gens = tuple( self.from_vector(b) for b in eigspace.basis() )
-                subalg = FiniteDimensionalEJASubalgebra(self,
-                                                        gens,
-                                                        check_axioms=False)
+                subalg = self.subalgebra(gens, check_axioms=False)
                 if eigval == 0:
                     J0 = subalg
                 elif eigval == 1:
@@ -1101,6 +1169,21 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         r"""
         The `r` polynomial coefficients of the "characteristic polynomial
         of" function.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import random_eja
+
+        TESTS:
+
+        The theory shows that these are all homogeneous polynomials of
+        a known degree::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
+            True
+
         """
         n = self.dimension()
         R = self.coordinate_polynomial_ring()
@@ -1110,7 +1193,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         def L_x_i_j(i,j):
             # From a result in my book, these are the entries of the
             # basis representation of L_x.
-            return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
+            return sum( vars[k]*self.gens()[k].operator().matrix()[i,j]
                         for k in range(n) )
 
         L_x = matrix(F, n, n, L_x_i_j)
@@ -1136,10 +1219,17 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         # The theory says that only the first "r" coefficients are
         # nonzero, and they actually live in the original polynomial
-        # ring and not the fraction field. We negate them because
-        # in the actual characteristic polynomial, they get moved
-        # to the other side where x^r lives.
-        return -A_rref.solve_right(E*b).change_ring(R)[:r]
+        # ring and not the fraction field. We negate them because in
+        # the actual characteristic polynomial, they get moved to the
+        # other side where x^r lives. We don't bother to trim A_rref
+        # down to a square matrix and solve the resulting system,
+        # because the upper-left r-by-r portion of A_rref is
+        # guaranteed to be the identity matrix, so e.g.
+        #
+        #   A_rref.solve_right(Y)
+        #
+        # would just be returning Y.
+        return (-E*b)[:r].change_ring(R)
 
     @cached_method
     def rank(self):
@@ -1200,7 +1290,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: set_random_seed()    # long time
             sage: J = random_eja()     # long time
-            sage: caches = J.rank()    # long time
+            sage: cached = J.rank()    # long time
             sage: J.rank.clear_cache() # long time
             sage: J.rank() == cached   # long time
             True
@@ -1209,6 +1299,14 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         return len(self._charpoly_coefficients())
 
 
+    def subalgebra(self, basis, **kwargs):
+        r"""
+        Create a subalgebra of this algebra from the given basis.
+        """
+        from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
+        return FiniteDimensionalEJASubalgebra(self, basis, **kwargs)
+
+
     def vector_space(self):
         """
         Return the vector space that underlies this algebra.
@@ -1227,7 +1325,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         return self.zero().to_vector().parent().ambient_vector_space()
 
 
-    Element = FiniteDimensionalEJAElement
 
 class RationalBasisEJA(FiniteDimensionalEJA):
     r"""
@@ -2336,7 +2433,11 @@ class HadamardEJA(ConcreteEJA):
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
         column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )
-        super().__init__(column_basis, jordan_product, inner_product, **kwargs)
+        super().__init__(column_basis,
+                         jordan_product,
+                         inner_product,
+                         associative=True,
+                         **kwargs)
         self.rank.set_cache(n)
 
         if n == 0:
@@ -2652,244 +2753,395 @@ class TrivialEJA(ConcreteEJA):
         # inappropriate for us.
         return cls(**kwargs)
 
-# class DirectSumEJA(ConcreteEJA):
-#     r"""
-#     The external (orthogonal) direct sum of two other Euclidean Jordan
-#     algebras. Essentially the Cartesian product of its two factors.
-#     Every Euclidean Jordan algebra decomposes into an orthogonal
-#     direct sum of simple Euclidean Jordan algebras, so no generality
-#     is lost by providing only this construction.
-
-#     SETUP::
-
-#         sage: from mjo.eja.eja_algebra import (random_eja,
-#         ....:                                  HadamardEJA,
-#         ....:                                  RealSymmetricEJA,
-#         ....:                                  DirectSumEJA)
-
-#     EXAMPLES::
-
-#         sage: J1 = HadamardEJA(2)
-#         sage: J2 = RealSymmetricEJA(3)
-#         sage: J = DirectSumEJA(J1,J2)
-#         sage: J.dimension()
-#         8
-#         sage: J.rank()
-#         5
-
-#     TESTS:
-
-#     The external direct sum construction is only valid when the two factors
-#     have the same base ring; an error is raised otherwise::
-
-#         sage: set_random_seed()
-#         sage: J1 = random_eja(field=AA)
-#         sage: J2 = random_eja(field=QQ,orthonormalize=False)
-#         sage: J = DirectSumEJA(J1,J2)
-#         Traceback (most recent call last):
-#         ...
-#         ValueError: algebras must share the same base field
-
-#     """
-#     def __init__(self, J1, J2, **kwargs):
-#         if J1.base_ring() != J2.base_ring():
-#             raise ValueError("algebras must share the same base field")
-#         field = J1.base_ring()
-
-#         self._factors = (J1, J2)
-#         n1 = J1.dimension()
-#         n2 = J2.dimension()
-#         n = n1+n2
-#         V = VectorSpace(field, n)
-#         mult_table = [ [ V.zero() for j in range(i+1) ]
-#                        for i in range(n) ]
-#         for i in range(n1):
-#             for j in range(i+1):
-#                 p = (J1.monomial(i)*J1.monomial(j)).to_vector()
-#                 mult_table[i][j] = V(p.list() + [field.zero()]*n2)
-
-#         for i in range(n2):
-#             for j in range(i+1):
-#                 p = (J2.monomial(i)*J2.monomial(j)).to_vector()
-#                 mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
-
-#         # TODO: build the IP table here from the two constituent IP
-#         # matrices (it'll be block diagonal, I think).
-#         ip_table = [ [ field.zero() for j in range(i+1) ]
-#                        for i in range(n) ]
-#         super(DirectSumEJA, self).__init__(field,
-#                                            mult_table,
-#                                            ip_table,
-#                                            check_axioms=False,
-#                                            **kwargs)
-#         self.rank.set_cache(J1.rank() + J2.rank())
-
-
-#     def factors(self):
-#         r"""
-#         Return the pair of this algebra's factors.
-
-#         SETUP::
-
-#             sage: from mjo.eja.eja_algebra import (HadamardEJA,
-#             ....:                                  JordanSpinEJA,
-#             ....:                                  DirectSumEJA)
-
-#         EXAMPLES::
-
-#             sage: J1 = HadamardEJA(2, field=QQ)
-#             sage: J2 = JordanSpinEJA(3, field=QQ)
-#             sage: J = DirectSumEJA(J1,J2)
-#             sage: J.factors()
-#             (Euclidean Jordan algebra of dimension 2 over Rational Field,
-#              Euclidean Jordan algebra of dimension 3 over Rational Field)
-
-#         """
-#         return self._factors
-
-#     def projections(self):
-#         r"""
-#         Return a pair of projections onto this algebra's factors.
-
-#         SETUP::
-
-#             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
-#             ....:                                  ComplexHermitianEJA,
-#             ....:                                  DirectSumEJA)
-
-#         EXAMPLES::
-
-#             sage: J1 = JordanSpinEJA(2)
-#             sage: J2 = ComplexHermitianEJA(2)
-#             sage: J = DirectSumEJA(J1,J2)
-#             sage: (pi_left, pi_right) = J.projections()
-#             sage: J.one().to_vector()
-#             (1, 0, 1, 0, 0, 1)
-#             sage: pi_left(J.one()).to_vector()
-#             (1, 0)
-#             sage: pi_right(J.one()).to_vector()
-#             (1, 0, 0, 1)
-
-#         """
-#         (J1,J2) = self.factors()
-#         m = J1.dimension()
-#         n = J2.dimension()
-#         V_basis = self.vector_space().basis()
-#         # Need to specify the dimensions explicitly so that we don't
-#         # wind up with a zero-by-zero matrix when we want e.g. a
-#         # zero-by-two matrix (important for composing things).
-#         P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
-#         P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
-#         pi_left = FiniteDimensionalEJAOperator(self,J1,P1)
-#         pi_right = FiniteDimensionalEJAOperator(self,J2,P2)
-#         return (pi_left, pi_right)
-
-#     def inclusions(self):
-#         r"""
-#         Return the pair of inclusion maps from our factors into us.
-
-#         SETUP::
-
-#             sage: from mjo.eja.eja_algebra import (random_eja,
-#             ....:                                  JordanSpinEJA,
-#             ....:                                  RealSymmetricEJA,
-#             ....:                                  DirectSumEJA)
-
-#         EXAMPLES::
-
-#             sage: J1 = JordanSpinEJA(3)
-#             sage: J2 = RealSymmetricEJA(2)
-#             sage: J = DirectSumEJA(J1,J2)
-#             sage: (iota_left, iota_right) = J.inclusions()
-#             sage: iota_left(J1.zero()) == J.zero()
-#             True
-#             sage: iota_right(J2.zero()) == J.zero()
-#             True
-#             sage: J1.one().to_vector()
-#             (1, 0, 0)
-#             sage: iota_left(J1.one()).to_vector()
-#             (1, 0, 0, 0, 0, 0)
-#             sage: J2.one().to_vector()
-#             (1, 0, 1)
-#             sage: iota_right(J2.one()).to_vector()
-#             (0, 0, 0, 1, 0, 1)
-#             sage: J.one().to_vector()
-#             (1, 0, 0, 1, 0, 1)
-
-#         TESTS:
-
-#         Composing a projection with the corresponding inclusion should
-#         produce the identity map, and mismatching them should produce
-#         the zero map::
-
-#             sage: set_random_seed()
-#             sage: J1 = random_eja()
-#             sage: J2 = random_eja()
-#             sage: J = DirectSumEJA(J1,J2)
-#             sage: (iota_left, iota_right) = J.inclusions()
-#             sage: (pi_left, pi_right) = J.projections()
-#             sage: pi_left*iota_left == J1.one().operator()
-#             True
-#             sage: pi_right*iota_right == J2.one().operator()
-#             True
-#             sage: (pi_left*iota_right).is_zero()
-#             True
-#             sage: (pi_right*iota_left).is_zero()
-#             True
-
-#         """
-#         (J1,J2) = self.factors()
-#         m = J1.dimension()
-#         n = J2.dimension()
-#         V_basis = self.vector_space().basis()
-#         # Need to specify the dimensions explicitly so that we don't
-#         # wind up with a zero-by-zero matrix when we want e.g. a
-#         # two-by-zero matrix (important for composing things).
-#         I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m])
-#         I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:])
-#         iota_left = FiniteDimensionalEJAOperator(J1,self,I1)
-#         iota_right = FiniteDimensionalEJAOperator(J2,self,I2)
-#         return (iota_left, iota_right)
-
-#     def inner_product(self, x, y):
-#         r"""
-#         The standard Cartesian inner-product.
-
-#         We project ``x`` and ``y`` onto our factors, and add up the
-#         inner-products from the subalgebras.
-
-#         SETUP::
-
-
-#             sage: from mjo.eja.eja_algebra import (HadamardEJA,
-#             ....:                                  QuaternionHermitianEJA,
-#             ....:                                  DirectSumEJA)
-
-#         EXAMPLE::
-
-#             sage: J1 = HadamardEJA(3,field=QQ)
-#             sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
-#             sage: J = DirectSumEJA(J1,J2)
-#             sage: x1 = J1.one()
-#             sage: x2 = x1
-#             sage: y1 = J2.one()
-#             sage: y2 = y1
-#             sage: x1.inner_product(x2)
-#             3
-#             sage: y1.inner_product(y2)
-#             2
-#             sage: J.one().inner_product(J.one())
-#             5
-
-#         """
-#         (pi_left, pi_right) = self.projections()
-#         x1 = pi_left(x)
-#         x2 = pi_right(x)
-#         y1 = pi_left(y)
-#         y2 = pi_right(y)
-
-#         return (x1.inner_product(y1) + x2.inner_product(y2))
+
+class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
+                          FiniteDimensionalEJA):
+    r"""
+    The external (orthogonal) direct sum of two or more Euclidean
+    Jordan algebras. Every Euclidean Jordan algebra decomposes into an
+    orthogonal direct sum of simple Euclidean Jordan algebras which is
+    then isometric to a Cartesian product, so no generality is lost by
+    providing only this construction.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import (random_eja,
+        ....:                                  CartesianProductEJA,
+        ....:                                  HadamardEJA,
+        ....:                                  JordanSpinEJA,
+        ....:                                  RealSymmetricEJA)
+
+    EXAMPLES:
+
+    The Jordan product is inherited from our factors and implemented by
+    our CombinatorialFreeModule Cartesian product superclass::
+
+        sage: set_random_seed()
+        sage: J1 = HadamardEJA(2)
+        sage: J2 = RealSymmetricEJA(2)
+        sage: J = cartesian_product([J1,J2])
+        sage: x,y = J.random_elements(2)
+        sage: x*y in J
+        True
+
+    The ability to retrieve the original factors is implemented by our
+    CombinatorialFreeModule Cartesian product superclass::
+
+        sage: J1 = HadamardEJA(2, field=QQ)
+        sage: J2 = JordanSpinEJA(3, field=QQ)
+        sage: J = cartesian_product([J1,J2])
+        sage: J.cartesian_factors()
+        (Euclidean Jordan algebra of dimension 2 over Rational Field,
+         Euclidean Jordan algebra of dimension 3 over Rational Field)
+
+    You can provide more than two factors::
+
+        sage: J1 = HadamardEJA(2)
+        sage: J2 = JordanSpinEJA(3)
+        sage: J3 = RealSymmetricEJA(3)
+        sage: cartesian_product([J1,J2,J3])
+        Euclidean Jordan algebra of dimension 2 over Algebraic Real
+        Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
+        Real Field (+) Euclidean Jordan algebra of dimension 6 over
+        Algebraic Real Field
+
+    Rank is additive on a Cartesian product::
+
+        sage: J1 = HadamardEJA(1)
+        sage: J2 = RealSymmetricEJA(2)
+        sage: J = cartesian_product([J1,J2])
+        sage: J1.rank.clear_cache()
+        sage: J2.rank.clear_cache()
+        sage: J.rank.clear_cache()
+        sage: J.rank()
+        3
+        sage: J.rank() == J1.rank() + J2.rank()
+        True
+
+    The same rank computation works over the rationals, with whatever
+    basis you like::
+
+        sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
+        sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
+        sage: J = cartesian_product([J1,J2])
+        sage: J1.rank.clear_cache()
+        sage: J2.rank.clear_cache()
+        sage: J.rank.clear_cache()
+        sage: J.rank()
+        3
+        sage: J.rank() == J1.rank() + J2.rank()
+        True
+
+    The product algebra will be associative if and only if all of its
+    components are associative::
+
+        sage: J1 = HadamardEJA(2)
+        sage: J1.is_associative()
+        True
+        sage: J2 = HadamardEJA(3)
+        sage: J2.is_associative()
+        True
+        sage: J3 = RealSymmetricEJA(3)
+        sage: J3.is_associative()
+        False
+        sage: CP1 = cartesian_product([J1,J2])
+        sage: CP1.is_associative()
+        True
+        sage: CP2 = cartesian_product([J1,J3])
+        sage: CP2.is_associative()
+        False
+
+    TESTS:
+
+    All factors must share the same base field::
+
+        sage: J1 = HadamardEJA(2, field=QQ)
+        sage: J2 = RealSymmetricEJA(2)
+        sage: CartesianProductEJA((J1,J2))
+        Traceback (most recent call last):
+        ...
+        ValueError: all factors must share the same base field
+
+    The cached unit element is the same one that would be computed::
+
+        sage: set_random_seed()              # long time
+        sage: J1 = random_eja()              # long time
+        sage: J2 = random_eja()              # long time
+        sage: J = cartesian_product([J1,J2]) # long time
+        sage: actual = J.one()               # long time
+        sage: J.one.clear_cache()            # long time
+        sage: expected = J.one()             # long time
+        sage: actual == expected             # long time
+        True
+
+    """
+    Element = FiniteDimensionalEJAElement
+
+
+    def __init__(self, algebras, **kwargs):
+        CombinatorialFreeModule_CartesianProduct.__init__(self,
+                                                          algebras,
+                                                          **kwargs)
+        field = algebras[0].base_ring()
+        if not all( J.base_ring() == field for J in algebras ):
+            raise ValueError("all factors must share the same base field")
+
+        associative = all( m.is_associative() for m in algebras )
+
+        # The definition of matrix_space() and self.basis() relies
+        # only on the stuff in the CFM_CartesianProduct class, which
+        # we've already initialized.
+        Js = self.cartesian_factors()
+        m = len(Js)
+        MS = self.matrix_space()
+        basis = tuple(
+            MS(tuple( self.cartesian_projection(i)(b).to_matrix()
+                      for i in range(m) ))
+            for b in self.basis()
+        )
+
+        # Define jordan/inner products that operate on that matrix_basis.
+        def jordan_product(x,y):
+            return MS(tuple(
+                (Js[i](x[i])*Js[i](y[i])).to_matrix() for i in range(m)
+            ))
+
+        def inner_product(x, y):
+            return sum(
+                Js[i](x[i]).inner_product(Js[i](y[i])) for i in range(m)
+            )
+
+        # There's no need to check the field since it already came
+        # from an EJA. Likewise the axioms are guaranteed to be
+        # satisfied, unless the guy writing this class sucks.
+        #
+        # If you want the basis to be orthonormalized, orthonormalize
+        # the factors.
+        FiniteDimensionalEJA.__init__(self,
+                                      basis,
+                                      jordan_product,
+                                      inner_product,
+                                      field=field,
+                                      orthonormalize=False,
+                                      associative=associative,
+                                      cartesian_product=True,
+                                      check_field=False,
+                                      check_axioms=False)
+
+        ones = tuple(J.one() for J in algebras)
+        self.one.set_cache(self._cartesian_product_of_elements(ones))
+        self.rank.set_cache(sum(J.rank() for J in algebras))
+
+    def matrix_space(self):
+        r"""
+        Return the space that our matrix basis lives in as a Cartesian
+        product.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            ....:                                  RealSymmetricEJA)
+
+        EXAMPLES::
+
+            sage: J1 = HadamardEJA(1)
+            sage: J2 = RealSymmetricEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: J.matrix_space()
+            The Cartesian product of (Full MatrixSpace of 1 by 1 dense
+            matrices over Algebraic Real Field, Full MatrixSpace of 2
+            by 2 dense matrices over Algebraic Real Field)
+
+        """
+        from sage.categories.cartesian_product import cartesian_product
+        return cartesian_product( [J.matrix_space()
+                                   for J in self.cartesian_factors()] )
+
+    @cached_method
+    def cartesian_projection(self, i):
+        r"""
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (random_eja,
+            ....:                                  JordanSpinEJA,
+            ....:                                  HadamardEJA,
+            ....:                                  RealSymmetricEJA,
+            ....:                                  ComplexHermitianEJA)
+
+        EXAMPLES:
+
+        The projection morphisms are Euclidean Jordan algebra
+        operators::
+
+            sage: J1 = HadamardEJA(2)
+            sage: J2 = RealSymmetricEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: J.cartesian_projection(0)
+            Linear operator between finite-dimensional Euclidean Jordan
+            algebras represented by the matrix:
+            [1 0 0 0 0]
+            [0 1 0 0 0]
+            Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
+            Real Field (+) Euclidean Jordan algebra of dimension 3 over
+            Algebraic Real Field
+            Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
+            Real Field
+            sage: J.cartesian_projection(1)
+            Linear operator between finite-dimensional Euclidean Jordan
+            algebras represented by the matrix:
+            [0 0 1 0 0]
+            [0 0 0 1 0]
+            [0 0 0 0 1]
+            Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
+            Real Field (+) Euclidean Jordan algebra of dimension 3 over
+            Algebraic Real Field
+            Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
+            Real Field
+
+        The projections work the way you'd expect on the vector
+        representation of an element::
+
+            sage: J1 = JordanSpinEJA(2)
+            sage: J2 = ComplexHermitianEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: pi_left = J.cartesian_projection(0)
+            sage: pi_right = J.cartesian_projection(1)
+            sage: pi_left(J.one()).to_vector()
+            (1, 0)
+            sage: pi_right(J.one()).to_vector()
+            (1, 0, 0, 1)
+            sage: J.one().to_vector()
+            (1, 0, 1, 0, 0, 1)
+
+        TESTS:
+
+        The answer never changes::
+
+            sage: set_random_seed()
+            sage: J1 = random_eja()
+            sage: J2 = random_eja()
+            sage: J = cartesian_product([J1,J2])
+            sage: P0 = J.cartesian_projection(0)
+            sage: P1 = J.cartesian_projection(0)
+            sage: P0 == P1
+            True
+
+        """
+        Ji = self.cartesian_factors()[i]
+        # Requires the fix on Trac 31421/31422 to work!
+        Pi = super().cartesian_projection(i)
+        return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix())
+
+    @cached_method
+    def cartesian_embedding(self, i):
+        r"""
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (random_eja,
+            ....:                                  JordanSpinEJA,
+            ....:                                  HadamardEJA,
+            ....:                                  RealSymmetricEJA)
+
+        EXAMPLES:
+
+        The embedding morphisms are Euclidean Jordan algebra
+        operators::
+
+            sage: J1 = HadamardEJA(2)
+            sage: J2 = RealSymmetricEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: J.cartesian_embedding(0)
+            Linear operator between finite-dimensional Euclidean Jordan
+            algebras represented by the matrix:
+            [1 0]
+            [0 1]
+            [0 0]
+            [0 0]
+            [0 0]
+            Domain: Euclidean Jordan algebra of dimension 2 over
+            Algebraic Real Field
+            Codomain: Euclidean Jordan algebra of dimension 2 over
+            Algebraic Real Field (+) Euclidean Jordan algebra of
+            dimension 3 over Algebraic Real Field
+            sage: J.cartesian_embedding(1)
+            Linear operator between finite-dimensional Euclidean Jordan
+            algebras represented by the matrix:
+            [0 0 0]
+            [0 0 0]
+            [1 0 0]
+            [0 1 0]
+            [0 0 1]
+            Domain: Euclidean Jordan algebra of dimension 3 over
+            Algebraic Real Field
+            Codomain: Euclidean Jordan algebra of dimension 2 over
+            Algebraic Real Field (+) Euclidean Jordan algebra of
+            dimension 3 over Algebraic Real Field
+
+        The embeddings work the way you'd expect on the vector
+        representation of an element::
+
+            sage: J1 = JordanSpinEJA(3)
+            sage: J2 = RealSymmetricEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: iota_left = J.cartesian_embedding(0)
+            sage: iota_right = J.cartesian_embedding(1)
+            sage: iota_left(J1.zero()) == J.zero()
+            True
+            sage: iota_right(J2.zero()) == J.zero()
+            True
+            sage: J1.one().to_vector()
+            (1, 0, 0)
+            sage: iota_left(J1.one()).to_vector()
+            (1, 0, 0, 0, 0, 0)
+            sage: J2.one().to_vector()
+            (1, 0, 1)
+            sage: iota_right(J2.one()).to_vector()
+            (0, 0, 0, 1, 0, 1)
+            sage: J.one().to_vector()
+            (1, 0, 0, 1, 0, 1)
+
+        TESTS:
+
+        The answer never changes::
+
+            sage: set_random_seed()
+            sage: J1 = random_eja()
+            sage: J2 = random_eja()
+            sage: J = cartesian_product([J1,J2])
+            sage: E0 = J.cartesian_embedding(0)
+            sage: E1 = J.cartesian_embedding(0)
+            sage: E0 == E1
+            True
+
+        Composing a projection with the corresponding inclusion should
+        produce the identity map, and mismatching them should produce
+        the zero map::
+
+            sage: set_random_seed()
+            sage: J1 = random_eja()
+            sage: J2 = random_eja()
+            sage: J = cartesian_product([J1,J2])
+            sage: iota_left = J.cartesian_embedding(0)
+            sage: iota_right = J.cartesian_embedding(1)
+            sage: pi_left = J.cartesian_projection(0)
+            sage: pi_right = J.cartesian_projection(1)
+            sage: pi_left*iota_left == J1.one().operator()
+            True
+            sage: pi_right*iota_right == J2.one().operator()
+            True
+            sage: (pi_left*iota_right).is_zero()
+            True
+            sage: (pi_right*iota_left).is_zero()
+            True
+
+        """
+        Ji = self.cartesian_factors()[i]
+        # Requires the fix on Trac 31421/31422 to work!
+        Ei = super().cartesian_embedding(i)
+        return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix())
+
 
 
+FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA
 
 random_eja = ConcreteEJA.random_instance