]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: factor out a class for real-embedded matrices.
[sage.d.git] / mjo / eja / eja_algebra.py
index a77bb871def9264ac0a82dfea4737b6c2b4411fc..5bf597565b2ae292032c3f0a532763d7889cd2a8 100644 (file)
@@ -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::
 
@@ -13,7 +57,6 @@ EXAMPLES::
 
     sage: random_eja()
     Euclidean Jordan algebra of dimension...
-
 """
 
 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.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
@@ -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
-from mjo.eja.eja_utils import _mat2vec
+from mjo.eja.eja_utils import _all2list, _mat2vec
 
 class FiniteDimensionalEJA(CombinatorialFreeModule):
     r"""
@@ -41,16 +83,50 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
     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
@@ -61,11 +137,13 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                  inner_product,
                  field=AA,
                  orthonormalize=True,
-                 associative=False,
+                 associative=None,
+                 cartesian_product=False,
                  check_field=True,
                  check_axioms=True,
-                 prefix='e',
-                 category=None):
+                 prefix="b"):
+
+        n = len(basis)
 
         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")
 
-        # 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
@@ -94,33 +168,49 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                 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.
-        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
-        # 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:
-            # 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."
@@ -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.
-            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))
@@ -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).
-        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:
@@ -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)
-                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:
@@ -226,6 +316,35 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
 
     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:
@@ -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
+
         """
         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.
@@ -296,9 +437,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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"""
@@ -317,6 +456,92 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                     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
@@ -326,11 +551,14 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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()):
@@ -340,12 +568,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                     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
 
@@ -359,7 +583,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            sage: from mjo.eja.eja_algebra import (random_eja,
+            ....:                                  JordanSpinEJA,
             ....:                                  HadamardEJA,
             ....:                                  RealSymmetricEJA)
 
@@ -381,29 +606,42 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             ...
             ValueError: not an element of this algebra
 
+        Tuples work as well, provided that the matrix basis for the
+        algebra consists of them::
+
+            sage: J1 = HadamardEJA(3)
+            sage: J2 = RealSymmetricEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
+            b1 + b5
+
         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: 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: J(x.to_vector().column()) == x
+            sage: J(x.to_matrix()) == x
             True
 
+        We cannot coerce elements between algebras just because their
+        matrix representations are compatible::
+
+            sage: J1 = HadamardEJA(3)
+            sage: J2 = JordanSpinEJA(3)
+            sage: J2(J1.one())
+            Traceback (most recent call last):
+            ...
+            ValueError: not an element of this algebra
+            sage: J1(J2.zero())
+            Traceback (most recent call last):
+            ...
+            ValueError: not an element of this algebra
         """
         msg = "not an element of this algebra"
-        if elt == 0:
-            # The superclass implementation of random_element()
-            # needs to be able to coerce "0" into the algebra.
-            return self.zero()
-        elif elt in self.base_ring():
+        if elt in self.base_ring():
             # Ensure that no base ring -> algebra coercion is performed
             # by this method. There's some stupidity in sage that would
             # otherwise propagate to this method; for example, sage thinks
@@ -411,9 +649,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             raise ValueError(msg)
 
         try:
+            # Try to convert a vector into a column-matrix...
             elt = elt.column()
         except (AttributeError, TypeError):
-            # Try to convert a vector into a column-matrix
+            # and ignore failure, because we weren't really expecting
+            # a vector as an argument anyway.
             pass
 
         if elt not in self.matrix_space():
@@ -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.
         #
+        # And, we also have to handle Cartesian product bases (when
+        # the matrix basis consists of tuples) here. The "good news"
+        # is that we're already converting everything to long vectors,
+        # and that strategy works for tuples as well.
+        #
         # We pass check=False because the matrix basis is "guaranteed"
         # to be linearly independent... right? Ha ha.
-        V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols())
-        W = V.span_of_basis( (_mat2vec(s) for s in self.matrix_basis()),
+        elt = _all2list(elt)
+        V = VectorSpace(self.base_ring(), len(elt))
+        W = V.span_of_basis( (V(_all2list(s)) for s in self.matrix_basis()),
                              check=False)
 
         try:
-            coords =  W.coordinate_vector(_mat2vec(elt))
+            coords = W.coordinate_vector(V(elt))
         except ArithmeticError:  # vector is not in free module
             raise ValueError(msg)
 
@@ -642,15 +888,15 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             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."
-        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)
@@ -704,7 +950,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             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]
@@ -715,7 +961,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             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]
@@ -731,12 +977,49 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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
-        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)
@@ -760,20 +1043,20 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             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()
-            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()
-            f0
+            c0
             sage: A.one().superalgebra_element()
-            e0 + e1 + e2 + e3 + e4
+            b0 + b1 + b2 + b3 + b4
 
         TESTS:
 
@@ -997,14 +1280,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
@@ -1014,9 +1295,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:
@@ -1235,6 +1514,14 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         return len(self._charpoly_coefficients())
 
 
+    def subalgebra(self, basis, **kwargs):
+        r"""
+        Create a subalgebra of this algebra from the given basis.
+        """
+        from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
+        return FiniteDimensionalEJASubalgebra(self, basis, **kwargs)
+
+
     def vector_space(self):
         """
         Return the vector space that underlies this algebra.
@@ -1253,7 +1540,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         return self.zero().to_vector().parent().ambient_vector_space()
 
 
-    Element = FiniteDimensionalEJAElement
 
 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")
 
+        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
@@ -1304,17 +1597,11 @@ class RationalBasisEJA(FiniteDimensionalEJA):
                                        jordan_product,
                                        inner_product,
                                        field=QQ,
+                                       associative=self.is_associative(),
                                        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"""
@@ -1438,6 +1725,21 @@ class ConcreteEJA(RationalBasisEJA):
 
 
 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"""
@@ -1483,9 +1785,6 @@ class MatrixEJA:
             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):
@@ -1494,29 +1793,11 @@ class MatrixEJA:
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
-            ....:                                  ComplexHermitianEJA,
+            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
             ....:                                  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)
@@ -1544,27 +1825,15 @@ class MatrixEJA:
             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
@@ -1577,19 +1846,19 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
     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::
 
-        sage: RealSymmetricEJA(2, field=RDF)
+        sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
         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
 
@@ -1630,7 +1899,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
 
     """
     @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.
 
@@ -1642,7 +1911,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
 
             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
 
@@ -1652,7 +1921,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
         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:
@@ -1673,26 +1942,32 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
         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
 
-        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)
-        idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
+        idV = self.matrix_space().one()
         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 = {}
@@ -1771,7 +2046,7 @@ class ComplexMatrixEJA(MatrixEJA):
             True
 
         """
-        super(ComplexMatrixEJA,cls).real_embed(M)
+        super().real_embed(M)
         n = M.nrows()
 
         # We don't need any adjoined elements...
@@ -1818,7 +2093,7 @@ class ComplexMatrixEJA(MatrixEJA):
             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())
@@ -1855,9 +2130,9 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
 
     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
-        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
 
@@ -1899,7 +2174,7 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
     """
 
     @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.
 
@@ -1917,15 +2192,14 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
 
             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
 
         """
-        field = ZZ
-        R = PolynomialRing(field, 'z')
+        R = PolynomialRing(ZZ, 'z')
         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:
@@ -1956,20 +2230,26 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
                 # "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 )
 
 
-    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
 
-        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().
@@ -1989,7 +2269,7 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
         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.
@@ -2052,7 +2332,7 @@ class QuaternionMatrixEJA(MatrixEJA):
             True
 
         """
-        super(QuaternionMatrixEJA,cls).real_embed(M)
+        super().real_embed(M)
         quaternions = M.base_ring()
         n = M.nrows()
 
@@ -2107,7 +2387,7 @@ class QuaternionMatrixEJA(MatrixEJA):
             True
 
         """
-        super(QuaternionMatrixEJA,cls).real_unembed(M)
+        super().real_unembed(M)
         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::
 
-        sage: QuaternionHermitianEJA(2, field=RDF)
+        sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
         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
 
@@ -2195,7 +2475,7 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
 
     """
     @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.
 
@@ -2213,12 +2493,11 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
 
             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
 
         """
-        field = ZZ
         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
 
-        # 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 )
 
 
-    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
 
-        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().
@@ -2318,19 +2604,19 @@ class HadamardEJA(ConcreteEJA):
     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
-        sage: e0*e2
+        sage: b0*b2
         0
-        sage: e1*e1
-        e1
-        sage: e1*e2
+        sage: b1*b1
+        b1
+        sage: b1*b2
         0
-        sage: e2*e2
-        e2
+        sage: b2*b2
+        b2
 
     TESTS:
 
@@ -2340,7 +2626,7 @@ class HadamardEJA(ConcreteEJA):
         (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
@@ -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
 
-        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:
@@ -2468,7 +2760,7 @@ class BilinearFormEJA(ConcreteEJA):
         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...
@@ -2496,11 +2788,20 @@ class BilinearFormEJA(ConcreteEJA):
             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
@@ -2560,20 +2861,20 @@ class JordanSpinEJA(BilinearFormEJA):
     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
-        sage: e1*e3
+        sage: b1*b3
         0
-        sage: e2*e3
+        sage: b2*b3
         0
 
     We can change the generator prefix::
@@ -2594,7 +2895,7 @@ class JordanSpinEJA(BilinearFormEJA):
             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)
@@ -2605,7 +2906,7 @@ class JordanSpinEJA(BilinearFormEJA):
 
         # 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():
@@ -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
 
-        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)
@@ -2679,8 +2982,7 @@ class TrivialEJA(ConcreteEJA):
         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
@@ -2690,7 +2992,8 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
 
     SETUP::
 
-        sage: from mjo.eja.eja_algebra import (CartesianProductEJA,
+        sage: from mjo.eja.eja_algebra import (random_eja,
+        ....:                                  CartesianProductEJA,
         ....:                                  HadamardEJA,
         ....:                                  JordanSpinEJA,
         ....:                                  RealSymmetricEJA)
@@ -2729,6 +3032,79 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
         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::
@@ -2740,32 +3116,62 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
         ...
         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")
 
-        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.
@@ -2775,11 +3181,55 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
                                       inner_product,
                                       field=field,
                                       orthonormalize=False,
+                                      associative=associative,
+                                      cartesian_product=True,
                                       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):
@@ -2851,9 +3301,12 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
             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
@@ -2959,90 +3412,64 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
             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())
 
 
-    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