]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: add subalgebra() method.
[sage.d.git] / mjo / eja / eja_algebra.py
index 39703dd1b65d416865af044f41164e6ac4c3ee48..8005412e65da8b6192072420dffe794c7975b5d0 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
@@ -29,7 +31,8 @@ from sage.modules.free_module import FreeModule, VectorSpace
 from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
                             PolynomialRing,
                             QuadraticField)
-from mjo.eja.eja_element import FiniteDimensionalEJAElement
+from mjo.eja.eja_element import (CartesianProductEJAElement,
+                                 FiniteDimensionalEJAElement)
 from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
 from mjo.eja.eja_utils import _mat2vec
 
@@ -39,7 +42,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,6 +71,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                  field=AA,
                  orthonormalize=True,
                  associative=False,
+                 cartesian_product=False,
                  check_field=True,
                  check_axioms=True,
                  prefix='e'):
@@ -73,7 +85,12 @@ 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.
+            basis = tuple( b.change_ring(field) for b in basis )
+
 
         if check_axioms:
             # Check commutativity of the Jordan and inner-products.
@@ -96,6 +113,8 @@ 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.
@@ -113,10 +132,17 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # we see in things like x = 1*e1 + 2*e2.
         vector_basis = basis
 
+        def flatten(b):
+            # flatten a vector, matrix, or cartesian product of those
+            # things into a long list.
+            if cartesian_product:
+                return sum(( b_i.list() for b_i in b ), [])
+            else:
+                return b.list()
+
         degree = 0
         if n > 0:
-            # Works on both column and square matrices...
-            degree = len(basis[0].list())
+            degree = len(flatten(basis[0]))
 
         # Build an ambient space that fits our matrix basis when
         # written out as "long vectors."
@@ -130,7 +156,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(flatten(b)) for b in basis )
 
             from mjo.eja.eja_utils import gram_schmidt
             basis = tuple(gram_schmidt(basis, inner_product))
@@ -142,7 +168,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(flatten(b)) for b in basis )
         W = V.span_of_basis( vector_basis, check=check_axioms)
 
         if orthonormalize:
@@ -174,7 +200,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(flatten(elt)))
                 self._multiplication_table[i][j] = self.from_vector(elt)
 
                 if not orthonormalize:
@@ -279,22 +305,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 +339,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 +367,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():
@@ -656,8 +692,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 +723,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.
 
@@ -993,14 +1029,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 +1044,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 +1133,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 +1157,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 +1183,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 +1254,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 +1263,17 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         return len(self._charpoly_coefficients())
 
 
+    def subalgebra(self, basis, **kwargs):
+        r"""
+        Create a subalgebra of this algebra from the given basis.
+
+        This is a simple wrapper around a subalgebra class constructor
+        that can be overridden in subclasses.
+        """
+        from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
+        return FiniteDimensionalEJASubalgebra(self, basis, **kwargs)
+
+
     def vector_space(self):
         """
         Return the vector space that underlies this algebra.
@@ -2336,7 +2401,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:
@@ -2653,262 +2722,427 @@ class TrivialEJA(ConcreteEJA):
         return cls(**kwargs)
 
 
-class DirectSumEJA(FiniteDimensionalEJA):
+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,
+        ....:                                  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: 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()
-        5
-        sage: J.matrix_space()
-        The Cartesian product of (Full MatrixSpace of 2 by 1 dense matrices
-        over Algebraic Real Field, Full MatrixSpace of 3 by 3 dense matrices
-        over Algebraic Real Field)
+        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:
 
-    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
+
+    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
 
     """
-    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()
-
-        M = J1.matrix_space().cartesian_product(J2.matrix_space())
-        self._cartprod_algebra = J1.cartesian_product(J2)
-
-        self._matrix_basis = tuple( [M((a,0)) for a in J1.matrix_basis()] +
-                                    [M((0,b)) for b in J2.matrix_basis()] )
-
-        n = len(self._matrix_basis)
-        self._sets = None
-        CombinatorialFreeModule.__init__(
-                         self,
-                         field,
-                         range(n),
-                         category=self._cartprod_algebra.category(),
-                         bracket=False,
-                         **kwargs)
-        self.rank.set_cache(J1.rank() + J2.rank())
+    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)
+            )
 
-    def product(self,x,y):
+        # 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 (JordanSpinEJA,
-            ....:                                  ComplexHermitianEJA,
-            ....:                                  DirectSumEJA)
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            ....:                                  RealSymmetricEJA)
 
-        TESTS::
+        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 = JordanSpinEJA(3, field=QQ)
-            sage: J2 = ComplexHermitianEJA(2, field=QQ, orthonormalize=False)
-            sage: J = DirectSumEJA(J1,J2)
-            sage: J.random_element()*J.random_element() in J
+            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
 
         """
-        xv = self._cartprod_algebra.from_vector(x.to_vector())
-        yv = self._cartprod_algebra.from_vector(y.to_vector())
-        return self.from_vector((xv*yv).to_vector())
+        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::
 
-    def cartesian_factors(self):
+            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 _element_constructor_(self, elt):
         r"""
-        Return the pair of this algebra's factors.
+        Construct an element of this algebra from an ordered tuple.
+
+        We just apply the element constructor from each of our factors
+        to the corresponding component of the tuple, and package up
+        the result.
 
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (HadamardEJA,
-            ....:                                  JordanSpinEJA,
-            ....:                                  DirectSumEJA)
+            ....:                                  RealSymmetricEJA)
 
         EXAMPLES::
 
-            sage: J1 = HadamardEJA(2, field=QQ)
-            sage: J2 = JordanSpinEJA(3, field=QQ)
-            sage: J = DirectSumEJA(J1,J2)
-            sage: J.cartesian_factors()
-            (Euclidean Jordan algebra of dimension 2 over Rational Field,
-             Euclidean Jordan algebra of dimension 3 over Rational Field)
-
-        """
-        return self._cartprod_algebra.cartesian_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))
+            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)
+        """
+        m = len(self.cartesian_factors())
+        try:
+            z = tuple( self.cartesian_factors()[i](elt[i]) for i in range(m) )
+            return self._cartesian_product_of_elements(z)
+        except:
+            raise ValueError("not an element of this algebra")
+
+    Element = CartesianProductEJAElement
 
 
+FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA
 
 random_eja = ConcreteEJA.random_instance
+#def random_eja(*args, **kwargs):
+#    from sage.categories.cartesian_product import cartesian_product
+#    J1 = HadamardEJA(1, **kwargs)
+#    J2 = RealSymmetricEJA(2, **kwargs)
+#    J =  cartesian_product([J1,J2])
+#    return J