]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
Revert "eja: factor out a separate class used by cartesian_product()."
[sage.d.git] / mjo / eja / eja_algebra.py
index 5910221ef8f8807dce7fee30eeb2da9045114399..55b8207ddc3f881e8d98f245bf0275bc4b947334 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
@@ -62,7 +64,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                  associative=False,
                  check_field=True,
                  check_axioms=True,
-                 prefix='e'):
+                 prefix='e',
+                 category=None):
 
         if check_field:
             if not field.is_subring(RR):
@@ -91,11 +94,12 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                 raise ValueError("inner-product is not commutative")
 
 
-        category = MagmaticAlgebras(field).FiniteDimensional()
-        category = category.WithBasis().Unital()
-        if associative:
-            # Element subalgebras can take advantage of this.
-            category = category.Associative()
+        if category is None:
+            category = MagmaticAlgebras(field).FiniteDimensional()
+            category = category.WithBasis().Unital()
+            if associative:
+                # Element subalgebras can take advantage of this.
+                category = category.Associative()
 
         # Call the superclass constructor so that we can use its from_vector()
         # method to build our multiplication table.
@@ -687,7 +691,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 +741,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
@@ -1101,6 +1105,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()
@@ -1136,10 +1155,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 +1226,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
@@ -1667,6 +1693,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
@@ -1721,9 +1779,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)
 
@@ -1762,26 +1821,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
@@ -1950,6 +1990,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
@@ -2054,8 +2113,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
@@ -2620,244 +2678,386 @@ 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 (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
+
+    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
+
+    """
+    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")
+
+        basis = tuple( b.to_vector().column() for b in self.basis() )
+
+        # Define jordan/inner products that operate on the basis.
+        def jordan_product(x_mat,y_mat):
+            x = self.from_vector(_mat2vec(x_mat))
+            y = self.from_vector(_mat2vec(y_mat))
+            return self.cartesian_jordan_product(x,y).to_vector().column()
+
+        def inner_product(x_mat, y_mat):
+            x = self.from_vector(_mat2vec(x_mat))
+            y = self.from_vector(_mat2vec(y_mat))
+            return self.cartesian_inner_product(x,y)
+
+        # Use whatever category the superclass came up with. Usually
+        # some join of the EJA and Cartesian product
+        # categories. There's no need to check the field since it
+        # already came from an EJA. Likewise the axioms are guaranteed
+        # to be satisfied. We can't orthonormalize by default because
+        # there's no way to pass "orthonormalize=False" to
+        # cartesian_product(...) when the base ring is QQ and
+        # orthonormalizing would give us irrational entries.
+        #
+        # TODO: create a separate constructor that is capable of
+        # orthonormalizing and is only used by the cartesian_product()
+        # thingy.
+        FiniteDimensionalEJA.__init__(self,
+                                      basis,
+                                      jordan_product,
+                                      inner_product,
+                                      field=field,
+                                      orthonormalize=False,
+                                      check_field=False,
+                                      check_axioms=False,
+                                      category=self.category())
+        # 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.
+        #
+
+        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 (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())
+
+
+    def cartesian_jordan_product(self, x, y):
+        r"""
+        The componentwise Jordan product.
+
+        We project ``x`` and ``y`` onto our factors, and add up the
+        Jordan products from the subalgebras. This may still be useful
+        after (if) the default Jordan product in the Cartesian product
+        algebra is overridden.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            ....:                                  JordanSpinEJA)
+
+        EXAMPLE::
+
+            sage: J1 = HadamardEJA(3)
+            sage: J2 = JordanSpinEJA(3)
+            sage: J = cartesian_product([J1,J2])
+            sage: x1 = J1.from_vector(vector(QQ,(1,2,1)))
+            sage: y1 = J1.from_vector(vector(QQ,(1,0,2)))
+            sage: x2 = J2.from_vector(vector(QQ,(1,2,3)))
+            sage: y2 = J2.from_vector(vector(QQ,(1,1,1)))
+            sage: z1 = J.from_vector(vector(QQ,(1,2,1,1,2,3)))
+            sage: z2 = J.from_vector(vector(QQ,(1,0,2,1,1,1)))
+            sage: (x1*y1).to_vector()
+            (1, 0, 2)
+            sage: (x2*y2).to_vector()
+            (6, 3, 4)
+            sage: J.cartesian_jordan_product(z1,z2).to_vector()
+            (1, 0, 2, 6, 3, 4)
+
+        """
+        m = len(self.cartesian_factors())
+        projections = ( self.cartesian_projection(i) for i in range(m) )
+        products = ( P(x)*P(y) for P in projections )
+        return self._cartesian_product_of_elements(tuple(products))
+
+    def cartesian_inner_product(self, x, y):
+        r"""
+        The standard componentwise Cartesian inner-product.
+
+        We project ``x`` and ``y`` onto our factors, and add up the
+        inner-products from the subalgebras. This may still be useful
+        after (if) the default inner product in the Cartesian product
+        algebra is overridden.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            ....:                                  QuaternionHermitianEJA)
+
+        EXAMPLE::
+
+            sage: J1 = HadamardEJA(3,field=QQ)
+            sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: J = cartesian_product([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: z1 = J._cartesian_product_of_elements((x1,y1))
+            sage: z2 = J._cartesian_product_of_elements((x2,y2))
+            sage: J.cartesian_inner_product(z1,z2)
+            5
+
+        """
+        m = len(self.cartesian_factors())
+        projections = ( self.cartesian_projection(i) for i in range(m) )
+        return sum( P(x).inner_product(P(y)) for P in projections )
+
+
+    Element = FiniteDimensionalEJAElement
 
 
+FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA
 random_eja = ConcreteEJA.random_instance