eja: factor out a class for real-embedded matrices.

[sage.d.git] / mjo / eja / eja_algebra.py

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index a77bb871def9264ac0a82dfea4737b6c2b4411fc..5bf597565b2ae292032c3f0a532763d7889cd2a8 100644 (file)
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -1,9 +1,53 @@
 """
 """

-Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
-specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
-are used in optimization, and have some additional nice methods beyond
-what can be supported in a general Jordan Algebra.
-
+Representations and constructions for Euclidean Jordan algebras.
+
+A Euclidean Jordan algebra is a Jordan algebra that has some
+additional properties:
+
+  1.   It is finite-dimensional.
+  2.   Its scalar field is the real numbers.
+  3a.  An inner product is defined on it, and...
+  3b.  That inner product is compatible with the Jordan product
+       in the sense that `<x*y,z> = <y,x*z>` for all elements
+       `x,y,z` in the algebra.
+
+Every Euclidean Jordan algebra is formally-real: for any two elements
+`x` and `y` in the algebra, `x^{2} + y^{2} = 0` implies that `x = y =
+0`. Conversely, every finite-dimensional formally-real Jordan algebra
+can be made into a Euclidean Jordan algebra with an appropriate choice
+of inner-product.
+
+Formally-real Jordan algebras were originally studied as a framework
+for quantum mechanics. Today, Euclidean Jordan algebras are crucial in
+symmetric cone optimization, since every symmetric cone arises as the
+cone of squares in some Euclidean Jordan algebra.
+
+It is known that every Euclidean Jordan algebra decomposes into an
+orthogonal direct sum (essentially, a Cartesian product) of simple
+algebras, and that moreover, up to Jordan-algebra isomorphism, there
+are only five families of simple algebras. We provide constructions
+for these simple algebras:
+
+  * :class:`BilinearFormEJA`
+  * :class:`RealSymmetricEJA`
+  * :class:`ComplexHermitianEJA`
+  * :class:`QuaternionHermitianEJA`
+
+Missing from this list is the algebra of three-by-three octononion
+Hermitian matrices, as there is (as of yet) no implementation of the
+octonions in SageMath. In addition to these, we provide two other
+example constructions,
+
+  * :class:`HadamardEJA`
+  * :class:`TrivialEJA`
+
+The Jordan spin algebra is a bilinear form algebra where the bilinear
+form is the identity. The Hadamard EJA is simply a Cartesian product
+of one-dimensional spin algebras. And last but not least, the trivial
+EJA is exactly what you think. Cartesian products of these are also
+supported using the usual ``cartesian_product()`` function; as a
+result, we support (up to isomorphism) all Euclidean Jordan algebras
+that don't involve octonions.

 
 SETUP::
 
 
 SETUP::
 


@@ -13,7 +57,6 @@ EXAMPLES::
 
     sage: random_eja()
     Euclidean Jordan algebra of dimension...
 
     sage: random_eja()
     Euclidean Jordan algebra of dimension...

-

 """
 
 from itertools import repeat
 """
 
 from itertools import repeat


@@ -21,8 +64,7 @@ from itertools import repeat
 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
 from sage.categories.magmatic_algebras import MagmaticAlgebras
 from sage.categories.sets_cat import cartesian_product
 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
 from sage.categories.magmatic_algebras import MagmaticAlgebras
 from sage.categories.sets_cat import cartesian_product

-from sage.combinat.free_module import (CombinatorialFreeModule,
-                                       CombinatorialFreeModule_CartesianProduct)
+from sage.combinat.free_module import CombinatorialFreeModule

 from sage.matrix.constructor import matrix
 from sage.matrix.matrix_space import MatrixSpace
 from sage.misc.cachefunc import cached_method
 from sage.matrix.constructor import matrix
 from sage.matrix.matrix_space import MatrixSpace
 from sage.misc.cachefunc import cached_method


@@ -33,7 +75,7 @@ from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
                             QuadraticField)
 from mjo.eja.eja_element import FiniteDimensionalEJAElement
 from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
                             QuadraticField)
 from mjo.eja.eja_element import FiniteDimensionalEJAElement
 from mjo.eja.eja_operator import FiniteDimensionalEJAOperator

-from mjo.eja.eja_utils import _mat2vec
+from mjo.eja.eja_utils import _all2list, _mat2vec

 
 class FiniteDimensionalEJA(CombinatorialFreeModule):
     r"""
 
 class FiniteDimensionalEJA(CombinatorialFreeModule):
     r"""


@@ -41,16 +83,50 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
     INPUT:
 
 
     INPUT:
 

-      - basis -- a tuple of basis elements in their matrix form.
+      - ``basis`` -- a tuple; 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`` -- a function; afunction of two ``basis``
+        elements (in matrix form) that returns their jordan product,
+        also in matrix form; this will be applied to ``basis`` to
+        compute a multiplication table for the algebra.
+
+      - ``inner_product`` -- a function; a function of two ``basis``
+        elements (in matrix form) that returns their inner
+        product. This will be applied to ``basis`` to compute an
+        inner-product table (basically a matrix) for this algebra.
+
+      - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
+        field for the algebra.
+
+      - ``orthonormalize`` -- boolean (default: ``True``); whether or
+        not to orthonormalize the basis. Doing so is expensive and
+        generally rules out using the rationals as your ``field``, but
+        is required for spectral decompositions.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import random_eja

 
 

-      - jordan_product -- function of two elements (in matrix form)
-        that returns their jordan product in this algebra; this will
-        be applied to ``basis`` to compute a multiplication table for
-        the algebra.
+    TESTS:
+
+    We should compute that an element subalgebra is associative even
+    if we circumvent the element method::

 
 

-      - inner_product -- function of two elements (in matrix form) that
-        returns their inner product. This will be applied to ``basis`` to
-        compute an inner-product table (basically a matrix) for this algebra.
+        sage: set_random_seed()
+        sage: J = random_eja(field=QQ,orthonormalize=False)
+        sage: x = J.random_element()
+        sage: A = x.subalgebra_generated_by(orthonormalize=False)
+        sage: basis = tuple(b.superalgebra_element() for b in A.basis())
+        sage: J.subalgebra(basis, orthonormalize=False).is_associative()
+        True

 
     """
     Element = FiniteDimensionalEJAElement
 
     """
     Element = FiniteDimensionalEJAElement


@@ -61,11 +137,13 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                  inner_product,
                  field=AA,
                  orthonormalize=True,
                  inner_product,
                  field=AA,
                  orthonormalize=True,

-                 associative=False,
+                 associative=None,
+                 cartesian_product=False,

                  check_field=True,
                  check_axioms=True,
                  check_field=True,
                  check_axioms=True,

-                 prefix='e',
-                 category=None):
+                 prefix="b"):
+
+        n = len(basis)

 
         if check_field:
             if not field.is_subring(RR):
 
         if check_field:
             if not field.is_subring(RR):


@@ -74,10 +152,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                 # we've specified a real embedding.
                 raise ValueError("scalar field is not real")
 
                 # we've specified a real embedding.
                 raise ValueError("scalar field is not real")
 

-        # 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 check_axioms:
             # Check commutativity of the Jordan and inner-products.
             # This has to be done before we build the multiplication
         if check_axioms:
             # Check commutativity of the Jordan and inner-products.
             # This has to be done before we build the multiplication


@@ -94,33 +168,49 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                 raise ValueError("inner-product is not commutative")
 
 
                 raise ValueError("inner-product is not commutative")
 
 

-        if category is None:
-            category = MagmaticAlgebras(field).FiniteDimensional()
-            category = category.WithBasis().Unital()
-            if associative:
-                # Element subalgebras can take advantage of this.
-                category = category.Associative()
+        category = MagmaticAlgebras(field).FiniteDimensional()
+        category = category.WithBasis().Unital().Commutative()
+
+        if associative is None:
+            # We should figure it out. As with check_axioms, we have to do
+            # this without the help of the _jordan_product_is_associative()
+            # method because we need to know the category before we
+            # initialize the algebra.
+            associative = all( jordan_product(jordan_product(bi,bj),bk)
+                               ==
+                               jordan_product(bi,jordan_product(bj,bk))
+                               for bi in basis
+                               for bj in basis
+                               for bk in basis)
+
+        if associative:
+            # Element subalgebras can take advantage of this.
+            category = category.Associative()
+        if cartesian_product:
+            # Use join() here because otherwise we only get the
+            # "Cartesian product of..." and not the things themselves.
+            category = category.join([category,
+                                      category.CartesianProducts()])

 
         # Call the superclass constructor so that we can use its from_vector()
         # method to build our multiplication table.
 
         # Call the superclass constructor so that we can use its from_vector()
         # method to build our multiplication table.

-        n = len(basis)
-        super().__init__(field,
-                         range(n),
-                         prefix=prefix,
-                         category=category,
-                         bracket=False)
+        CombinatorialFreeModule.__init__(self,
+                                         field,
+                                         range(n),
+                                         prefix=prefix,
+                                         category=category,
+                                         bracket=False)

 
         # Now comes all of the hard work. We'll be constructing an
         # ambient vector space V that our (vectorized) basis lives in,
         # as well as a subspace W of V spanned by those (vectorized)
         # basis elements. The W-coordinates are the coefficients that
 
         # Now comes all of the hard work. We'll be constructing an
         # ambient vector space V that our (vectorized) basis lives in,
         # as well as a subspace W of V spanned by those (vectorized)
         # basis elements. The W-coordinates are the coefficients that

-        # we see in things like x = 1*e1 + 2*e2.
+        # we see in things like x = 1*b1 + 2*b2.

         vector_basis = basis
 
         degree = 0
         if n > 0:
         vector_basis = basis
 
         degree = 0
         if n > 0:

-            # Works on both column and square matrices...
-            degree = len(basis[0].list())
+            degree = len(_all2list(basis[0]))

 
         # Build an ambient space that fits our matrix basis when
         # written out as "long vectors."
 
         # Build an ambient space that fits our matrix basis when
         # written out as "long vectors."


@@ -134,7 +224,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.
             # Save a copy of the un-orthonormalized basis for later.
             # Convert it to ambient V (vector) coordinates while we're
             # at it, because we'd have to do it later anyway.

-            deortho_vector_basis = tuple( V(b.list()) for b in basis )
+            deortho_vector_basis = tuple( V(_all2list(b)) for b in basis )

 
             from mjo.eja.eja_utils import gram_schmidt
             basis = tuple(gram_schmidt(basis, inner_product))
 
             from mjo.eja.eja_utils import gram_schmidt
             basis = tuple(gram_schmidt(basis, inner_product))


@@ -146,7 +236,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).
         # Now create the vector space for the algebra, which will have
         # its own set of non-ambient coordinates (in terms of the
         # supplied basis).

-        vector_basis = tuple( V(b.list()) for b in basis )
+        vector_basis = tuple( V(_all2list(b)) for b in basis )

         W = V.span_of_basis( vector_basis, check=check_axioms)
 
         if orthonormalize:
         W = V.span_of_basis( vector_basis, check=check_axioms)
 
         if orthonormalize:


@@ -178,7 +268,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)
                 # The jordan product returns a matrixy answer, so we
                 # have to convert it to the algebra coordinates.
                 elt = jordan_product(q_i, q_j)

-                elt = W.coordinate_vector(V(elt.list()))
+                elt = W.coordinate_vector(V(_all2list(elt)))

                 self._multiplication_table[i][j] = self.from_vector(elt)
 
                 if not orthonormalize:
                 self._multiplication_table[i][j] = self.from_vector(elt)
 
                 if not orthonormalize:


@@ -226,6 +316,35 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
 
     def product_on_basis(self, i, j):
 
 
     def product_on_basis(self, i, j):

+        r"""
+        Returns the Jordan product of the `i` and `j`th basis elements.
+
+        This completely defines the Jordan product on the algebra, and
+        is used direclty by our superclass machinery to implement
+        :meth:`product`.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import random_eja
+
+        TESTS::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: n = J.dimension()
+            sage: bi = J.zero()
+            sage: bj = J.zero()
+            sage: bi_bj = J.zero()*J.zero()
+            sage: if n > 0:
+            ....:     i = ZZ.random_element(n)
+            ....:     j = ZZ.random_element(n)
+            ....:     bi = J.monomial(i)
+            ....:     bj = J.monomial(j)
+            ....:     bi_bj = J.product_on_basis(i,j)
+            sage: bi*bj == bi_bj
+            True
+
+        """

         # We only stored the lower-triangular portion of the
         # multiplication table.
         if j <= i:
         # We only stored the lower-triangular portion of the
         # multiplication table.
         if j <= i:


@@ -283,11 +402,33 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: y = J.random_element()
             sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
             True
             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())
 
 
         """
         B = self._inner_product_matrix
         return (B*x.to_vector()).inner_product(y.to_vector())
 
 

+    def is_associative(self):
+        r"""
+        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
+
+        """
+        return "Associative" in self.category().axioms()
+

     def _is_commutative(self):
         r"""
         Whether or not this algebra's multiplication table is commutative.
     def _is_commutative(self):
         r"""
         Whether or not this algebra's multiplication table is commutative.


@@ -296,9 +437,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         this algebra was constructed with ``check_axioms=False`` and
         passed an invalid multiplication table.
         """
         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 all( x*y == y*x for x in self.gens() for y in self.gens() )

 
     def _is_jordanian(self):
         r"""
 
     def _is_jordanian(self):
         r"""


@@ -317,6 +456,92 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                     for i in range(self.dimension())
                     for j in range(self.dimension()) )
 
                     for i in range(self.dimension())
                     for j in range(self.dimension()) )
 

+    def _jordan_product_is_associative(self):
+        r"""
+        Return whether or not this algebra's Jordan product is
+        associative; that is, whether or not `x*(y*z) = (x*y)*z`
+        for all `x,y,x`.
+
+        This method should agree with :meth:`is_associative` unless
+        you lied about the value of the ``associative`` parameter
+        when you constructed the algebra.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (random_eja,
+            ....:                                  RealSymmetricEJA,
+            ....:                                  ComplexHermitianEJA,
+            ....:                                  QuaternionHermitianEJA)
+
+        EXAMPLES::
+
+            sage: J = RealSymmetricEJA(4, orthonormalize=False)
+            sage: J._jordan_product_is_associative()
+            False
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by()
+            sage: A._jordan_product_is_associative()
+            True
+
+        ::
+
+            sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: J._jordan_product_is_associative()
+            False
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
+            sage: A._jordan_product_is_associative()
+            True
+
+        ::
+
+            sage: J = QuaternionHermitianEJA(2)
+            sage: J._jordan_product_is_associative()
+            False
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by()
+            sage: A._jordan_product_is_associative()
+            True
+
+        TESTS:
+
+        The values we've presupplied to the constructors agree with
+        the computation::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: J.is_associative() == J._jordan_product_is_associative()
+            True
+
+        """
+        R = self.base_ring()
+
+        # Used to check whether or not something is zero.
+        epsilon = R.zero()
+        if not R.is_exact():
+            # I don't know of any examples that make this magnitude
+            # necessary because I don't know how to make an
+            # associative algebra when the element subalgebra
+            # construction is unreliable (as it is over RDF; we can't
+            # find the degree of an element because we can't compute
+            # the rank of a matrix). But even multiplication of floats
+            # is non-associative, so *some* epsilon is needed... let's
+            # just take the one from _inner_product_is_associative?
+            epsilon = 1e-15
+
+        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)
+                    diff = (x*y)*z - x*(y*z)
+
+                    if diff.norm() > epsilon:
+                        return False
+
+        return True
+

     def _inner_product_is_associative(self):
         r"""
         Return whether or not this algebra's inner product `B` is
     def _inner_product_is_associative(self):
         r"""
         Return whether or not this algebra's inner product `B` is


@@ -326,11 +551,14 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         this algebra was constructed with ``check_axioms=False`` and
         passed an invalid Jordan or inner-product.
         """
         this algebra was constructed with ``check_axioms=False`` and
         passed an invalid Jordan or inner-product.
         """

+        R = self.base_ring()

 
 

-        # Used to check whether or not something is zero in an inexact
-        # ring. This number is sufficient to allow the construction of
-        # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
-        epsilon = 1e-16
+        # Used to check whether or not something is zero.
+        epsilon = R.zero()
+        if not R.is_exact():
+            # This choice is sufficient to allow the construction of
+            # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
+            epsilon = 1e-15

 
         for i in range(self.dimension()):
             for j in range(self.dimension()):
 
         for i in range(self.dimension()):
             for j in range(self.dimension()):


@@ -340,12 +568,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                     z = self.monomial(k)
                     diff = (x*y).inner_product(z) - x.inner_product(y*z)
 
                     z = self.monomial(k)
                     diff = (x*y).inner_product(z) - x.inner_product(y*z)
 

-                    if self.base_ring().is_exact():
-                        if diff != 0:
-                            return False
-                    else:
-                        if diff.abs() > epsilon:
-                            return False
+                    if diff.abs() > epsilon:
+                        return False

 
         return True
 
 
         return True
 


@@ -359,7 +583,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         SETUP::
 
 
         SETUP::
 

-            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            sage: from mjo.eja.eja_algebra import (random_eja,
+            ....:                                  JordanSpinEJA,

             ....:                                  HadamardEJA,
             ....:                                  RealSymmetricEJA)
 
             ....:                                  HadamardEJA,
             ....:                                  RealSymmetricEJA)
 


@@ -381,29 +606,42 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             ...
             ValueError: not an element of this algebra
 
             ...
             ValueError: not an element of this algebra
 

+        Tuples work as well, provided that the matrix basis for the
+        algebra consists of them::
+
+            sage: J1 = HadamardEJA(3)
+            sage: J2 = RealSymmetricEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
+            b1 + b5
+

         TESTS:
 
         TESTS:
 

-        Ensure that we can convert any element of the two non-matrix
-        simple algebras (whose matrix representations are columns)
-        back and forth faithfully::
+        Ensure that we can convert any element back and forth
+        faithfully between its matrix and algebra representations::

 
             sage: set_random_seed()
 
             sage: set_random_seed()

-            sage: J = HadamardEJA.random_instance()
-            sage: x = J.random_element()
-            sage: J(x.to_vector().column()) == x
-            True
-            sage: J = JordanSpinEJA.random_instance()
+            sage: J = random_eja()

             sage: x = J.random_element()
             sage: x = J.random_element()

-            sage: J(x.to_vector().column()) == x
+            sage: J(x.to_matrix()) == x

             True
 
             True
 

+        We cannot coerce elements between algebras just because their
+        matrix representations are compatible::
+
+            sage: J1 = HadamardEJA(3)
+            sage: J2 = JordanSpinEJA(3)
+            sage: J2(J1.one())
+            Traceback (most recent call last):
+            ...
+            ValueError: not an element of this algebra
+            sage: J1(J2.zero())
+            Traceback (most recent call last):
+            ...
+            ValueError: not an element of this algebra

         """
         msg = "not an element of this algebra"
         """
         msg = "not an element of this algebra"

-        if elt == 0:
-            # The superclass implementation of random_element()
-            # needs to be able to coerce "0" into the algebra.
-            return self.zero()
-        elif elt in self.base_ring():
+        if elt in self.base_ring():

             # Ensure that no base ring -> algebra coercion is performed
             # by this method. There's some stupidity in sage that would
             # otherwise propagate to this method; for example, sage thinks
             # Ensure that no base ring -> algebra coercion is performed
             # by this method. There's some stupidity in sage that would
             # otherwise propagate to this method; for example, sage thinks


@@ -411,9 +649,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             raise ValueError(msg)
 
         try:
             raise ValueError(msg)
 
         try:

+            # Try to convert a vector into a column-matrix...

             elt = elt.column()
         except (AttributeError, TypeError):
             elt = elt.column()
         except (AttributeError, TypeError):

-            # Try to convert a vector into a column-matrix
+            # and ignore failure, because we weren't really expecting
+            # a vector as an argument anyway.

             pass
 
         if elt not in self.matrix_space():
             pass
 
         if elt not in self.matrix_space():


@@ -426,14 +666,20 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # closure whereas the base ring of the 3-by-3 identity matrix
         # could be QQ instead of QQbar.
         #
         # closure whereas the base ring of the 3-by-3 identity matrix
         # could be QQ instead of QQbar.
         #

+        # And, we also have to handle Cartesian product bases (when
+        # the matrix basis consists of tuples) here. The "good news"
+        # is that we're already converting everything to long vectors,
+        # and that strategy works for tuples as well.
+        #

         # We pass check=False because the matrix basis is "guaranteed"
         # to be linearly independent... right? Ha ha.
         # We pass check=False because the matrix basis is "guaranteed"
         # to be linearly independent... right? Ha ha.

-        V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols())
-        W = V.span_of_basis( (_mat2vec(s) for s in self.matrix_basis()),
+        elt = _all2list(elt)
+        V = VectorSpace(self.base_ring(), len(elt))
+        W = V.span_of_basis( (V(_all2list(s)) for s in self.matrix_basis()),

                              check=False)
 
         try:
                              check=False)
 
         try:

-            coords =  W.coordinate_vector(_mat2vec(elt))
+            coords = W.coordinate_vector(V(elt))

         except ArithmeticError:  # vector is not in free module
             raise ValueError(msg)
 
         except ArithmeticError:  # vector is not in free module
             raise ValueError(msg)
 


@@ -642,15 +888,15 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: J = JordanSpinEJA(4)
             sage: J.multiplication_table()
             +----++----+----+----+----+
             sage: J = JordanSpinEJA(4)
             sage: J.multiplication_table()
             +----++----+----+----+----+

-            | *  || e0 | e1 | e2 | e3 |
+            | *  || b0 | b1 | b2 | b3 |

             +====++====+====+====+====+
             +====++====+====+====+====+

-            | e0 || e0 | e1 | e2 | e3 |
+            | b0 || b0 | b1 | b2 | b3 |

             +----++----+----+----+----+
             +----++----+----+----+----+

-            | e1 || e1 | e0 | 0  | 0  |
+            | b1 || b1 | b0 | 0  | 0  |

             +----++----+----+----+----+
             +----++----+----+----+----+

-            | e2 || e2 | 0  | e0 | 0  |
+            | b2 || b2 | 0  | b0 | 0  |

             +----++----+----+----+----+
             +----++----+----+----+----+

-            | e3 || e3 | 0  | 0  | e0 |
+            | b3 || b3 | 0  | 0  | b0 |

             +----++----+----+----+----+
 
         """
             +----++----+----+----+----+
 
         """


@@ -660,8 +906,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         # And to each subsequent row, prepend an entry that belongs to
         # the left-side "header column."
 
         # 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.monomial(i)] + [ self.monomial(i)*self.monomial(j)
+                                    for j in range(n) ]

                for i in range(n) ]
 
         return table(M, header_row=True, header_column=True, frame=True)
                for i in range(n) ]
 
         return table(M, header_row=True, header_column=True, frame=True)


@@ -704,7 +950,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: J = RealSymmetricEJA(2)
             sage: J.basis()
 
             sage: J = RealSymmetricEJA(2)
             sage: J.basis()

-            Finite family {0: e0, 1: e1, 2: e2}
+            Finite family {0: b0, 1: b1, 2: b2}

             sage: J.matrix_basis()
             (
             [1 0]  [                  0 0.7071067811865475?]  [0 0]
             sage: J.matrix_basis()
             (
             [1 0]  [                  0 0.7071067811865475?]  [0 0]


@@ -715,7 +961,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: J = JordanSpinEJA(2)
             sage: J.basis()
 
             sage: J = JordanSpinEJA(2)
             sage: J.basis()

-            Finite family {0: e0, 1: e1}
+            Finite family {0: b0, 1: b1}

             sage: J.matrix_basis()
             (
             [1]  [0]
             sage: J.matrix_basis()
             (
             [1]  [0]


@@ -731,12 +977,49 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         we think of them as matrices (including column vectors of the
         appropriate size).
 
         we think of them as matrices (including column vectors of the
         appropriate size).
 

-        Generally this will be an `n`-by-`1` column-vector space,
+        "By default" this will be an `n`-by-`1` column-matrix space,

         except when the algebra is trivial. There it's `n`-by-`n`
         (where `n` is zero), to ensure that two elements of the matrix
         except when the algebra is trivial. There it's `n`-by-`n`
         (where `n` is zero), to ensure that two elements of the matrix

-        space (empty matrices) can be multiplied.
+        space (empty matrices) can be multiplied. For algebras of
+        matrices, this returns the space in which their
+        real embeddings live.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
+            ....:                                  JordanSpinEJA,
+            ....:                                  QuaternionHermitianEJA,
+            ....:                                  TrivialEJA)
+
+        EXAMPLES:
+
+        By default, the matrix representation is just a column-matrix
+        equivalent to the vector representation::
+
+            sage: J = JordanSpinEJA(3)
+            sage: J.matrix_space()
+            Full MatrixSpace of 3 by 1 dense matrices over Algebraic
+            Real Field
+
+        The matrix representation in the trivial algebra is
+        zero-by-zero instead of the usual `n`-by-one::
+
+            sage: J = TrivialEJA()
+            sage: J.matrix_space()
+            Full MatrixSpace of 0 by 0 dense matrices over Algebraic
+            Real Field
+
+        The matrix space for complex/quaternion Hermitian matrix EJA
+        is the space in which their real-embeddings live, not the
+        original complex/quaternion matrix space::
+
+            sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: J.matrix_space()
+            Full MatrixSpace of 4 by 4 dense matrices over Rational Field
+            sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
+            sage: J.matrix_space()
+            Full MatrixSpace of 4 by 4 dense matrices over Rational Field

 
 

-        Matrix algebras override this with something more useful.

         """
         if self.is_trivial():
             return MatrixSpace(self.base_ring(), 0)
         """
         if self.is_trivial():
             return MatrixSpace(self.base_ring(), 0)


@@ -760,20 +1043,20 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: J = HadamardEJA(5)
             sage: J.one()
 
             sage: J = HadamardEJA(5)
             sage: J.one()

-            e0 + e1 + e2 + e3 + e4
+            b0 + b1 + b2 + b3 + b4

 
         The unit element in the Hadamard EJA is inherited in the
         subalgebras generated by its elements::
 
             sage: J = HadamardEJA(5)
             sage: J.one()
 
         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
+            b0 + b1 + b2 + b3 + b4

             sage: x = sum(J.gens())
             sage: A = x.subalgebra_generated_by(orthonormalize=False)
             sage: A.one()
             sage: x = sum(J.gens())
             sage: A = x.subalgebra_generated_by(orthonormalize=False)
             sage: A.one()

-            f0
+            c0

             sage: A.one().superalgebra_element()
             sage: A.one().superalgebra_element()

-            e0 + e1 + e2 + e3 + e4
+            b0 + b1 + b2 + b3 + b4

 
         TESTS:
 
 
         TESTS:
 


@@ -997,14 +1280,12 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         if not c.is_idempotent():
             raise ValueError("element is not idempotent: %s" % c)
 
         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).
         # 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
         J0 = trivial                          # eigenvalue zero
         J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
         J1 = trivial                          # eigenvalue one


@@ -1014,9 +1295,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                 J5 = eigspace
             else:
                 gens = tuple( self.from_vector(b) for b in eigspace.basis() )
                 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:
                 if eigval == 0:
                     J0 = subalg
                 elif eigval == 1:


@@ -1235,6 +1514,14 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         return len(self._charpoly_coefficients())
 
 
         return len(self._charpoly_coefficients())
 
 

+    def subalgebra(self, basis, **kwargs):
+        r"""
+        Create a subalgebra of this algebra from the given basis.
+        """
+        from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
+        return FiniteDimensionalEJASubalgebra(self, basis, **kwargs)
+
+

     def vector_space(self):
         """
         Return the vector space that underlies this algebra.
     def vector_space(self):
         """
         Return the vector space that underlies this algebra.


@@ -1253,7 +1540,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         return self.zero().to_vector().parent().ambient_vector_space()
 
 
         return self.zero().to_vector().parent().ambient_vector_space()
 
 

-    Element = FiniteDimensionalEJAElement

 
 class RationalBasisEJA(FiniteDimensionalEJA):
     r"""
 
 class RationalBasisEJA(FiniteDimensionalEJA):
     r"""


@@ -1291,6 +1577,13 @@ 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")
 
             if not all( all(b_i in QQ for b_i in b.list()) for b in basis ):
                 raise TypeError("basis not rational")
 

+        super().__init__(basis,
+                         jordan_product,
+                         inner_product,
+                         field=field,
+                         check_field=check_field,
+                         **kwargs)
+

         self._rational_algebra = None
         if field is not QQ:
             # There's no point in constructing the extra algebra if this
         self._rational_algebra = None
         if field is not QQ:
             # There's no point in constructing the extra algebra if this


@@ -1304,17 +1597,11 @@ class RationalBasisEJA(FiniteDimensionalEJA):
                                        jordan_product,
                                        inner_product,
                                        field=QQ,
                                        jordan_product,
                                        inner_product,
                                        field=QQ,

+                                       associative=self.is_associative(),

                                        orthonormalize=False,
                                        check_field=False,
                                        check_axioms=False)
 
                                        orthonormalize=False,
                                        check_field=False,
                                        check_axioms=False)
 

-        super().__init__(basis,
-                         jordan_product,
-                         inner_product,
-                         field=field,
-                         check_field=check_field,
-                         **kwargs)
-

     @cached_method
     def _charpoly_coefficients(self):
         r"""
     @cached_method
     def _charpoly_coefficients(self):
         r"""


@@ -1438,6 +1725,21 @@ class ConcreteEJA(RationalBasisEJA):
 
 
 class MatrixEJA:
 
 
 class MatrixEJA:

+    @staticmethod
+    def jordan_product(X,Y):
+        return (X*Y + Y*X)/2
+
+    @staticmethod
+    def trace_inner_product(X,Y):
+        r"""
+        A trace inner-product for matrices that aren't embedded in the
+        reals.
+        """
+        # We take the norm (absolute value) because Octonions() isn't
+        # smart enough yet to coerce its one() into the base field.
+        return (X*Y).trace().abs()
+
+class RealEmbeddedMatrixEJA(MatrixEJA):

     @staticmethod
     def dimension_over_reals():
         r"""
     @staticmethod
     def dimension_over_reals():
         r"""


@@ -1483,9 +1785,6 @@ class MatrixEJA:
             raise ValueError("the matrix 'M' must be a real embedding")
         return M
 
             raise ValueError("the matrix 'M' must be a real embedding")
         return M
 

-    @staticmethod
-    def jordan_product(X,Y):
-        return (X*Y + Y*X)/2

 
     @classmethod
     def trace_inner_product(cls,X,Y):
 
     @classmethod
     def trace_inner_product(cls,X,Y):


@@ -1494,29 +1793,11 @@ class MatrixEJA:
 
         SETUP::
 
 
         SETUP::
 

-            sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
-            ....:                                  ComplexHermitianEJA,
+            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,

             ....:                                  QuaternionHermitianEJA)
 
         EXAMPLES::
 
             ....:                                  QuaternionHermitianEJA)
 
         EXAMPLES::
 

-        This gives the same answer as it would if we computed the trace
-        from the unembedded (original) matrices::
-
-            sage: set_random_seed()
-            sage: J = RealSymmetricEJA.random_instance()
-            sage: x,y = J.random_elements(2)
-            sage: Xe = x.to_matrix()
-            sage: Ye = y.to_matrix()
-            sage: X = J.real_unembed(Xe)
-            sage: Y = J.real_unembed(Ye)
-            sage: expected = (X*Y).trace()
-            sage: actual = J.trace_inner_product(Xe,Ye)
-            sage: actual == expected
-            True
-
-        ::
-

             sage: set_random_seed()
             sage: J = ComplexHermitianEJA.random_instance()
             sage: x,y = J.random_elements(2)
             sage: set_random_seed()
             sage: J = ComplexHermitianEJA.random_instance()
             sage: x,y = J.random_elements(2)


@@ -1544,27 +1825,15 @@ class MatrixEJA:
             True
 
         """
             True
 
         """

-        Xu = cls.real_unembed(X)
-        Yu = cls.real_unembed(Y)
-        tr = (Xu*Yu).trace()
-
-        try:
-            # Works in QQ, AA, RDF, et cetera.
-            return tr.real()
-        except AttributeError:
-            # A quaternion doesn't have a real() method, but does
-            # have coefficient_tuple() method that returns the
-            # coefficients of 1, i, j, and k -- in that order.
-            return tr.coefficient_tuple()[0]
+        # This does in fact compute the real part of the trace.
+        # If we compute the trace of e.g. a complex matrix M,
+        # then we do so by adding up its diagonal entries --
+        # call them z_1 through z_n. The real embedding of z_1
+        # will be a 2-by-2 REAL matrix [a, b; -b, a] whose trace
+        # as a REAL matrix will be 2*a = 2*Re(z_1). And so forth.
+        return (X*Y).trace()/cls.dimension_over_reals()

 
 

-
-class RealMatrixEJA(MatrixEJA):
-    @staticmethod
-    def dimension_over_reals():
-        return 1
-
-
-class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
+class RealSymmetricEJA(ConcreteEJA, MatrixEJA):

     """
     The rank-n simple EJA consisting of real symmetric n-by-n
     matrices, the usual symmetric Jordan product, and the trace inner
     """
     The rank-n simple EJA consisting of real symmetric n-by-n
     matrices, the usual symmetric Jordan product, and the trace inner


@@ -1577,19 +1846,19 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
     EXAMPLES::
 
         sage: J = RealSymmetricEJA(2)
     EXAMPLES::
 
         sage: J = RealSymmetricEJA(2)

-        sage: e0, e1, e2 = J.gens()
-        sage: e0*e0
-        e0
-        sage: e1*e1
-        1/2*e0 + 1/2*e2
-        sage: e2*e2
-        e2
+        sage: b0, b1, b2 = J.gens()
+        sage: b0*b0
+        b0
+        sage: b1*b1
+        1/2*b0 + 1/2*b2
+        sage: b2*b2
+        b2

 
     In theory, our "field" can be any subfield of the reals::
 
 
     In theory, our "field" can be any subfield of the reals::
 

-        sage: RealSymmetricEJA(2, field=RDF)
+        sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)

         Euclidean Jordan algebra of dimension 3 over Real Double Field
         Euclidean Jordan algebra of dimension 3 over Real Double Field

-        sage: RealSymmetricEJA(2, field=RR)
+        sage: RealSymmetricEJA(2, field=RR, check_axioms=True)

         Euclidean Jordan algebra of dimension 3 over Real Field with
         53 bits of precision
 
         Euclidean Jordan algebra of dimension 3 over Real Field with
         53 bits of precision
 


@@ -1630,7 +1899,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
 
     """
     @classmethod
 
     """
     @classmethod

-    def _denormalized_basis(cls, n):
+    def _denormalized_basis(cls, n, field):

         """
         Return a basis for the space of real symmetric n-by-n matrices.
 
         """
         Return a basis for the space of real symmetric n-by-n matrices.
 


@@ -1642,7 +1911,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)

-            sage: B = RealSymmetricEJA._denormalized_basis(n)
+            sage: B = RealSymmetricEJA._denormalized_basis(n,ZZ)

             sage: all( M.is_symmetric() for M in  B)
             True
 
             sage: all( M.is_symmetric() for M in  B)
             True
 


@@ -1652,7 +1921,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
         S = []
         for i in range(n):
             for j in range(i+1):
         S = []
         for i in range(n):
             for j in range(i+1):

-                Eij = matrix(ZZ, n, lambda k,l: k==i and l==j)
+                Eij = matrix(field, n, lambda k,l: k==i and l==j)

                 if i == j:
                     Sij = Eij
                 else:
                 if i == j:
                     Sij = Eij
                 else:


@@ -1673,26 +1942,32 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
         return cls(n, **kwargs)
 
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
         return cls(n, **kwargs)
 

-    def __init__(self, n, **kwargs):
+    def __init__(self, n, field=AA, **kwargs):

         # 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
 
         # 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
 

-        super(RealSymmetricEJA, self).__init__(self._denormalized_basis(n),
-                                               self.jordan_product,
-                                               self.trace_inner_product,
-                                               **kwargs)
+        associative = False
+        if n <= 1:
+            associative = True
+
+        super().__init__(self._denormalized_basis(n,field),
+                         self.jordan_product,
+                         self.trace_inner_product,
+                         field=field,
+                         associative=associative,
+                         **kwargs)

 
         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         self.rank.set_cache(n)
 
         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         self.rank.set_cache(n)

-        idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
+        idV = self.matrix_space().one()

         self.one.set_cache(self(idV))
 
 
 
         self.one.set_cache(self(idV))
 
 
 

-class ComplexMatrixEJA(MatrixEJA):
+class ComplexMatrixEJA(RealEmbeddedMatrixEJA):

     # A manual dictionary-cache for the complex_extension() method,
     # since apparently @classmethods can't also be @cached_methods.
     _complex_extension = {}
     # A manual dictionary-cache for the complex_extension() method,
     # since apparently @classmethods can't also be @cached_methods.
     _complex_extension = {}


@@ -1771,7 +2046,7 @@ class ComplexMatrixEJA(MatrixEJA):
             True
 
         """
             True
 
         """

-        super(ComplexMatrixEJA,cls).real_embed(M)
+        super().real_embed(M)

         n = M.nrows()
 
         # We don't need any adjoined elements...
         n = M.nrows()
 
         # We don't need any adjoined elements...


@@ -1818,7 +2093,7 @@ class ComplexMatrixEJA(MatrixEJA):
             True
 
         """
             True
 
         """

-        super(ComplexMatrixEJA,cls).real_unembed(M)
+        super().real_unembed(M)

         n = ZZ(M.nrows())
         d = cls.dimension_over_reals()
         F = cls.complex_extension(M.base_ring())
         n = ZZ(M.nrows())
         d = cls.dimension_over_reals()
         F = cls.complex_extension(M.base_ring())


@@ -1855,9 +2130,9 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
 
     In theory, our "field" can be any subfield of the reals::
 
 
     In theory, our "field" can be any subfield of the reals::
 

-        sage: ComplexHermitianEJA(2, field=RDF)
+        sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True)

         Euclidean Jordan algebra of dimension 4 over Real Double Field
         Euclidean Jordan algebra of dimension 4 over Real Double Field

-        sage: ComplexHermitianEJA(2, field=RR)
+        sage: ComplexHermitianEJA(2, field=RR, check_axioms=True)

         Euclidean Jordan algebra of dimension 4 over Real Field with
         53 bits of precision
 
         Euclidean Jordan algebra of dimension 4 over Real Field with
         53 bits of precision
 


@@ -1899,7 +2174,7 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
     """
 
     @classmethod
     """
 
     @classmethod

-    def _denormalized_basis(cls, n):
+    def _denormalized_basis(cls, n, field):

         """
         Returns a basis for the space of complex Hermitian n-by-n matrices.
 
         """
         Returns a basis for the space of complex Hermitian n-by-n matrices.
 


@@ -1917,15 +2192,14 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)

-            sage: B = ComplexHermitianEJA._denormalized_basis(n)
+            sage: B = ComplexHermitianEJA._denormalized_basis(n,ZZ)

             sage: all( M.is_symmetric() for M in  B)
             True
 
         """
             sage: all( M.is_symmetric() for M in  B)
             True
 
         """

-        field = ZZ
-        R = PolynomialRing(field, 'z')
+        R = PolynomialRing(ZZ, 'z')

         z = R.gen()
         z = R.gen()

-        F = field.extension(z**2 + 1, 'I')
+        F = ZZ.extension(z**2 + 1, 'I')

         I = F.gen(1)
 
         # This is like the symmetric case, but we need to be careful:
         I = F.gen(1)
 
         # This is like the symmetric case, but we need to be careful:


@@ -1956,20 +2230,26 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
                 # "erase" E_ij
                 Eij[i,j] = 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".
+        # Since we embedded the entries, we can drop back to the
+        # desired real "field" instead of the extension "F".

         return tuple( s.change_ring(field) for s in S )
 
 
         return tuple( s.change_ring(field) for s in S )
 
 

-    def __init__(self, n, **kwargs):
+    def __init__(self, n, field=AA, **kwargs):

         # 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
 
         # 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
 

-        super(ComplexHermitianEJA, self).__init__(self._denormalized_basis(n),
-                                                  self.jordan_product,
-                                                  self.trace_inner_product,
-                                                  **kwargs)
+        associative = False
+        if n <= 1:
+            associative = True
+
+        super().__init__(self._denormalized_basis(n,field),
+                         self.jordan_product,
+                         self.trace_inner_product,
+                         field=field,
+                         associative=associative,
+                         **kwargs)

         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().


@@ -1989,7 +2269,7 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
         return cls(n, **kwargs)
 
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
         return cls(n, **kwargs)
 

-class QuaternionMatrixEJA(MatrixEJA):
+class QuaternionMatrixEJA(RealEmbeddedMatrixEJA):

 
     # A manual dictionary-cache for the quaternion_extension() method,
     # since apparently @classmethods can't also be @cached_methods.
 
     # A manual dictionary-cache for the quaternion_extension() method,
     # since apparently @classmethods can't also be @cached_methods.


@@ -2052,7 +2332,7 @@ class QuaternionMatrixEJA(MatrixEJA):
             True
 
         """
             True
 
         """

-        super(QuaternionMatrixEJA,cls).real_embed(M)
+        super().real_embed(M)

         quaternions = M.base_ring()
         n = M.nrows()
 
         quaternions = M.base_ring()
         n = M.nrows()
 


@@ -2107,7 +2387,7 @@ class QuaternionMatrixEJA(MatrixEJA):
             True
 
         """
             True
 
         """

-        super(QuaternionMatrixEJA,cls).real_unembed(M)
+        super().real_unembed(M)

         n = ZZ(M.nrows())
         d = cls.dimension_over_reals()
 
         n = ZZ(M.nrows())
         d = cls.dimension_over_reals()
 


@@ -2152,9 +2432,9 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
 
     In theory, our "field" can be any subfield of the reals::
 
 
     In theory, our "field" can be any subfield of the reals::
 

-        sage: QuaternionHermitianEJA(2, field=RDF)
+        sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)

         Euclidean Jordan algebra of dimension 6 over Real Double Field
         Euclidean Jordan algebra of dimension 6 over Real Double Field

-        sage: QuaternionHermitianEJA(2, field=RR)
+        sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)

         Euclidean Jordan algebra of dimension 6 over Real Field with
         53 bits of precision
 
         Euclidean Jordan algebra of dimension 6 over Real Field with
         53 bits of precision
 


@@ -2195,7 +2475,7 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
 
     """
     @classmethod
 
     """
     @classmethod

-    def _denormalized_basis(cls, n):
+    def _denormalized_basis(cls, n, field):

         """
         Returns a basis for the space of quaternion Hermitian n-by-n matrices.
 
         """
         Returns a basis for the space of quaternion Hermitian n-by-n matrices.
 


@@ -2213,12 +2493,11 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)

-            sage: B = QuaternionHermitianEJA._denormalized_basis(n)
+            sage: B = QuaternionHermitianEJA._denormalized_basis(n,ZZ)

             sage: all( M.is_symmetric() for M in B )
             True
 
         """
             sage: all( M.is_symmetric() for M in B )
             True
 
         """

-        field = ZZ

         Q = QuaternionAlgebra(QQ,-1,-1)
         I,J,K = Q.gens()
 
         Q = QuaternionAlgebra(QQ,-1,-1)
         I,J,K = Q.gens()
 


@@ -2262,20 +2541,27 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
                 # "erase" E_ij
                 Eij[i,j] = 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".
+        # Since we embedded the entries, we can drop back to the
+        # desired real "field" instead of the quaternion algebra "Q".

         return tuple( s.change_ring(field) for s in S )
 
 
         return tuple( s.change_ring(field) for s in S )
 
 

-    def __init__(self, n, **kwargs):
+    def __init__(self, n, field=AA, **kwargs):

         # 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
 
         # 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
 

-        super(QuaternionHermitianEJA, self).__init__(self._denormalized_basis(n),
-                                                     self.jordan_product,
-                                                     self.trace_inner_product,
-                                                     **kwargs)
+        associative = False
+        if n <= 1:
+            associative = True
+
+        super().__init__(self._denormalized_basis(n,field),
+                         self.jordan_product,
+                         self.trace_inner_product,
+                         field=field,
+                         associative=associative,
+                         **kwargs)
+

         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().


@@ -2318,19 +2604,19 @@ class HadamardEJA(ConcreteEJA):
     This multiplication table can be verified by hand::
 
         sage: J = HadamardEJA(3)
     This multiplication table can be verified by hand::
 
         sage: J = HadamardEJA(3)

-        sage: e0,e1,e2 = J.gens()
-        sage: e0*e0
-        e0
-        sage: e0*e1
+        sage: b0,b1,b2 = J.gens()
+        sage: b0*b0
+        b0
+        sage: b0*b1

         0
         0

-        sage: e0*e2
+        sage: b0*b2

         0
         0

-        sage: e1*e1
-        e1
-        sage: e1*e2
+        sage: b1*b1
+        b1
+        sage: b1*b2

         0
         0

-        sage: e2*e2
-        e2
+        sage: b2*b2
+        b2

 
     TESTS:
 
 
     TESTS:
 


@@ -2340,7 +2626,7 @@ class HadamardEJA(ConcreteEJA):
         (r0, r1, r2)
 
     """
         (r0, r1, r2)
 
     """

-    def __init__(self, n, **kwargs):
+    def __init__(self, n, field=AA, **kwargs):

         if n == 0:
             jordan_product = lambda x,y: x
             inner_product = lambda x,y: x
         if n == 0:
             jordan_product = lambda x,y: x
             inner_product = lambda x,y: x


@@ -2361,8 +2647,14 @@ class HadamardEJA(ConcreteEJA):
         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
         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)
+        column_basis = tuple( b.column()
+                              for b in FreeModule(field, n).basis() )
+        super().__init__(column_basis,
+                         jordan_product,
+                         inner_product,
+                         field=field,
+                         associative=True,
+                         **kwargs)

         self.rank.set_cache(n)
 
         if n == 0:
         self.rank.set_cache(n)
 
         if n == 0:


@@ -2468,7 +2760,7 @@ class BilinearFormEJA(ConcreteEJA):
         True
 
     """
         True
 
     """

-    def __init__(self, B, **kwargs):
+    def __init__(self, B, field=AA, **kwargs):

         # 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...
         # 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...


@@ -2496,11 +2788,20 @@ class BilinearFormEJA(ConcreteEJA):
             return P([z0] + zbar.list())
 
         n = B.nrows()
             return P([z0] + zbar.list())
 
         n = B.nrows()

-        column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )
-        super(BilinearFormEJA, self).__init__(column_basis,
-                                              jordan_product,
-                                              inner_product,
-                                              **kwargs)
+        column_basis = tuple( b.column()
+                              for b in FreeModule(field, n).basis() )
+
+        # TODO: I haven't actually checked this, but it seems legit.
+        associative = False
+        if n <= 2:
+            associative = True
+
+        super().__init__(column_basis,
+                         jordan_product,
+                         inner_product,
+                         field=field,
+                         associative=associative,
+                         **kwargs)

 
         # The rank of this algebra is two, unless we're in a
         # one-dimensional ambient space (because the rank is bounded
 
         # The rank of this algebra is two, unless we're in a
         # one-dimensional ambient space (because the rank is bounded


@@ -2560,20 +2861,20 @@ class JordanSpinEJA(BilinearFormEJA):
     This multiplication table can be verified by hand::
 
         sage: J = JordanSpinEJA(4)
     This multiplication table can be verified by hand::
 
         sage: J = JordanSpinEJA(4)

-        sage: e0,e1,e2,e3 = J.gens()
-        sage: e0*e0
-        e0
-        sage: e0*e1
-        e1
-        sage: e0*e2
-        e2
-        sage: e0*e3
-        e3
-        sage: e1*e2
+        sage: b0,b1,b2,b3 = J.gens()
+        sage: b0*b0
+        b0
+        sage: b0*b1
+        b1
+        sage: b0*b2
+        b2
+        sage: b0*b3
+        b3
+        sage: b1*b2

         0
         0

-        sage: e1*e3
+        sage: b1*b3

         0
         0

-        sage: e2*e3
+        sage: b2*b3

         0
 
     We can change the generator prefix::
         0
 
     We can change the generator prefix::


@@ -2594,7 +2895,7 @@ class JordanSpinEJA(BilinearFormEJA):
             True
 
     """
             True
 
     """

-    def __init__(self, n, **kwargs):
+    def __init__(self, n, *args, **kwargs):

         # This is a special case of the BilinearFormEJA with the
         # identity matrix as its bilinear form.
         B = matrix.identity(ZZ, n)
         # This is a special case of the BilinearFormEJA with the
         # identity matrix as its bilinear form.
         B = matrix.identity(ZZ, n)


@@ -2605,7 +2906,7 @@ class JordanSpinEJA(BilinearFormEJA):
 
         # But also don't pass check_field=False here, because the user
         # can pass in a field!
 
         # But also don't pass check_field=False here, because the user
         # can pass in a field!

-        super(JordanSpinEJA, self).__init__(B, **kwargs)
+        super().__init__(B, *args, **kwargs)

 
     @staticmethod
     def _max_random_instance_size():
 
     @staticmethod
     def _max_random_instance_size():


@@ -2663,10 +2964,12 @@ class TrivialEJA(ConcreteEJA):
         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 

-        super(TrivialEJA, self).__init__(basis,
-                                         jordan_product,
-                                         inner_product,
-                                         **kwargs)
+        super().__init__(basis,
+                         jordan_product,
+                         inner_product,
+                         associative=True,
+                         **kwargs)
+

         # The rank is zero using my definition, namely the dimension of the
         # largest subalgebra generated by any element.
         self.rank.set_cache(0)
         # The rank is zero using my definition, namely the dimension of the
         # largest subalgebra generated by any element.
         self.rank.set_cache(0)


@@ -2679,8 +2982,7 @@ class TrivialEJA(ConcreteEJA):
         return cls(**kwargs)
 
 
         return cls(**kwargs)
 
 

-class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
-                          FiniteDimensionalEJA):
+class CartesianProductEJA(FiniteDimensionalEJA):

     r"""
     The external (orthogonal) direct sum of two or more Euclidean
     Jordan algebras. Every Euclidean Jordan algebra decomposes into an
     r"""
     The external (orthogonal) direct sum of two or more Euclidean
     Jordan algebras. Every Euclidean Jordan algebra decomposes into an


@@ -2690,7 +2992,8 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
 
     SETUP::
 
 
     SETUP::
 

-        sage: from mjo.eja.eja_algebra import (CartesianProductEJA,
+        sage: from mjo.eja.eja_algebra import (random_eja,
+        ....:                                  CartesianProductEJA,

         ....:                                  HadamardEJA,
         ....:                                  JordanSpinEJA,
         ....:                                  RealSymmetricEJA)
         ....:                                  HadamardEJA,
         ....:                                  JordanSpinEJA,
         ....:                                  RealSymmetricEJA)


@@ -2729,6 +3032,79 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
         Real Field (+) Euclidean Jordan algebra of dimension 6 over
         Algebraic Real Field
 
         Real Field (+) Euclidean Jordan algebra of dimension 6 over
         Algebraic Real Field
 

+    Rank is additive on a Cartesian product::
+
+        sage: J1 = HadamardEJA(1)
+        sage: J2 = RealSymmetricEJA(2)
+        sage: J = cartesian_product([J1,J2])
+        sage: J1.rank.clear_cache()
+        sage: J2.rank.clear_cache()
+        sage: J.rank.clear_cache()
+        sage: J.rank()
+        3
+        sage: J.rank() == J1.rank() + J2.rank()
+        True
+
+    The same rank computation works over the rationals, with whatever
+    basis you like::
+
+        sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
+        sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
+        sage: J = cartesian_product([J1,J2])
+        sage: J1.rank.clear_cache()
+        sage: J2.rank.clear_cache()
+        sage: J.rank.clear_cache()
+        sage: J.rank()
+        3
+        sage: J.rank() == J1.rank() + J2.rank()
+        True
+
+    The product algebra will be associative if and only if all of its
+    components are associative::
+
+        sage: J1 = HadamardEJA(2)
+        sage: J1.is_associative()
+        True
+        sage: J2 = HadamardEJA(3)
+        sage: J2.is_associative()
+        True
+        sage: J3 = RealSymmetricEJA(3)
+        sage: J3.is_associative()
+        False
+        sage: CP1 = cartesian_product([J1,J2])
+        sage: CP1.is_associative()
+        True
+        sage: CP2 = cartesian_product([J1,J3])
+        sage: CP2.is_associative()
+        False
+
+    Cartesian products of Cartesian products work::
+
+        sage: J1 = JordanSpinEJA(1)
+        sage: J2 = JordanSpinEJA(1)
+        sage: J3 = JordanSpinEJA(1)
+        sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
+        sage: J.multiplication_table()
+        +----++----+----+----+
+        | *  || b0 | b1 | b2 |
+        +====++====+====+====+
+        | b0 || b0 | 0  | 0  |
+        +----++----+----+----+
+        | b1 || 0  | b1 | 0  |
+        +----++----+----+----+
+        | b2 || 0  | 0  | b2 |
+        +----++----+----+----+
+        sage: HadamardEJA(3).multiplication_table()
+        +----++----+----+----+
+        | *  || b0 | b1 | b2 |
+        +====++====+====+====+
+        | b0 || b0 | 0  | 0  |
+        +----++----+----+----+
+        | b1 || 0  | b1 | 0  |
+        +----++----+----+----+
+        | b2 || 0  | 0  | b2 |
+        +----++----+----+----+
+

     TESTS:
 
     All factors must share the same base field::
     TESTS:
 
     All factors must share the same base field::


@@ -2740,32 +3116,62 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
         ...
         ValueError: all factors 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, modules, **kwargs):
-        CombinatorialFreeModule_CartesianProduct.__init__(self,
-                                                          modules,
-                                                          **kwargs)
-        field = modules[0].base_ring()
-        if not all( J.base_ring() == field for J in modules ):
+    Element = FiniteDimensionalEJAElement
+
+
+    def __init__(self, factors, **kwargs):
+        m = len(factors)
+        if m == 0:
+            return TrivialEJA()
+
+        self._sets = factors
+
+        field = factors[0].base_ring()
+        if not all( J.base_ring() == field for J in factors ):

             raise ValueError("all factors must share the same base field")
 
             raise ValueError("all factors must share the same base field")
 

-        basis = tuple( b.to_vector().column() for b in self.basis() )
+        associative = all( f.is_associative() for f in factors )

 
 

-        # 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()
+        MS = self.matrix_space()
+        basis = []
+        zero = MS.zero()
+        for i in range(m):
+            for b in factors[i].matrix_basis():
+                z = list(zero)
+                z[i] = b
+                basis.append(z)

 
 

-        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)
+        basis = tuple( MS(b) for b in basis )

 
 

-        # 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.
+        # Define jordan/inner products that operate on that matrix_basis.
+        def jordan_product(x,y):
+            return MS(tuple(
+                (factors[i](x[i])*factors[i](y[i])).to_matrix()
+                for i in range(m)
+            ))
+
+        def inner_product(x, y):
+            return sum(
+                factors[i](x[i]).inner_product(factors[i](y[i]))
+                for i in range(m)
+            )
+
+        # There's no need to check the field since it already came
+        # from an EJA. Likewise the axioms are guaranteed to be
+        # satisfied, unless the guy writing this class sucks.

         #
         # If you want the basis to be orthonormalized, orthonormalize
         # the factors.
         #
         # If you want the basis to be orthonormalized, orthonormalize
         # the factors.


@@ -2775,11 +3181,55 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
                                       inner_product,
                                       field=field,
                                       orthonormalize=False,
                                       inner_product,
                                       field=field,
                                       orthonormalize=False,

+                                      associative=associative,
+                                      cartesian_product=True,

                                       check_field=False,
                                       check_field=False,

-                                      check_axioms=False,
-                                      category=self.category())
+                                      check_axioms=False)
+
+        ones = tuple(J.one().to_matrix() for J in factors)
+        self.one.set_cache(self(ones))
+        self.rank.set_cache(sum(J.rank() for J in factors))
+
+    def cartesian_factors(self):
+        # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
+        return self._sets
+
+    def cartesian_factor(self, i):
+        r"""
+        Return the ``i``th factor of this algebra.
+        """
+        return self._sets[i]
+
+    def _repr_(self):
+        # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
+        from sage.categories.cartesian_product import cartesian_product
+        return cartesian_product.symbol.join("%s" % factor
+                                             for factor in self._sets)

 
 

-        self.rank.set_cache(sum(J.rank() for J in modules))
+    def matrix_space(self):
+        r"""
+        Return the space that our matrix basis lives in as a Cartesian
+        product.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            ....:                                  RealSymmetricEJA)
+
+        EXAMPLES::
+
+            sage: J1 = HadamardEJA(1)
+            sage: J2 = RealSymmetricEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: J.matrix_space()
+            The Cartesian product of (Full MatrixSpace of 1 by 1 dense
+            matrices over Algebraic Real Field, Full MatrixSpace of 2
+            by 2 dense matrices over Algebraic Real Field)
+
+        """
+        from sage.categories.cartesian_product import cartesian_product
+        return cartesian_product( [J.matrix_space()
+                                   for J in self.cartesian_factors()] )

 
     @cached_method
     def cartesian_projection(self, i):
 
     @cached_method
     def cartesian_projection(self, i):


@@ -2851,9 +3301,12 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
             True
 
         """
             True
 
         """

-        Ji = self.cartesian_factors()[i]
-        # Requires the fix on Trac 31421/31422 to work!
-        Pi = super().cartesian_projection(i)
+        offset = sum( self.cartesian_factor(k).dimension()
+                      for k in range(i) )
+        Ji = self.cartesian_factor(i)
+        Pi = self._module_morphism(lambda j: Ji.monomial(j - offset),
+                                   codomain=Ji)
+

         return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix())
 
     @cached_method
         return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix())
 
     @cached_method


@@ -2959,90 +3412,64 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
             True
 
         """
             True
 
         """

-        Ji = self.cartesian_factors()[i]
-        # Requires the fix on Trac 31421/31422 to work!
-        Ei = super().cartesian_embedding(i)
+        offset = sum( self.cartesian_factor(k).dimension()
+                      for k in range(i) )
+        Ji = self.cartesian_factor(i)
+        Ei = Ji._module_morphism(lambda j: self.monomial(j + offset),
+                                 codomain=self)

         return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix())
 
 
         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.
+FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA

 
 

-        SETUP::
+class RationalBasisCartesianProductEJA(CartesianProductEJA,
+                                       RationalBasisEJA):
+    r"""
+    A separate class for products of algebras for which we know a
+    rational basis.

 
 

-            sage: from mjo.eja.eja_algebra import (HadamardEJA,
-            ....:                                  QuaternionHermitianEJA)
+    SETUP::

 
 

-        EXAMPLE::
+        sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+        ....:                                  RealSymmetricEJA)

 
 

-            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
+    EXAMPLES:

 
 

-        """
-        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 )
+    This gives us fast characteristic polynomial computations in
+    product algebras, too::

 
 
 
 

-    Element = FiniteDimensionalEJAElement
+        sage: J1 = JordanSpinEJA(2)
+        sage: J2 = RealSymmetricEJA(3)
+        sage: J = cartesian_product([J1,J2])
+        sage: J.characteristic_polynomial_of().degree()
+        5
+        sage: J.rank()
+        5

 
 

+    """
+    def __init__(self, algebras, **kwargs):
+        CartesianProductEJA.__init__(self, algebras, **kwargs)

 
 

-FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA
-random_eja = ConcreteEJA.random_instance
+        self._rational_algebra = None
+        if self.vector_space().base_field() is not QQ:
+            self._rational_algebra = cartesian_product([
+                r._rational_algebra for r in algebras
+            ])
+
+
+RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
+
+def random_eja(*args, **kwargs):
+    J1 = ConcreteEJA.random_instance(*args, **kwargs)
+
+    # This might make Cartesian products appear roughly as often as
+    # any other ConcreteEJA.
+    if ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1) == 0:
+        # Use random_eja() again so we can get more than two factors.
+        J2 = random_eja(*args, **kwargs)
+        J = cartesian_product([J1,J2])
+        return J
+    else:
+        return J1