]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: begin rework of Cartesian product EJA.
[sage.d.git] / mjo / eja / eja_algebra.py
index ea62da8ca796031bd99382c5886bb1091d145b78..c7c0df8de0e67aef03b391099ea6a72374dce13a 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
@@ -133,7 +135,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             deortho_vector_basis = tuple( V(b.list()) for b in basis )
 
             from mjo.eja.eja_utils import gram_schmidt
-            basis = gram_schmidt(basis, inner_product)
+            basis = tuple(gram_schmidt(basis, inner_product))
 
         # Save the (possibly orthonormalized) matrix basis for
         # later...
@@ -158,7 +160,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         # Now we actually compute the multiplication and inner-product
         # tables/matrices using the possibly-orthonormalized basis.
-        self._inner_product_matrix = matrix.zero(field, n)
+        self._inner_product_matrix = matrix.identity(field, n)
         self._multiplication_table = [ [0 for j in range(i+1)]
                                        for i in range(n) ]
 
@@ -171,15 +173,20 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                 q_i = basis[i]
                 q_j = basis[j]
 
-                elt = jordan_product(q_i, q_j)
-                ip = inner_product(q_i, q_j)
-
                 # 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()))
                 self._multiplication_table[i][j] = self.from_vector(elt)
-                self._inner_product_matrix[i,j] = ip
-                self._inner_product_matrix[j,i] = ip
+
+                if not orthonormalize:
+                    # If we're orthonormalizing the basis with respect
+                    # to an inner-product, then the inner-product
+                    # matrix with respect to the resulting basis is
+                    # just going to be the identity.
+                    ip = inner_product(q_i, q_j)
+                    self._inner_product_matrix[i,j] = ip
+                    self._inner_product_matrix[j,i] = ip
 
         self._inner_product_matrix._cache = {'hermitian': True}
         self._inner_product_matrix.set_immutable()
@@ -315,7 +322,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         This method should of course always return ``True``, unless
         this algebra was constructed with ``check_axioms=False`` and
-        passed an invalid multiplication table.
+        passed an invalid Jordan or inner-product.
         """
 
         # Used to check whether or not something is zero in an inexact
@@ -682,7 +689,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.
 
@@ -732,7 +739,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
@@ -745,23 +752,57 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: from mjo.eja.eja_algebra import (HadamardEJA,
             ....:                                  random_eja)
 
-        EXAMPLES::
+        EXAMPLES:
+
+        We can compute unit element in the Hadamard EJA::
+
+            sage: J = HadamardEJA(5)
+            sage: J.one()
+            e0 + e1 + e2 + e3 + e4
+
+        The unit element in the Hadamard EJA is inherited in the
+        subalgebras generated by its elements::
 
             sage: J = HadamardEJA(5)
             sage: J.one()
             e0 + e1 + e2 + e3 + e4
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
+            sage: A.one()
+            f0
+            sage: A.one().superalgebra_element()
+            e0 + e1 + e2 + e3 + e4
 
         TESTS:
 
-        The identity element acts like the identity::
+        The identity element acts like the identity, regardless of
+        whether or not we orthonormalize::
 
             sage: set_random_seed()
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: J.one()*x == x and x*J.one() == x
             True
+            sage: A = x.subalgebra_generated_by()
+            sage: y = A.random_element()
+            sage: A.one()*y == y and y*A.one() == y
+            True
+
+        ::
+
+            sage: set_random_seed()
+            sage: J = random_eja(field=QQ, orthonormalize=False)
+            sage: x = J.random_element()
+            sage: J.one()*x == x and x*J.one() == x
+            True
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
+            sage: y = A.random_element()
+            sage: A.one()*y == y and y*A.one() == y
+            True
 
-        The matrix of the unit element's operator is the identity::
+        The matrix of the unit element's operator is the identity,
+        regardless of the base field and whether or not we
+        orthonormalize::
 
             sage: set_random_seed()
             sage: J = random_eja()
@@ -769,6 +810,27 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: expected = matrix.identity(J.base_ring(), J.dimension())
             sage: actual == expected
             True
+            sage: x = J.random_element()
+            sage: A = x.subalgebra_generated_by()
+            sage: actual = A.one().operator().matrix()
+            sage: expected = matrix.identity(A.base_ring(), A.dimension())
+            sage: actual == expected
+            True
+
+        ::
+
+            sage: set_random_seed()
+            sage: J = random_eja(field=QQ, orthonormalize=False)
+            sage: actual = J.one().operator().matrix()
+            sage: expected = matrix.identity(J.base_ring(), J.dimension())
+            sage: actual == expected
+            True
+            sage: x = J.random_element()
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
+            sage: actual = A.one().operator().matrix()
+            sage: expected = matrix.identity(A.base_ring(), A.dimension())
+            sage: actual == expected
+            True
 
         Ensure that the cached unit element (often precomputed by
         hand) agrees with the computed one::
@@ -780,6 +842,15 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: J.one() == cached
             True
 
+        ::
+
+            sage: set_random_seed()
+            sage: J = random_eja(field=QQ, orthonormalize=False)
+            sage: cached = J.one()
+            sage: J.one.clear_cache()
+            sage: J.one() == cached
+            True
+
         """
         # We can brute-force compute the matrices of the operators
         # that correspond to the basis elements of this algebra.
@@ -1032,6 +1103,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()
@@ -1067,10 +1153,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):
@@ -1131,7 +1224,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
@@ -1187,9 +1280,7 @@ class RationalBasisEJA(FiniteDimensionalEJA):
                  jordan_product,
                  inner_product,
                  field=AA,
-                 orthonormalize=True,
                  check_field=True,
-                 check_axioms=True,
                  **kwargs):
 
         if check_field:
@@ -1198,6 +1289,7 @@ class RationalBasisEJA(FiniteDimensionalEJA):
             if not all( all(b_i in QQ for b_i in b.list()) for b in basis ):
                 raise TypeError("basis not rational")
 
+        self._rational_algebra = None
         if field is not QQ:
             # There's no point in constructing the extra algebra if this
             # one is already rational.
@@ -1212,15 +1304,13 @@ class RationalBasisEJA(FiniteDimensionalEJA):
                                        field=QQ,
                                        orthonormalize=False,
                                        check_field=False,
-                                       check_axioms=False,
-                                       **kwargs)
+                                       check_axioms=False)
 
         super().__init__(basis,
                          jordan_product,
                          inner_product,
                          field=field,
                          check_field=check_field,
-                         check_axioms=check_axioms,
                          **kwargs)
 
     @cached_method
@@ -1260,7 +1350,14 @@ class RationalBasisEJA(FiniteDimensionalEJA):
         a = ( a_i.change_ring(self.base_ring())
               for a_i in self._rational_algebra._charpoly_coefficients() )
 
-        # Now convert the coordinate variables back to the
+        if self._deortho_matrix is None:
+            # This can happen if our base ring was, say, AA and we
+            # chose not to (or didn't need to) orthonormalize. It's
+            # still faster to do the computations over QQ even if
+            # the numbers in the boxes stay the same.
+            return tuple(a)
+
+        # Otherwise, convert the coordinate variables back to the
         # deorthonormalized ones.
         R = self.coordinate_polynomial_ring()
         from sage.modules.free_module_element import vector
@@ -1594,6 +1691,38 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
 
 
 class ComplexMatrixEJA(MatrixEJA):
+    # A manual dictionary-cache for the complex_extension() method,
+    # since apparently @classmethods can't also be @cached_methods.
+    _complex_extension = {}
+
+    @classmethod
+    def complex_extension(cls,field):
+        r"""
+        The complex field that we embed/unembed, as an extension
+        of the given ``field``.
+        """
+        if field in cls._complex_extension:
+            return cls._complex_extension[field]
+
+        # Sage doesn't know how to adjoin the complex "i" (the root of
+        # x^2 + 1) to a field in a general way. Here, we just enumerate
+        # all of the cases that I have cared to support so far.
+        if field is AA:
+            # Sage doesn't know how to embed AA into QQbar, i.e. how
+            # to adjoin sqrt(-1) to AA.
+            F = QQbar
+        elif not field.is_exact():
+            # RDF or RR
+            F = field.complex_field()
+        else:
+            # Works for QQ and... maybe some other fields.
+            R = PolynomialRing(field, 'z')
+            z = R.gen()
+            F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
+
+        cls._complex_extension[field] = F
+        return F
+
     @staticmethod
     def dimension_over_reals():
         return 2
@@ -1648,9 +1777,10 @@ class ComplexMatrixEJA(MatrixEJA):
 
         blocks = []
         for z in M.list():
-            a = z.list()[0] # real part, I guess
-            b = z.list()[1] # imag part, I guess
-            blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
+            a = z.real()
+            b = z.imag()
+            blocks.append(matrix(field, 2, [ [ a, b],
+                                             [-b, a] ]))
 
         return matrix.block(field, n, blocks)
 
@@ -1689,26 +1819,7 @@ class ComplexMatrixEJA(MatrixEJA):
         super(ComplexMatrixEJA,cls).real_unembed(M)
         n = ZZ(M.nrows())
         d = cls.dimension_over_reals()
-
-        # If "M" was normalized, its base ring might have roots
-        # adjoined and they can stick around after unembedding.
-        field = M.base_ring()
-        R = PolynomialRing(field, 'z')
-        z = R.gen()
-
-        # Sage doesn't know how to adjoin the complex "i" (the root of
-        # x^2 + 1) to a field in a general way. Here, we just enumerate
-        # all of the cases that I have cared to support so far.
-        if field is AA:
-            # Sage doesn't know how to embed AA into QQbar, i.e. how
-            # to adjoin sqrt(-1) to AA.
-            F = QQbar
-        elif not field.is_exact():
-            # RDF or RR
-            F = field.complex_field()
-        else:
-            # Works for QQ and... maybe some other fields.
-            F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
+        F = cls.complex_extension(M.base_ring())
         i = F.gen()
 
         # Go top-left to bottom-right (reading order), converting every
@@ -1804,7 +1915,6 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
-            sage: field = QuadraticField(2, 'sqrt2')
             sage: B = ComplexHermitianEJA._denormalized_basis(n)
             sage: all( M.is_symmetric() for M in  B)
             True
@@ -1822,18 +1932,27 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
         #   * The diagonal will (as a result) be real.
         #
         S = []
+        Eij = matrix.zero(F,n)
         for i in range(n):
             for j in range(i+1):
-                Eij = matrix(F, n, lambda k,l: k==i and l==j)
+                # "build" E_ij
+                Eij[i,j] = 1
                 if i == j:
                     Sij = cls.real_embed(Eij)
                     S.append(Sij)
                 else:
                     # The second one has a minus because it's conjugated.
-                    Sij_real = cls.real_embed(Eij + Eij.transpose())
+                    Eij[j,i] = 1 # Eij = Eij + Eij.transpose()
+                    Sij_real = cls.real_embed(Eij)
                     S.append(Sij_real)
-                    Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose())
+                    # Eij = I*Eij - I*Eij.transpose()
+                    Eij[i,j] = I
+                    Eij[j,i] = -I
+                    Sij_imag = cls.real_embed(Eij)
                     S.append(Sij_imag)
+                    Eij[j,i] = 0
+                # "erase" E_ij
+                Eij[i,j] = 0
 
         # Since we embedded these, we can drop back to the "field" that we
         # started with instead of the complex extension "F".
@@ -1869,6 +1988,25 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
         return cls(n, **kwargs)
 
 class QuaternionMatrixEJA(MatrixEJA):
+
+    # A manual dictionary-cache for the quaternion_extension() method,
+    # since apparently @classmethods can't also be @cached_methods.
+    _quaternion_extension = {}
+
+    @classmethod
+    def quaternion_extension(cls,field):
+        r"""
+        The quaternion field that we embed/unembed, as an extension
+        of the given ``field``.
+        """
+        if field in cls._quaternion_extension:
+            return cls._quaternion_extension[field]
+
+        Q = QuaternionAlgebra(field,-1,-1)
+
+        cls._quaternion_extension[field] = Q
+        return Q
+
     @staticmethod
     def dimension_over_reals():
         return 4
@@ -1973,8 +2111,7 @@ class QuaternionMatrixEJA(MatrixEJA):
 
         # Use the base ring of the matrix to ensure that its entries can be
         # multiplied by elements of the quaternion algebra.
-        field = M.base_ring()
-        Q = QuaternionAlgebra(field,-1,-1)
+        Q = cls.quaternion_extension(M.base_ring())
         i,j,k = Q.gens()
 
         # Go top-left to bottom-right (reading order), converting every
@@ -2089,23 +2226,39 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
         #   * The diagonal will (as a result) be real.
         #
         S = []
+        Eij = matrix.zero(Q,n)
         for i in range(n):
             for j in range(i+1):
-                Eij = matrix(Q, n, lambda k,l: k==i and l==j)
+                # "build" E_ij
+                Eij[i,j] = 1
                 if i == j:
                     Sij = cls.real_embed(Eij)
                     S.append(Sij)
                 else:
                     # The second, third, and fourth ones have a minus
                     # because they're conjugated.
-                    Sij_real = cls.real_embed(Eij + Eij.transpose())
+                    # Eij = Eij + Eij.transpose()
+                    Eij[j,i] = 1
+                    Sij_real = cls.real_embed(Eij)
                     S.append(Sij_real)
-                    Sij_I = cls.real_embed(I*Eij - I*Eij.transpose())
+                    # Eij = I*(Eij - Eij.transpose())
+                    Eij[i,j] = I
+                    Eij[j,i] = -I
+                    Sij_I = cls.real_embed(Eij)
                     S.append(Sij_I)
-                    Sij_J = cls.real_embed(J*Eij - J*Eij.transpose())
+                    # Eij = J*(Eij - Eij.transpose())
+                    Eij[i,j] = J
+                    Eij[j,i] = -J
+                    Sij_J = cls.real_embed(Eij)
                     S.append(Sij_J)
-                    Sij_K = cls.real_embed(K*Eij - K*Eij.transpose())
+                    # Eij = K*(Eij - Eij.transpose())
+                    Eij[i,j] = K
+                    Eij[j,i] = -K
+                    Sij_K = cls.real_embed(Eij)
                     S.append(Sij_K)
+                    Eij[j,i] = 0
+                # "erase" E_ij
+                Eij[i,j] = 0
 
         # Since we embedded these, we can drop back to the "field" that we
         # started with instead of the quaternion algebra "Q".
@@ -2311,31 +2464,38 @@ class BilinearFormEJA(ConcreteEJA):
         ....:              for j in range(n-1) ]
         sage: actual == expected
         True
+
     """
     def __init__(self, B, **kwargs):
-        if not B.is_positive_definite():
-            raise ValueError("bilinear form is not positive-definite")
+        # The matrix "B" is supplied by the user in most cases,
+        # so it makes sense to check whether or not its positive-
+        # definite unless we are specifically asked not to...
+        if ("check_axioms" not in kwargs) or kwargs["check_axioms"]:
+            if not B.is_positive_definite():
+                raise ValueError("bilinear form is not positive-definite")
+
+        # However, all of the other data for this EJA is computed
+        # by us in manner that guarantees the axioms are
+        # satisfied. So, again, unless we are specifically asked to
+        # verify things, we'll skip the rest of the checks.
+        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
         def inner_product(x,y):
-            return (B*x).inner_product(y)
+            return (y.T*B*x)[0,0]
 
         def jordan_product(x,y):
             P = x.parent()
-            x0 = x[0]
-            xbar = x[1:]
-            y0 = y[0]
-            ybar = y[1:]
-            z0 = inner_product(x,y)
+            x0 = x[0,0]
+            xbar = x[1:,0]
+            y0 = y[0,0]
+            ybar = y[1:,0]
+            z0 = inner_product(y,x)
             zbar = y0*xbar + x0*ybar
-            return P((z0,) + tuple(zbar))
-
-        # We know this is a valid EJA, but will double-check
-        # if the user passes check_axioms=True.
-        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
+            return P([z0] + zbar.list())
 
         n = B.nrows()
-        standard_basis = FreeModule(ZZ, n).basis()
-        super(BilinearFormEJA, self).__init__(standard_basis,
+        column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )
+        super(BilinearFormEJA, self).__init__(column_basis,
                                               jordan_product,
                                               inner_product,
                                               **kwargs)
@@ -2516,243 +2676,361 @@ class TrivialEJA(ConcreteEJA):
         # inappropriate for us.
         return cls(**kwargs)
 
-class DirectSumEJA(ConcreteEJA):
+
+class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
+                          FiniteDimensionalEJA):
     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.
+    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,
+        sage: from mjo.eja.eja_algebra import (CartesianProductEJA,
         ....:                                  HadamardEJA,
-        ....:                                  RealSymmetricEJA,
-        ....:                                  DirectSumEJA)
+        ....:                                  JordanSpinEJA,
+        ....:                                  RealSymmetricEJA)
 
-    EXAMPLES::
+    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(3)
-        sage: J = DirectSumEJA(J1,J2)
-        sage: J.dimension()
-        8
-        sage: J.rank()
-        5
+        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)
 
     TESTS:
 
-    The external direct sum construction is only valid when the two factors
-    have the same base ring; an error is raised otherwise::
+    All factors must share the same base field::
 
-        sage: set_random_seed()
-        sage: J1 = random_eja(field=AA)
-        sage: J2 = random_eja(field=QQ,orthonormalize=False)
-        sage: J = DirectSumEJA(J1,J2)
+        sage: J1 = HadamardEJA(2, field=QQ)
+        sage: J2 = RealSymmetricEJA(2)
+        sage: CartesianProductEJA((J1,J2))
         Traceback (most recent call last):
         ...
-        ValueError: algebras must share the same base field
-
+        ValueError: all factors 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)
+    def __init__(self, modules, **kwargs):
+        CombinatorialFreeModule_CartesianProduct.__init__(self, modules, **kwargs)
+        field = modules[0].base_ring()
+        if not all( J.base_ring() == field for J in modules ):
+            raise ValueError("all factors must share the same base field")
+
+        M = cartesian_product( [J.matrix_space() for J in modules] )
+
+        m = len(modules)
+        W = VectorSpace(field,m)
+        self._matrix_basis = []
+        for k in range(m):
+            for a in modules[k].matrix_basis():
+                v = W.zero().list()
+                v[k] = a
+                self._matrix_basis.append(M(v))
+
+        self._matrix_basis = tuple(self._matrix_basis)
+
+        n = len(self._matrix_basis)
+        # TODO:
+        #
+        # Initialize the FDEJA class, too. Does this override the
+        # initialization that we did for the
+        # CombinatorialFreeModule_CartesianProduct class? If not, we
+        # will probably have to duplicate some of the work (i.e. one
+        # of the constructors).  Since the CartesianProduct one is
+        # smaller, that makes the most sense to copy/paste if it comes
+        # down to that.
+        #
 
-        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.
+        self.rank.set_cache(sum(J.rank() for J in modules))
 
+    @cached_method
+    def cartesian_projection(self, i):
+        r"""
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import (HadamardEJA,
-            ....:                                  JordanSpinEJA,
-            ....:                                  DirectSumEJA)
+            sage: from mjo.eja.eja_algebra import (random_eja,
+            ....:                                  HadamardEJA,
+            ....:                                  RealSymmetricEJA)
 
         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
+            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
 
-    def projections(self):
-        r"""
-        Return a pair of projections onto this algebra's factors.
+        TESTS:
 
-        SETUP::
+        The answer never changes::
 
-            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
-            ....:                                  ComplexHermitianEJA,
-            ....:                                  DirectSumEJA)
+            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
 
-        EXAMPLES::
+        """
+        Ji = self.cartesian_factors()[i]
+        # We reimplement the CombinatorialFreeModule superclass method
+        # because if we don't, something gets messed up with the caching
+        # and the answer changes the second time you run it. See the TESTS.
+        Pi = self._module_morphism(lambda j_t: Ji.monomial(j_t[1])
+                                   if i == j_t[0] else Ji.zero(),
+                                   codomain=Ji)
+        return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix())
 
-            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):
+    @cached_method
+    def cartesian_embedding(self, i):
         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)
+            ....:                                  HadamardEJA,
+            ....:                                  RealSymmetricEJA)
 
         EXAMPLES::
 
-            sage: J1 = JordanSpinEJA(3)
+            sage: J1 = HadamardEJA(2)
             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)
+            sage: J = cartesian_product([J1,J2])
+            sage: J
+            foo
+            sage: J.cartesian_embedding
+            bar
+            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
 
         TESTS:
 
-        Composing a projection with the corresponding inclusion should
-        produce the identity map, and mismatching them should produce
-        the zero map::
+        The answer never changes::
 
             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()
+            sage: J = cartesian_product([J1,J2])
+            sage: E0 = J.cartesian_embedding(0)
+            sage: E1 = J.cartesian_embedding(0)
+            sage: E0 == E1
             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))
+        Ji = self.cartesian_factors()[i]
+        # We reimplement the CombinatorialFreeModule superclass method
+        # because if we don't, something gets messed up with the caching
+        # and the answer changes the second time you run it. See the TESTS.
+        Ei = Ji._module_morphism(lambda t: self.monomial((i, t)), codomain=self)
+        return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix())
+
+
+FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA
+
+
+#     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))