]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
COPYING,LICENSE: add (AGPL-3.0+)
[sage.d.git] / mjo / eja / eja_algebra.py
index 631d695aa75208a65cb457a4ca11280f44823765..adcc3436b1302e09cd20007d0525aee08e32a48f 100644 (file)
@@ -1,4 +1,4 @@
-"""
+r"""
 Representations and constructions for Euclidean Jordan algebras.
 
 A Euclidean Jordan algebra is a Jordan algebra that has some
 Representations and constructions for Euclidean Jordan algebras.
 
 A Euclidean Jordan algebra is a Jordan algebra that has some
@@ -32,22 +32,118 @@ for these simple algebras:
   * :class:`RealSymmetricEJA`
   * :class:`ComplexHermitianEJA`
   * :class:`QuaternionHermitianEJA`
   * :class:`RealSymmetricEJA`
   * :class:`ComplexHermitianEJA`
   * :class:`QuaternionHermitianEJA`
+  * :class:`OctonionHermitianEJA`
 
 
-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,
+In addition to these, we provide a few other example constructions,
 
 
+  * :class:`JordanSpinEJA`
   * :class:`HadamardEJA`
   * :class:`HadamardEJA`
+  * :class:`AlbertEJA`
   * :class:`TrivialEJA`
   * :class:`TrivialEJA`
+  * :class:`ComplexSkewSymmetricEJA`
 
 The Jordan spin algebra is a bilinear form algebra where the bilinear
 form is the identity. The Hadamard EJA is simply a Cartesian product
 
 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.
+of one-dimensional spin algebras. The Albert EJA is simply a special
+case of the :class:`OctonionHermitianEJA` where the matrices are
+three-by-three and the resulting space has dimension 27. And
+last/least, the trivial EJA is exactly what you think it is; it could
+also be obtained by constructing a dimension-zero instance of any of
+the other algebras. 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.
+
+At a minimum, the following are required to construct a Euclidean
+Jordan algebra:
+
+  * A basis of matrices, column vectors, or MatrixAlgebra elements
+  * A Jordan product defined on the basis
+  * Its inner product defined on the basis
+
+The real numbers form a Euclidean Jordan algebra when both the Jordan
+and inner products are the usual multiplication. We use this as our
+example, and demonstrate a few ways to construct an EJA.
+
+First, we can use one-by-one SageMath matrices with algebraic real
+entries to represent real numbers. We define the Jordan and inner
+products to be essentially real-number multiplication, with the only
+difference being that the Jordan product again returns a one-by-one
+matrix, whereas the inner product must return a scalar. Our basis for
+the one-by-one matrices is of course the set consisting of a single
+matrix with its sole entry non-zero::
+
+    sage: from mjo.eja.eja_algebra import EJA
+    sage: jp = lambda X,Y: X*Y
+    sage: ip = lambda X,Y: X[0,0]*Y[0,0]
+    sage: b1 = matrix(AA, [[1]])
+    sage: J1 = EJA((b1,), jp, ip)
+    sage: J1
+    Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
+
+In fact, any positive scalar multiple of that inner-product would work::
+
+    sage: ip2 = lambda X,Y: 16*ip(X,Y)
+    sage: J2 = EJA((b1,), jp, ip2)
+    sage: J2
+    Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
+
+But beware that your basis will be orthonormalized _with respect to the
+given inner-product_ unless you pass ``orthonormalize=False`` to the
+constructor. For example::
+
+    sage: J3 = EJA((b1,), jp, ip2, orthonormalize=False)
+    sage: J3
+    Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
+
+To see the difference, you can take the first and only basis element
+of the resulting algebra, and ask for it to be converted back into
+matrix form::
+
+    sage: J1.basis()[0].to_matrix()
+    [1]
+    sage: J2.basis()[0].to_matrix()
+    [1/4]
+    sage: J3.basis()[0].to_matrix()
+    [1]
+
+Since square roots are used in that process, the default scalar field
+that we use is the field of algebraic real numbers, ``AA``. You can
+also Use rational numbers, but only if you either pass
+``orthonormalize=False`` or know that orthonormalizing your basis
+won't stray beyond the rational numbers. The example above would
+have worked only because ``sqrt(16) == 4`` is rational.
+
+Another option for your basis is to use elemebts of a
+:class:`MatrixAlgebra`::
+
+    sage: from mjo.matrix_algebra import MatrixAlgebra
+    sage: A = MatrixAlgebra(1,AA,AA)
+    sage: J4 = EJA(A.gens(), jp, ip)
+    sage: J4
+    Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
+    sage: J4.basis()[0].to_matrix()
+    +---+
+    | 1 |
+    +---+
+
+An easier way to view the entire EJA basis in its original (but
+perhaps orthonormalized) matrix form is to use the ``matrix_basis``
+method::
+
+    sage: J4.matrix_basis()
+    (+---+
+    | 1 |
+     +---+,)
+
+In particular, a :class:`MatrixAlgebra` is needed to work around the
+fact that matrices in SageMath must have entries in the same
+(commutative and associative) ring as its scalars. There are many
+Euclidean Jordan algebras whose elements are matrices that violate
+those assumptions. The complex, quaternion, and octonion Hermitian
+matrices all have entries in a ring (the complex numbers, quaternions,
+or octonions...) that differs from the algebra's scalar ring (the real
+numbers). Quaternions are also non-commutative; the octonions are
+neither commutative nor associative.
 
 SETUP::
 
 
 SETUP::
 
@@ -59,13 +155,10 @@ EXAMPLES::
     Euclidean Jordan algebra of dimension...
 """
 
     Euclidean Jordan algebra of dimension...
 """
 
-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
@@ -74,107 +167,145 @@ from sage.modules.free_module import FreeModule, VectorSpace
 from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
                             PolynomialRing,
                             QuadraticField)
 from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
                             PolynomialRing,
                             QuadraticField)
-from mjo.eja.eja_element import FiniteDimensionalEJAElement
-from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
-from mjo.eja.eja_utils import _all2list, _mat2vec
+from mjo.eja.eja_element import (CartesianProductEJAElement,
+                                 EJAElement)
+from mjo.eja.eja_operator import EJAOperator
+from mjo.eja.eja_utils import _all2list
+
+def EuclideanJordanAlgebras(field):
+    r"""
+    The category of Euclidean Jordan algebras over ``field``, which
+    must be a subfield of the real numbers. For now this is just a
+    convenient wrapper around all of the other category axioms that
+    apply to all EJAs.
+    """
+    category = MagmaticAlgebras(field).FiniteDimensional()
+    category = category.WithBasis().Unital().Commutative()
+    return category
 
 
-class FiniteDimensionalEJA(CombinatorialFreeModule):
+class EJA(CombinatorialFreeModule):
     r"""
     A finite-dimensional Euclidean Jordan algebra.
 
     INPUT:
 
     r"""
     A finite-dimensional Euclidean Jordan algebra.
 
     INPUT:
 
-      - basis -- a tuple of basis elements in "matrix form," which
-        must be the same form as the arguments to ``jordan_product``
-        and ``inner_product``. In reality, "matrix form" can be either
-        vectors, matrices, or a Cartesian product (ordered tuple)
-        of vectors or matrices. All of these would ideally be vector
-        spaces in sage with no special-casing needed; but in reality
-        we turn vectors into column-matrices and Cartesian products
-        `(a,b)` into column matrices `(a,b)^{T}` after converting
-        `a` and `b` themselves.
-
-      - jordan_product -- function of two ``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 -- 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.
+      - ``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.
+
+      - ``matrix_space`` -- the space that your matrix basis lives in,
+        or ``None`` (the default). So long as your basis does not have
+        length zero you can omit this. But in trivial algebras, it is
+        required.
+
+      - ``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
+
+    TESTS:
+
+    We should compute that an element subalgebra is associative even
+    if we circumvent the element method::
+
+        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 = EJAElement
+
+    @staticmethod
+    def _check_input_field(field):
+        if not field.is_subring(RR):
+            # Note: this does return true for the real algebraic
+            # field, the rationals, and any quadratic field where
+            # we've specified a real embedding.
+            raise ValueError("scalar field is not real")
+
+    @staticmethod
+    def _check_input_axioms(basis, jordan_product, inner_product):
+        if not all( jordan_product(bi,bj) == jordan_product(bj,bi)
+                    for bi in basis
+                    for bj in basis ):
+            raise ValueError("Jordan product is not commutative")
+
+        if not all( inner_product(bi,bj) == inner_product(bj,bi)
+                    for bi in basis
+                    for bj in basis ):
+            raise ValueError("inner-product is not commutative")
 
     def __init__(self,
                  basis,
                  jordan_product,
                  inner_product,
                  field=AA,
 
     def __init__(self,
                  basis,
                  jordan_product,
                  inner_product,
                  field=AA,
+                 matrix_space=None,
                  orthonormalize=True,
                  orthonormalize=True,
-                 associative=False,
-                 cartesian_product=False,
+                 associative=None,
                  check_field=True,
                  check_axioms=True,
                  check_field=True,
                  check_axioms=True,
-                 prefix='e'):
-
-        # Keep track of whether or not the matrix basis consists of
-        # tuples, since we need special cases for them damned near
-        # everywhere.  This is INDEPENDENT of whether or not the
-        # algebra is a cartesian product, since a subalgebra of a
-        # cartesian product will have a basis of tuples, but will not
-        # in general itself be a cartesian product algebra.
-        self._matrix_basis_is_cartesian = False
+                 prefix="b"):
+
         n = len(basis)
         n = len(basis)
-        if n > 0:
-            if hasattr(basis[0], 'cartesian_factors'):
-                self._matrix_basis_is_cartesian = True
 
         if check_field:
 
         if check_field:
-            if not field.is_subring(RR):
-                # Note: this does return true for the real algebraic
-                # field, the rationals, and any quadratic field where
-                # 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.
-        if not cartesian_product:
-            # The field for a cartesian product algebra comes from one
-            # of its factors and is the same for all factors, so
-            # there's no need to "reapply" it on product algebras.
-            if self._matrix_basis_is_cartesian:
-                # OK since if n == 0, the basis does not consist of tuples.
-                P = basis[0].parent()
-                basis = tuple( P(tuple(b_i.change_ring(field) for b_i in b))
-                               for b in basis )
-            else:
-                basis = tuple( b.change_ring(field) for b in basis )
-
+            self._check_input_field(field)
 
         if check_axioms:
             # Check commutativity of the Jordan and inner-products.
             # This has to be done before we build the multiplication
             # and inner-product tables/matrices, because we take
             # advantage of symmetry in the process.
 
         if check_axioms:
             # Check commutativity of the Jordan and inner-products.
             # This has to be done before we build the multiplication
             # and inner-product tables/matrices, because we take
             # advantage of symmetry in the process.
-            if not all( jordan_product(bi,bj) == jordan_product(bj,bi)
-                        for bi in basis
-                        for bj in basis ):
-                raise ValueError("Jordan product is not commutative")
-
-            if not all( inner_product(bi,bj) == inner_product(bj,bi)
-                        for bi in basis
-                        for bj in basis ):
-                raise ValueError("inner-product is not commutative")
+            self._check_input_axioms(basis, jordan_product, inner_product)
+
+        if n <= 1:
+            # All zero- and one-dimensional algebras are just the real
+            # numbers with (some positive multiples of) the usual
+            # multiplication as its Jordan and inner-product.
+            associative = True
+        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)
+
+        category = EuclideanJordanAlgebras(field)
 
 
-
-        category = MagmaticAlgebras(field).FiniteDimensional()
-        category = category.WithBasis().Unital()
         if associative:
             # Element subalgebras can take advantage of this.
             category = category.Associative()
         if associative:
             # Element subalgebras can take advantage of this.
             category = category.Associative()
-        if cartesian_product:
-            category = category.CartesianProducts()
 
         # Call the superclass constructor so that we can use its from_vector()
         # method to build our multiplication table.
 
         # Call the superclass constructor so that we can use its from_vector()
         # method to build our multiplication table.
@@ -189,8 +320,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # 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
         # 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.
-        vector_basis = basis
+        # we see in things like x = 1*b1 + 2*b2.
 
         degree = 0
         if n > 0:
 
         degree = 0
         if n > 0:
@@ -200,9 +330,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # written out as "long vectors."
         V = VectorSpace(field, degree)
 
         # written out as "long vectors."
         V = VectorSpace(field, degree)
 
-        # The matrix that will hole the orthonormal -> unorthonormal
-        # coordinate transformation.
-        self._deortho_matrix = None
+        # The matrix that will hold the orthonormal -> unorthonormal
+        # coordinate transformation. Default to an identity matrix of
+        # the appropriate size to avoid special cases for None
+        # everywhere.
+        self._deortho_matrix = matrix.identity(field,n)
 
         if orthonormalize:
             # Save a copy of the un-orthonormalized basis for later.
 
         if orthonormalize:
             # Save a copy of the un-orthonormalized basis for later.
@@ -214,30 +346,42 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             basis = tuple(gram_schmidt(basis, inner_product))
 
         # Save the (possibly orthonormalized) matrix basis for
             basis = tuple(gram_schmidt(basis, inner_product))
 
         # Save the (possibly orthonormalized) matrix basis for
-        # later...
+        # later, as well as the space that its elements live in.
+        # In most cases we can deduce the matrix space, but when
+        # n == 0 (that is, there are no basis elements) we cannot.
         self._matrix_basis = basis
         self._matrix_basis = basis
+        if matrix_space is None:
+            self._matrix_space = self._matrix_basis[0].parent()
+        else:
+            self._matrix_space = matrix_space
 
         # 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(_all2list(b)) for b in 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(_all2list(b)) for b in basis )
-        W = V.span_of_basis( vector_basis, check=check_axioms)
+
+        # Save the span of our matrix basis (when written out as long
+        # vectors) because otherwise we'll have to reconstruct it
+        # every time we want to coerce a matrix into the algebra.
+        self._matrix_span = V.span_of_basis( vector_basis, check=check_axioms)
 
         if orthonormalize:
 
         if orthonormalize:
-            # Now "W" is the vector space of our algebra coordinates. The
-            # variables "X1", "X2",...  refer to the entries of vectors in
-            # W. Thus to convert back and forth between the orthonormal
-            # coordinates and the given ones, we need to stick the original
-            # basis in W.
+            # Now "self._matrix_span" is the vector space of our
+            # algebra coordinates. The variables "X0", "X1",...  refer
+            # to the entries of vectors in self._matrix_span. Thus to
+            # convert back and forth between the orthonormal
+            # coordinates and the given ones, we need to stick the
+            # original basis in self._matrix_span.
             U = V.span_of_basis( deortho_vector_basis, check=check_axioms)
             U = V.span_of_basis( deortho_vector_basis, check=check_axioms)
-            self._deortho_matrix = matrix( U.coordinate_vector(q)
-                                           for q in vector_basis )
+            self._deortho_matrix = matrix.column( U.coordinate_vector(q)
+                                                  for q in vector_basis )
 
 
         # Now we actually compute the multiplication and inner-product
         # tables/matrices using the possibly-orthonormalized basis.
         self._inner_product_matrix = matrix.identity(field, n)
 
 
         # Now we actually compute the multiplication and inner-product
         # tables/matrices using the possibly-orthonormalized basis.
         self._inner_product_matrix = matrix.identity(field, n)
-        self._multiplication_table = [ [0 for j in range(i+1)]
+        zed = self.zero()
+        self._multiplication_table = [ [zed for j in range(i+1)]
                                        for i in range(n) ]
 
         # Note: the Jordan and inner-products are defined in terms
                                        for i in range(n) ]
 
         # Note: the Jordan and inner-products are defined in terms
@@ -252,7 +396,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(_all2list(elt)))
+                elt = self._matrix_span.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:
@@ -288,7 +432,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         TESTS::
 
 
         TESTS::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: J(1)
             Traceback (most recent call last):
             sage: J = random_eja()
             sage: J(1)
             Traceback (most recent call last):
@@ -313,19 +456,18 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         TESTS::
 
 
         TESTS::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: n = J.dimension()
             sage: J = random_eja()
             sage: n = J.dimension()
-            sage: ei = J.zero()
-            sage: ej = J.zero()
-            sage: ei_ej = J.zero()*J.zero()
+            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)
             sage: if n > 0:
             ....:     i = ZZ.random_element(n)
             ....:     j = ZZ.random_element(n)
-            ....:     ei = J.gens()[i]
-            ....:     ej = J.gens()[j]
-            ....:     ei_ej = J.product_on_basis(i,j)
-            sage: ei*ej == ei_ej
+            ....:     bi = J.monomial(i)
+            ....:     bj = J.monomial(j)
+            ....:     bi_bj = J.product_on_basis(i,j)
+            sage: bi*bj == bi_bj
             True
 
         """
             True
 
         """
@@ -355,7 +497,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         Our inner product is "associative," which means the following for
         a symmetric bilinear form::
 
         Our inner product is "associative," which means the following for
         a symmetric bilinear form::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: x,y,z = J.random_elements(3)
             sage: (x*y).inner_product(z) == y.inner_product(x*z)
             sage: J = random_eja()
             sage: x,y,z = J.random_elements(3)
             sage: (x*y).inner_product(z) == y.inner_product(x*z)
@@ -366,7 +507,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         Ensure that this is the usual inner product for the algebras
         over `R^n`::
 
         Ensure that this is the usual inner product for the algebras
         over `R^n`::
 
-            sage: set_random_seed()
             sage: J = HadamardEJA.random_instance()
             sage: x,y = J.random_elements(2)
             sage: actual = x.inner_product(y)
             sage: J = HadamardEJA.random_instance()
             sage: x,y = J.random_elements(2)
             sage: actual = x.inner_product(y)
@@ -379,7 +519,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         one). This is in Faraut and Koranyi, and also my "On the
         symmetry..." paper::
 
         one). This is in Faraut and Koranyi, and also my "On the
         symmetry..." paper::
 
-            sage: set_random_seed()
             sage: J = BilinearFormEJA.random_instance()
             sage: n = J.dimension()
             sage: x = J.random_element()
             sage: J = BilinearFormEJA.random_instance()
             sage: n = J.dimension()
             sage: x = J.random_element()
@@ -413,6 +552,16 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         """
         return "Associative" in self.category().axioms()
 
         """
         return "Associative" in self.category().axioms()
 
+    def _is_commutative(self):
+        r"""
+        Whether or not this algebra's multiplication table is commutative.
+
+        This method should of course always return ``True``, unless
+        this algebra was constructed with ``check_axioms=False`` and
+        passed an invalid multiplication table.
+        """
+        return all( x*y == y*x for x in self.gens() for y in self.gens() )
+
     def _is_jordanian(self):
         r"""
         Whether or not this algebra's multiplication table respects the
     def _is_jordanian(self):
         r"""
         Whether or not this algebra's multiplication table respects the
@@ -420,16 +569,101 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         We only check one arrangement of `x` and `y`, so for a
         ``True`` result to be truly true, you should also check
 
         We only check one arrangement of `x` and `y`, so for a
         ``True`` result to be truly true, you should also check
-        :meth:`is_commutative`. This method should of course always
+        :meth:`_is_commutative`. This method should of course always
         return ``True``, unless this algebra was constructed with
         ``check_axioms=False`` and passed an invalid multiplication table.
         """
         return ``True``, unless this algebra was constructed with
         ``check_axioms=False`` and passed an invalid multiplication table.
         """
-        return all( (self.gens()[i]**2)*(self.gens()[i]*self.gens()[j])
+        return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
                     ==
                     ==
-                    (self.gens()[i])*((self.gens()[i]**2)*self.gens()[j])
+                    (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j))
                     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: 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
@@ -439,40 +673,40 @@ 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 k in range(self.dimension()):
 
         for i in range(self.dimension()):
             for j in range(self.dimension()):
                 for k in range(self.dimension()):
-                    x = self.gens()[i]
-                    y = self.gens()[j]
-                    z = self.gens()[k]
+                    x = self.monomial(i)
+                    y = self.monomial(j)
+                    z = self.monomial(k)
                     diff = (x*y).inner_product(z) - x.inner_product(y*z)
 
                     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
 
     def _element_constructor_(self, elt):
         """
 
         return True
 
     def _element_constructor_(self, elt):
         """
-        Construct an element of this algebra from its vector or matrix
-        representation.
+        Construct an element of this algebra or a subalgebra from its
+        EJA element, vector, or matrix representation.
 
         This gets called only after the parent element _call_ method
         fails to find a coercion for the argument.
 
         SETUP::
 
 
         This gets called only after the parent element _call_ method
         fails to find a coercion for the argument.
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            sage: from mjo.eja.eja_algebra import (random_eja,
+            ....:                                  JordanSpinEJA,
             ....:                                  HadamardEJA,
             ....:                                  RealSymmetricEJA)
 
             ....:                                  HadamardEJA,
             ....:                                  RealSymmetricEJA)
 
@@ -501,22 +735,26 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: J2 = RealSymmetricEJA(2)
             sage: J = cartesian_product([J1,J2])
             sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
             sage: J2 = RealSymmetricEJA(2)
             sage: J = cartesian_product([J1,J2])
             sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
-            e(0, 1) + e(1, 2)
+            b1 + b5
+
+        Subalgebra elements are embedded into the superalgebra::
+
+            sage: J = JordanSpinEJA(3)
+            sage: J.one()
+            b0
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by()
+            sage: J(A.one())
+            b0
 
         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: 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
 
         We cannot coerce elements between algebras just because their
             True
 
         We cannot coerce elements between algebras just because their
@@ -542,13 +780,16 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             # that the integer 3 belongs to the space of 2-by-2 matrices.
             raise ValueError(msg)
 
             # that the integer 3 belongs to the space of 2-by-2 matrices.
             raise ValueError(msg)
 
-        try:
-            # Try to convert a vector into a column-matrix...
+        if hasattr(elt, 'superalgebra_element'):
+            # Handle subalgebra elements
+            if elt.parent().superalgebra() == self:
+                return elt.superalgebra_element()
+
+        if hasattr(elt, 'sparse_vector'):
+            # Convert a vector into a column-matrix. We check for
+            # "sparse_vector" and not "column" because matrices also
+            # have a "column" method.
             elt = elt.column()
             elt = elt.column()
-        except (AttributeError, TypeError):
-            # and ignore failure, because we weren't really expecting
-            # a vector as an argument anyway.
-            pass
 
         if elt not in self.matrix_space():
             raise ValueError(msg)
 
         if elt not in self.matrix_space():
             raise ValueError(msg)
@@ -565,15 +806,10 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # is that we're already converting everything to long vectors,
         # and that strategy works for tuples as well.
         #
         # 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.
-        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)
+        elt = self._matrix_span.ambient_vector_space()(_all2list(elt))
 
         try:
 
         try:
-            coords = W.coordinate_vector(V(elt))
+            coords = self._matrix_span.coordinate_vector(elt)
         except ArithmeticError:  # vector is not in free module
             raise ValueError(msg)
 
         except ArithmeticError:  # vector is not in free module
             raise ValueError(msg)
 
@@ -628,7 +864,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: J = JordanSpinEJA(3)
             sage: p = J.characteristic_polynomial_of(); p
 
             sage: J = JordanSpinEJA(3)
             sage: p = J.characteristic_polynomial_of(); p
-            X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
+            X0^2 - X1^2 - X2^2 + (-2*t)*X0 + t^2
             sage: xvec = J.one().to_vector()
             sage: p(*xvec)
             t^2 - 2*t + 1
             sage: xvec = J.one().to_vector()
             sage: p(*xvec)
             t^2 - 2*t + 1
@@ -677,13 +913,13 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: J = HadamardEJA(2)
             sage: J.coordinate_polynomial_ring()
 
             sage: J = HadamardEJA(2)
             sage: J.coordinate_polynomial_ring()
-            Multivariate Polynomial Ring in X1, X2...
+            Multivariate Polynomial Ring in X0, X1...
             sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
             sage: J.coordinate_polynomial_ring()
             sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
             sage: J.coordinate_polynomial_ring()
-            Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
+            Multivariate Polynomial Ring in X0, X1, X2, X3, X4, X5...
 
         """
 
         """
-        var_names = tuple( "X%d" % z for z in range(1, self.dimension()+1) )
+        var_names = tuple( "X%d" % z for z in range(self.dimension()) )
         return PolynomialRing(self.base_ring(), var_names)
 
     def inner_product(self, x, y):
         return PolynomialRing(self.base_ring(), var_names)
 
     def inner_product(self, x, y):
@@ -705,7 +941,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         Our inner product is "associative," which means the following for
         a symmetric bilinear form::
 
         Our inner product is "associative," which means the following for
         a symmetric bilinear form::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: x,y,z = J.random_elements(3)
             sage: (x*y).inner_product(z) == y.inner_product(x*z)
             sage: J = random_eja()
             sage: x,y,z = J.random_elements(3)
             sage: (x*y).inner_product(z) == y.inner_product(x*z)
@@ -716,7 +951,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         Ensure that this is the usual inner product for the algebras
         over `R^n`::
 
         Ensure that this is the usual inner product for the algebras
         over `R^n`::
 
-            sage: set_random_seed()
             sage: J = HadamardEJA.random_instance()
             sage: x,y = J.random_elements(2)
             sage: actual = x.inner_product(y)
             sage: J = HadamardEJA.random_instance()
             sage: x,y = J.random_elements(2)
             sage: actual = x.inner_product(y)
@@ -729,7 +963,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         one). This is in Faraut and Koranyi, and also my "On the
         symmetry..." paper::
 
         one). This is in Faraut and Koranyi, and also my "On the
         symmetry..." paper::
 
-            sage: set_random_seed()
             sage: J = BilinearFormEJA.random_instance()
             sage: n = J.dimension()
             sage: x = J.random_element()
             sage: J = BilinearFormEJA.random_instance()
             sage: n = J.dimension()
             sage: x = J.random_element()
@@ -782,15 +1015,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 |
             +----++----+----+----+----+
 
         """
             +----++----+----+----+----+
 
         """
@@ -800,7 +1033,7 @@ 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.gens()[i]] + [ self.product_on_basis(i,j)
+        M += [ [self.monomial(i)] + [ self.monomial(i)*self.monomial(j)
                                     for j in range(n) ]
                for i in range(n) ]
 
                                     for j in range(n) ]
                for i in range(n) ]
 
@@ -844,7 +1077,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]
@@ -855,7 +1088,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]
@@ -909,16 +1142,16 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
             sage: J.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
+            Module of 2 by 2 matrices with entries in Algebraic Field over
+            the scalar ring Rational Field
             sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
             sage: J.matrix_space()
             sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
             sage: J.matrix_space()
-            Full MatrixSpace of 4 by 4 dense matrices over Rational Field
+            Module of 1 by 1 matrices with entries in Quaternion
+            Algebra (-1, -1) with base ring Rational Field over
+            the scalar ring Rational Field
 
         """
 
         """
-        if self.is_trivial():
-            return MatrixSpace(self.base_ring(), 0)
-        else:
-            return self.matrix_basis()[0].parent()
+        return self._matrix_space
 
 
     @cached_method
 
 
     @cached_method
@@ -937,39 +1170,37 @@ 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:
 
         The identity element acts like the identity, regardless of
         whether or not we orthonormalize::
 
 
         TESTS:
 
         The identity element acts like the identity, regardless of
         whether or not we orthonormalize::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: J.one()*x == x and x*J.one() == x
             True
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: J.one()*x == x and x*J.one() == x
             True
-            sage: A = x.subalgebra_generated_by()
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
             sage: y = A.random_element()
             sage: A.one()*y == y and y*A.one() == y
             True
 
         ::
 
             sage: y = A.random_element()
             sage: A.one()*y == y and y*A.one() == y
             True
 
         ::
 
-            sage: set_random_seed()
             sage: J = random_eja(field=QQ, orthonormalize=False)
             sage: x = J.random_element()
             sage: J.one()*x == x and x*J.one() == x
             sage: J = random_eja(field=QQ, orthonormalize=False)
             sage: x = J.random_element()
             sage: J.one()*x == x and x*J.one() == x
@@ -983,14 +1214,13 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         regardless of the base field and whether or not we
         orthonormalize::
 
         regardless of the base field and whether or not we
         orthonormalize::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: actual = J.one().operator().matrix()
             sage: expected = matrix.identity(J.base_ring(), J.dimension())
             sage: actual == expected
             True
             sage: x = J.random_element()
             sage: J = random_eja()
             sage: actual = J.one().operator().matrix()
             sage: expected = matrix.identity(J.base_ring(), J.dimension())
             sage: actual == expected
             True
             sage: x = J.random_element()
-            sage: A = x.subalgebra_generated_by()
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
             sage: actual = A.one().operator().matrix()
             sage: expected = matrix.identity(A.base_ring(), A.dimension())
             sage: actual == expected
             sage: actual = A.one().operator().matrix()
             sage: expected = matrix.identity(A.base_ring(), A.dimension())
             sage: actual == expected
@@ -998,7 +1228,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         ::
 
 
         ::
 
-            sage: set_random_seed()
             sage: J = random_eja(field=QQ, orthonormalize=False)
             sage: actual = J.one().operator().matrix()
             sage: expected = matrix.identity(J.base_ring(), J.dimension())
             sage: J = random_eja(field=QQ, orthonormalize=False)
             sage: actual = J.one().operator().matrix()
             sage: expected = matrix.identity(J.base_ring(), J.dimension())
@@ -1014,7 +1243,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         Ensure that the cached unit element (often precomputed by
         hand) agrees with the computed one::
 
         Ensure that the cached unit element (often precomputed by
         hand) agrees with the computed one::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: cached = J.one()
             sage: J.one.clear_cache()
             sage: J = random_eja()
             sage: cached = J.one()
             sage: J.one.clear_cache()
@@ -1023,7 +1251,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         ::
 
 
         ::
 
-            sage: set_random_seed()
             sage: J = random_eja(field=QQ, orthonormalize=False)
             sage: cached = J.one()
             sage: J.one.clear_cache()
             sage: J = random_eja(field=QQ, orthonormalize=False)
             sage: cached = J.one()
             sage: J.one.clear_cache()
@@ -1039,7 +1266,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         #
         # Of course, matrices aren't vectors in sage, so we have to
         # appeal to the "long vectors" isometry.
         #
         # Of course, matrices aren't vectors in sage, so we have to
         # appeal to the "long vectors" isometry.
-        oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]
+
+        V = VectorSpace(self.base_ring(), self.dimension()**2)
+        oper_vecs = [ V(g.operator().matrix().list()) for g in self.gens() ]
 
         # Now we use basic linear algebra to find the coefficients,
         # of the matrices-as-vectors-linear-combination, which should
 
         # Now we use basic linear algebra to find the coefficients,
         # of the matrices-as-vectors-linear-combination, which should
@@ -1049,7 +1278,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # We used the isometry on the left-hand side already, but we
         # still need to do it for the right-hand side. Recall that we
         # wanted something that summed to the identity matrix.
         # We used the isometry on the left-hand side already, but we
         # still need to do it for the right-hand side. Recall that we
         # wanted something that summed to the identity matrix.
-        b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) )
+        b = V( matrix.identity(self.base_ring(), self.dimension()).list() )
 
         # Now if there's an identity element in the algebra, this
         # should work. We solve on the left to avoid having to
 
         # Now if there's an identity element in the algebra, this
         # should work. We solve on the left to avoid having to
@@ -1134,7 +1363,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         Every algebra decomposes trivially with respect to its identity
         element::
 
         Every algebra decomposes trivially with respect to its identity
         element::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: J0,J5,J1 = J.peirce_decomposition(J.one())
             sage: J0.dimension() == 0 and J5.dimension() == 0
             sage: J = random_eja()
             sage: J0,J5,J1 = J.peirce_decomposition(J.one())
             sage: J0.dimension() == 0 and J5.dimension() == 0
@@ -1147,7 +1375,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         elements in the two subalgebras are the projections onto their
         respective subspaces of the superalgebra's identity element::
 
         elements in the two subalgebras are the projections onto their
         respective subspaces of the superalgebra's identity element::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: if not J.is_trivial():
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: if not J.is_trivial():
@@ -1179,7 +1406,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # 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).
         # 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 = self.subalgebra(())
+        trivial = self.subalgebra((), check_axioms=False)
         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
@@ -1221,26 +1448,13 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # For a general base ring... maybe we can trust this to do the
         # right thing? Unlikely, but.
         V = self.vector_space()
         # For a general base ring... maybe we can trust this to do the
         # right thing? Unlikely, but.
         V = self.vector_space()
-        v = V.random_element()
-
-        if self.base_ring() is AA:
-            # The "random element" method of the algebraic reals is
-            # stupid at the moment, and only returns integers between
-            # -2 and 2, inclusive:
-            #
-            #   https://trac.sagemath.org/ticket/30875
-            #
-            # Instead, we implement our own "random vector" method,
-            # and then coerce that into the algebra. We use the vector
-            # space degree here instead of the dimension because a
-            # subalgebra could (for example) be spanned by only two
-            # vectors, each with five coordinates.  We need to
-            # generate all five coordinates.
-            if thorough:
-                v *= QQbar.random_element().real()
-            else:
-                v *= QQ.random_element()
+        if self.base_ring() is AA and not thorough:
+            # Now that AA generates actually random random elements
+            # (post Trac 30875), we only need to de-thorough the
+            # randomness when asked to.
+            V = V.change_ring(QQ)
 
 
+        v = V.random_element()
         return self.from_vector(V.coordinate_vector(v))
 
     def random_elements(self, count, thorough=False):
         return self.from_vector(V.coordinate_vector(v))
 
     def random_elements(self, count, thorough=False):
@@ -1273,6 +1487,64 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                       for idx in range(count) )
 
 
                       for idx in range(count) )
 
 
+    def operator_polynomial_matrix(self):
+        r"""
+        Return the matrix of polynomials (over this algebra's
+        :meth:`coordinate_polynomial_ring`) that, when evaluated at
+        the basis coordinates of an element `x`, produces the basis
+        representation of `L_{x}`.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            ....:                                  JordanSpinEJA)
+
+        EXAMPLES::
+
+            sage: J = HadamardEJA(4)
+            sage: L_x = J.operator_polynomial_matrix()
+            sage: L_x
+            [X0  0  0  0]
+            [ 0 X1  0  0]
+            [ 0  0 X2  0]
+            [ 0  0  0 X3]
+            sage: x = J.one()
+            sage: d = zip(J.coordinate_polynomial_ring().gens(), x.to_vector())
+            sage: L_x.subs(dict(d))
+            [1 0 0 0]
+            [0 1 0 0]
+            [0 0 1 0]
+            [0 0 0 1]
+
+        ::
+
+            sage: J = JordanSpinEJA(4)
+            sage: L_x = J.operator_polynomial_matrix()
+            sage: L_x
+            [X0 X1 X2 X3]
+            [X1 X0  0  0]
+            [X2  0 X0  0]
+            [X3  0  0 X0]
+            sage: x = J.one()
+            sage: d = zip(J.coordinate_polynomial_ring().gens(), x.to_vector())
+            sage: L_x.subs(dict(d))
+            [1 0 0 0]
+            [0 1 0 0]
+            [0 0 1 0]
+            [0 0 0 1]
+
+        """
+        R = self.coordinate_polynomial_ring()
+
+        def L_x_i_j(i,j):
+            # From a result in my book, these are the entries of the
+            # basis representation of L_x.
+            return sum( v*self.monomial(k).operator().matrix()[i,j]
+                        for (k,v) in enumerate(R.gens()) )
+
+        n = self.dimension()
+        return matrix(R, n, n, L_x_i_j)
+
     @cached_method
     def _charpoly_coefficients(self):
         r"""
     @cached_method
     def _charpoly_coefficients(self):
         r"""
@@ -1288,7 +1560,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         The theory shows that these are all homogeneous polynomials of
         a known degree::
 
         The theory shows that these are all homogeneous polynomials of
         a known degree::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
             True
             sage: J = random_eja()
             sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
             True
@@ -1296,16 +1567,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         """
         n = self.dimension()
         R = self.coordinate_polynomial_ring()
         """
         n = self.dimension()
         R = self.coordinate_polynomial_ring()
-        vars = R.gens()
         F = R.fraction_field()
 
         F = R.fraction_field()
 
-        def L_x_i_j(i,j):
-            # From a result in my book, these are the entries of the
-            # basis representation of L_x.
-            return sum( vars[k]*self.gens()[k].operator().matrix()[i,j]
-                        for k in range(n) )
-
-        L_x = matrix(F, n, n, L_x_i_j)
+        L_x = self.operator_polynomial_matrix()
 
         r = None
         if self.rank.is_in_cache():
 
         r = None
         if self.rank.is_in_cache():
@@ -1386,7 +1650,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         positive integer rank, unless the algebra is trivial in
         which case its rank will be zero::
 
         positive integer rank, unless the algebra is trivial in
         which case its rank will be zero::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: r = J.rank()
             sage: r in ZZ
             sage: J = random_eja()
             sage: r = J.rank()
             sage: r in ZZ
@@ -1397,7 +1660,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         Ensure that computing the rank actually works, since the ranks
         of all simple algebras are known and will be cached by default::
 
         Ensure that computing the rank actually works, since the ranks
         of all simple algebras are known and will be cached by default::
 
-            sage: set_random_seed()    # long time
             sage: J = random_eja()     # long time
             sage: cached = J.rank()    # long time
             sage: J.rank.clear_cache() # long time
             sage: J = random_eja()     # long time
             sage: cached = J.rank()    # long time
             sage: J.rank.clear_cache() # long time
@@ -1412,8 +1674,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         r"""
         Create a subalgebra of this algebra from the given basis.
         """
         r"""
         Create a subalgebra of this algebra from the given basis.
         """
-        from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
-        return FiniteDimensionalEJASubalgebra(self, basis, **kwargs)
+        from mjo.eja.eja_subalgebra import EJASubalgebra
+        return EJASubalgebra(self, basis, **kwargs)
 
 
     def vector_space(self):
 
 
     def vector_space(self):
@@ -1435,9 +1697,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
 
 
 
 
 
-class RationalBasisEJA(FiniteDimensionalEJA):
+class RationalBasisEJA(EJA):
     r"""
     r"""
-    New class for algebras whose supplied basis elements have all rational entries.
+    Algebras whose supplied basis elements have all rational entries.
 
     SETUP::
 
 
     SETUP::
 
@@ -1468,9 +1730,20 @@ class RationalBasisEJA(FiniteDimensionalEJA):
         if check_field:
             # Abuse the check_field parameter to check that the entries of
             # out basis (in ambient coordinates) are in the field QQ.
         if check_field:
             # Abuse the check_field parameter to check that the entries of
             # out basis (in ambient coordinates) are in the field QQ.
-            if not all( all(b_i in QQ for b_i in b.list()) for b in basis ):
+            # Use _all2list to get the vector coordinates of octonion
+            # entries and not the octonions themselves (which are not
+            # rational).
+            if not all( all(b_i in QQ for b_i in _all2list(b))
+                        for b in basis ):
                 raise TypeError("basis not rational")
 
                 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
@@ -1479,21 +1752,25 @@ class RationalBasisEJA(FiniteDimensionalEJA):
             # Note: the same Jordan and inner-products work here,
             # because they are necessarily defined with respect to
             # ambient coordinates and not any particular basis.
             # Note: the same Jordan and inner-products work here,
             # because they are necessarily defined with respect to
             # ambient coordinates and not any particular basis.
-            self._rational_algebra = FiniteDimensionalEJA(
+            self._rational_algebra = EJA(
                                        basis,
                                        jordan_product,
                                        inner_product,
                                        field=QQ,
                                        basis,
                                        jordan_product,
                                        inner_product,
                                        field=QQ,
+                                       matrix_space=self.matrix_space(),
+                                       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)
+    def rational_algebra(self):
+        # Using None as a flag here (rather than just assigning "self"
+        # to self._rational_algebra by default) feels a little bit
+        # more sane to me in a garbage-collected environment.
+        if self._rational_algebra is None:
+            return self
+        else:
+            return self._rational_algebra
 
     @cached_method
     def _charpoly_coefficients(self):
 
     @cached_method
     def _charpoly_coefficients(self):
@@ -1511,7 +1788,7 @@ class RationalBasisEJA(FiniteDimensionalEJA):
 
             sage: J = JordanSpinEJA(3)
             sage: J._charpoly_coefficients()
 
             sage: J = JordanSpinEJA(3)
             sage: J._charpoly_coefficients()
-            (X1^2 - X2^2 - X3^2, -2*X1)
+            (X0^2 - X1^2 - X2^2, -2*X0)
             sage: a0 = J._charpoly_coefficients()[0]
             sage: J.base_ring()
             Algebraic Real Field
             sage: a0 = J._charpoly_coefficients()[0]
             sage: J.base_ring()
             Algebraic Real Field
@@ -1519,28 +1796,17 @@ class RationalBasisEJA(FiniteDimensionalEJA):
             Algebraic Real Field
 
         """
             Algebraic Real Field
 
         """
-        if self._rational_algebra is None:
-            # There's no need to construct *another* algebra over the
-            # rationals if this one is already over the
-            # rationals. Likewise, if we never orthonormalized our
-            # basis, we might as well just use the given one.
+        if self.rational_algebra() is self:
+            # Bypass the hijinks if they won't benefit us.
             return super()._charpoly_coefficients()
 
             return super()._charpoly_coefficients()
 
-        # Do the computation over the rationals. The answer will be
-        # the same, because all we've done is a change of basis.
-        # Then, change back from QQ to our real base ring
+        # Do the computation over the rationals.
         a = ( a_i.change_ring(self.base_ring())
         a = ( a_i.change_ring(self.base_ring())
-              for a_i in self._rational_algebra._charpoly_coefficients() )
+              for a_i in self.rational_algebra()._charpoly_coefficients() )
 
 
-        if self._deortho_matrix is None:
-            # This can happen if our base ring was, say, AA and we
-            # chose not to (or didn't need to) orthonormalize. It's
-            # still faster to do the computations over QQ even if
-            # the numbers in the boxes stay the same.
-            return tuple(a)
-
-        # Otherwise, convert the coordinate variables back to the
-        # deorthonormalized ones.
+        # Convert our coordinate variables into deorthonormalized ones
+        # and substitute them into the deorthonormalized charpoly
+        # coefficients.
         R = self.coordinate_polynomial_ring()
         from sage.modules.free_module_element import vector
         X = vector(R, R.gens())
         R = self.coordinate_polynomial_ring()
         from sage.modules.free_module_element import vector
         X = vector(R, R.gens())
@@ -1549,7 +1815,7 @@ class RationalBasisEJA(FiniteDimensionalEJA):
         subs_dict = { X[i]: BX[i] for i in range(len(X)) }
         return tuple( a_i.subs(subs_dict) for a_i in a )
 
         subs_dict = { X[i]: BX[i] for i in range(len(X)) }
         return tuple( a_i.subs(subs_dict) for a_i in a )
 
-class ConcreteEJA(RationalBasisEJA):
+class ConcreteEJA(EJA):
     r"""
     A class for the Euclidean Jordan algebras that we know by name.
 
     r"""
     A class for the Euclidean Jordan algebras that we know by name.
 
@@ -1567,7 +1833,6 @@ class ConcreteEJA(RationalBasisEJA):
     Our basis is normalized with respect to the algebra's inner
     product, unless we specify otherwise::
 
     Our basis is normalized with respect to the algebra's inner
     product, unless we specify otherwise::
 
-        sage: set_random_seed()
         sage: J = ConcreteEJA.random_instance()
         sage: all( b.norm() == 1 for b in J.gens() )
         True
         sage: J = ConcreteEJA.random_instance()
         sage: all( b.norm() == 1 for b in J.gens() )
         True
@@ -1578,7 +1843,6 @@ class ConcreteEJA(RationalBasisEJA):
     natural->EJA basis representation is an isometry and within the
     EJA the operator is self-adjoint by the Jordan axiom::
 
     natural->EJA basis representation is an isometry and within the
     EJA the operator is self-adjoint by the Jordan axiom::
 
-        sage: set_random_seed()
         sage: J = ConcreteEJA.random_instance()
         sage: x = J.random_element()
         sage: x.operator().is_self_adjoint()
         sage: J = ConcreteEJA.random_instance()
         sage: x = J.random_element()
         sage: x.operator().is_self_adjoint()
@@ -1586,27 +1850,62 @@ class ConcreteEJA(RationalBasisEJA):
     """
 
     @staticmethod
     """
 
     @staticmethod
-    def _max_random_instance_size():
+    def _max_random_instance_dimension():
+        r"""
+        The maximum dimension of any random instance. Ten dimensions seems
+        to be about the point where everything takes a turn for the
+        worse. And dimension ten (but not nine) allows the 4-by-4 real
+        Hermitian matrices, the 2-by-2 quaternion Hermitian matrices,
+        and the 2-by-2 octonion Hermitian matrices.
+        """
+        return 10
+
+    @staticmethod
+    def _max_random_instance_size(max_dimension):
         """
         Return an integer "size" that is an upper bound on the size of
         """
         Return an integer "size" that is an upper bound on the size of
-        this algebra when it is used in a random test
-        case. Unfortunately, the term "size" is ambiguous -- when
-        dealing with `R^n` under either the Hadamard or Jordan spin
-        product, the "size" refers to the dimension `n`. When dealing
-        with a matrix algebra (real symmetric or complex/quaternion
-        Hermitian), it refers to the size of the matrix, which is far
-        less than the dimension of the underlying vector space.
+        this algebra when it is used in a random test case. This size
+        (which can be passed to the algebra's constructor) is itself
+        based on the ``max_dimension`` parameter.
 
         This method must be implemented in each subclass.
         """
         raise NotImplementedError
 
     @classmethod
 
         This method must be implemented in each subclass.
         """
         raise NotImplementedError
 
     @classmethod
-    def random_instance(cls, *args, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         """
-        Return a random instance of this type of algebra.
+        Return a random instance of this type of algebra whose dimension
+        is less than or equal to the lesser of ``max_dimension`` and
+        the value returned by ``_max_random_instance_dimension()``. If
+        the dimension bound is omitted, then only the
+        ``_max_random_instance_dimension()`` is used as a bound.
 
         This method should be implemented in each subclass.
 
         This method should be implemented in each subclass.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import ConcreteEJA
+
+        TESTS:
+
+        Both the class bound and the ``max_dimension`` argument are upper
+        bounds on the dimension of the algebra returned::
+
+            sage: from sage.misc.prandom import choice
+            sage: eja_class = choice(ConcreteEJA.__subclasses__())
+            sage: class_max_d = eja_class._max_random_instance_dimension()
+            sage: J = eja_class.random_instance(max_dimension=20,
+            ....:                               field=QQ,
+            ....:                               orthonormalize=False)
+            sage: J.dimension() <= class_max_d
+            True
+            sage: J = eja_class.random_instance(max_dimension=2,
+            ....:                               field=QQ,
+            ....:                               orthonormalize=False)
+            sage: J.dimension() <= 2
+            True
+
         """
         from sage.misc.prandom import choice
         eja_class = choice(cls.__subclasses__())
         """
         from sage.misc.prandom import choice
         eja_class = choice(cls.__subclasses__())
@@ -1614,137 +1913,162 @@ class ConcreteEJA(RationalBasisEJA):
         # These all bubble up to the RationalBasisEJA superclass
         # constructor, so any (kw)args valid there are also valid
         # here.
         # These all bubble up to the RationalBasisEJA superclass
         # constructor, so any (kw)args valid there are also valid
         # here.
-        return eja_class.random_instance(*args, **kwargs)
+        return eja_class.random_instance(max_dimension, *args, **kwargs)
 
 
 
 
-class MatrixEJA:
+class HermitianMatrixEJA(EJA):
     @staticmethod
     @staticmethod
-    def dimension_over_reals():
-        r"""
-        The dimension of this matrix's base ring over the reals.
+    def _denormalized_basis(A):
+        """
+        Returns a basis for the given Hermitian matrix space.
 
 
-        The reals are dimension one over themselves, obviously; that's
-        just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
-        have dimension two. Finally, the quaternions have dimension
-        four over the reals.
+        Why do we embed these? Basically, because all of numerical linear
+        algebra assumes that you're working with vectors consisting of `n`
+        entries from a field and scalars from the same field. There's no way
+        to tell SageMath that (for example) the vectors contain complex
+        numbers, while the scalar field is real.
 
 
-        This is used to determine the size of the matrix returned from
-        :meth:`real_embed`, among other things.
-        """
-        raise NotImplementedError
+        SETUP::
 
 
-    @classmethod
-    def real_embed(cls,M):
-        """
-        Embed the matrix ``M`` into a space of real matrices.
+            sage: from mjo.hurwitz import (ComplexMatrixAlgebra,
+            ....:                          QuaternionMatrixAlgebra,
+            ....:                          OctonionMatrixAlgebra)
+            sage: from mjo.eja.eja_algebra import HermitianMatrixEJA
 
 
-        The matrix ``M`` can have entries in any field at the moment:
-        the real numbers, complex numbers, or quaternions. And although
-        they are not a field, we can probably support octonions at some
-        point, too. This function returns a real matrix that "acts like"
-        the original with respect to matrix multiplication; i.e.
+        TESTS::
 
 
-          real_embed(M*N) = real_embed(M)*real_embed(N)
+            sage: n = ZZ.random_element(1,5)
+            sage: A = MatrixSpace(QQ, n)
+            sage: B = HermitianMatrixEJA._denormalized_basis(A)
+            sage: all( M.is_hermitian() for M in  B)
+            True
 
 
-        """
-        if M.ncols() != M.nrows():
-            raise ValueError("the matrix 'M' must be square")
-        return M
+        ::
 
 
+            sage: n = ZZ.random_element(1,5)
+            sage: A = ComplexMatrixAlgebra(n, scalars=QQ)
+            sage: B = HermitianMatrixEJA._denormalized_basis(A)
+            sage: all( M.is_hermitian() for M in  B)
+            True
+
+        ::
+
+            sage: n = ZZ.random_element(1,5)
+            sage: A = QuaternionMatrixAlgebra(n, scalars=QQ)
+            sage: B = HermitianMatrixEJA._denormalized_basis(A)
+            sage: all( M.is_hermitian() for M in B )
+            True
+
+        ::
+
+            sage: n = ZZ.random_element(1,5)
+            sage: A = OctonionMatrixAlgebra(n, scalars=QQ)
+            sage: B = HermitianMatrixEJA._denormalized_basis(A)
+            sage: all( M.is_hermitian() for M in B )
+            True
 
 
-    @classmethod
-    def real_unembed(cls,M):
-        """
-        The inverse of :meth:`real_embed`.
         """
         """
-        if M.ncols() != M.nrows():
-            raise ValueError("the matrix 'M' must be square")
-        if not ZZ(M.nrows()).mod(cls.dimension_over_reals()).is_zero():
-            raise ValueError("the matrix 'M' must be a real embedding")
-        return M
+        # These work for real MatrixSpace, whose monomials only have
+        # two coordinates (because the last one would always be "1").
+        es = A.base_ring().gens()
+        gen = lambda A,m: A.monomial(m[:2])
+
+        if hasattr(A, 'entry_algebra_gens'):
+            # We've got a MatrixAlgebra, and its monomials will have
+            # three coordinates.
+            es = A.entry_algebra_gens()
+            gen = lambda A,m: A.monomial(m)
+
+        basis = []
+        for i in range(A.nrows()):
+            for j in range(i+1):
+                if i == j:
+                    E_ii = gen(A, (i,j,es[0]))
+                    basis.append(E_ii)
+                else:
+                    for e in es:
+                        E_ij  = gen(A, (i,j,e))
+                        E_ij += E_ij.conjugate_transpose()
+                        basis.append(E_ij)
+
+        return tuple( basis )
 
     @staticmethod
     def jordan_product(X,Y):
         return (X*Y + Y*X)/2
 
 
     @staticmethod
     def jordan_product(X,Y):
         return (X*Y + Y*X)/2
 
-    @classmethod
-    def trace_inner_product(cls,X,Y):
+    @staticmethod
+    def trace_inner_product(X,Y):
         r"""
         r"""
-        Compute the trace inner-product of two real-embeddings.
+        A trace inner-product for matrices that aren't embedded in the
+        reals. It takes MATRICES as arguments, not EJA elements.
 
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
             ....:                                  ComplexHermitianEJA,
 
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
             ....:                                  ComplexHermitianEJA,
-            ....:                                  QuaternionHermitianEJA)
+            ....:                                  QuaternionHermitianEJA,
+            ....:                                  OctonionHermitianEJA)
 
         EXAMPLES::
 
 
         EXAMPLES::
 
-        This gives the same answer as it would if we computed the trace
-        from the unembedded (original) matrices::
+            sage: J = RealSymmetricEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2
 
 
-            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: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2
 
         ::
 
 
         ::
 
-            sage: set_random_seed()
-            sage: J = ComplexHermitianEJA.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().real()
-            sage: actual = J.trace_inner_product(Xe,Ye)
-            sage: actual == expected
-            True
+            sage: J = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2
 
         ::
 
 
         ::
 
-            sage: set_random_seed()
-            sage: J = QuaternionHermitianEJA.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().coefficient_tuple()[0]
-            sage: actual = J.trace_inner_product(Xe,Ye)
-            sage: actual == expected
-            True
+            sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2
 
         """
 
         """
-        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.
+        tr = (X*Y).trace()
+        if hasattr(tr, 'coefficient'):
+            # Works for octonions, and has to come first because they
+            # also have a "real()" method that doesn't return an
+            # element of the scalar ring.
+            return tr.coefficient(0)
+        elif hasattr(tr, 'coefficient_tuple'):
+            # Works for quaternions.
             return tr.coefficient_tuple()[0]
 
             return tr.coefficient_tuple()[0]
 
+        # Works for real and complex numbers.
+        return tr.real()
 
 
-class RealMatrixEJA(MatrixEJA):
-    @staticmethod
-    def dimension_over_reals():
-        return 1
 
 
+    def __init__(self, matrix_space, **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().__init__(self._denormalized_basis(matrix_space),
+                         self.jordan_product,
+                         self.trace_inner_product,
+                         field=matrix_space.base_ring(),
+                         matrix_space=matrix_space,
+                         **kwargs)
+
+        self.rank.set_cache(matrix_space.nrows())
+        self.one.set_cache( self(matrix_space.one()) )
 
 
-class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
+class RealSymmetricEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA):
     """
     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
@@ -1757,19 +2081,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
 
@@ -1777,16 +2101,14 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
 
     The dimension of this algebra is `(n^2 + n) / 2`::
 
 
     The dimension of this algebra is `(n^2 + n) / 2`::
 
-        sage: set_random_seed()
-        sage: n_max = RealSymmetricEJA._max_random_instance_size()
-        sage: n = ZZ.random_element(1, n_max)
+        sage: d = RealSymmetricEJA._max_random_instance_dimension()
+        sage: n = RealSymmetricEJA._max_random_instance_size(d)
         sage: J = RealSymmetricEJA(n)
         sage: J.dimension() == (n^2 + n)/2
         True
 
     The Jordan multiplication is what we think it is::
 
         sage: J = RealSymmetricEJA(n)
         sage: J.dimension() == (n^2 + n)/2
         True
 
     The Jordan multiplication is what we think it is::
 
-        sage: set_random_seed()
         sage: J = RealSymmetricEJA.random_instance()
         sage: x,y = J.random_elements(2)
         sage: actual = (x*y).to_matrix()
         sage: J = RealSymmetricEJA.random_instance()
         sage: x,y = J.random_elements(2)
         sage: actual = (x*y).to_matrix()
@@ -1809,218 +2131,36 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """
-    @classmethod
-    def _denormalized_basis(cls, n):
-        """
-        Return a basis for the space of real symmetric n-by-n matrices.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import RealSymmetricEJA
-
-        TESTS::
-
-            sage: set_random_seed()
-            sage: n = ZZ.random_element(1,5)
-            sage: B = RealSymmetricEJA._denormalized_basis(n)
-            sage: all( M.is_symmetric() for M in  B)
-            True
-
-        """
-        # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
-        # coordinates.
-        S = []
-        for i in range(n):
-            for j in range(i+1):
-                Eij = matrix(ZZ, n, lambda k,l: k==i and l==j)
-                if i == j:
-                    Sij = Eij
-                else:
-                    Sij = Eij + Eij.transpose()
-                S.append(Sij)
-        return tuple(S)
-
-
     @staticmethod
     @staticmethod
-    def _max_random_instance_size():
-        return 4 # Dimension 10
+    def _max_random_instance_size(max_dimension):
+        # Obtained by solving d = (n^2 + n)/2.
+        # The ZZ-int-ZZ thing is just "floor."
+        return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/2 - 1/2))
 
     @classmethod
 
     @classmethod
-    def random_instance(cls, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this type of algebra.
         """
         """
         Return a random instance of this type of algebra.
         """
-        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
         return cls(n, **kwargs)
 
-    def __init__(self, n, **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)
-
-        # 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)
-        self.one.set_cache(self(idV))
+    def __init__(self, n, field=AA, **kwargs):
+        A = MatrixSpace(field, n)
+        super().__init__(A, **kwargs)
 
 
+        from mjo.eja.eja_cache import real_symmetric_eja_coeffs
+        a = real_symmetric_eja_coeffs(self)
+        if a is not None:
+            self.rational_algebra()._charpoly_coefficients.set_cache(a)
 
 
 
 
-class ComplexMatrixEJA(MatrixEJA):
-    # A manual dictionary-cache for the complex_extension() method,
-    # since apparently @classmethods can't also be @cached_methods.
-    _complex_extension = {}
 
 
-    @classmethod
-    def complex_extension(cls,field):
-        r"""
-        The complex field that we embed/unembed, as an extension
-        of the given ``field``.
-        """
-        if field in cls._complex_extension:
-            return cls._complex_extension[field]
-
-        # Sage doesn't know how to adjoin the complex "i" (the root of
-        # x^2 + 1) to a field in a general way. Here, we just enumerate
-        # all of the cases that I have cared to support so far.
-        if field is AA:
-            # Sage doesn't know how to embed AA into QQbar, i.e. how
-            # to adjoin sqrt(-1) to AA.
-            F = QQbar
-        elif not field.is_exact():
-            # RDF or RR
-            F = field.complex_field()
-        else:
-            # Works for QQ and... maybe some other fields.
-            R = PolynomialRing(field, 'z')
-            z = R.gen()
-            F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
-
-        cls._complex_extension[field] = F
-        return F
-
-    @staticmethod
-    def dimension_over_reals():
-        return 2
-
-    @classmethod
-    def real_embed(cls,M):
-        """
-        Embed the n-by-n complex matrix ``M`` into the space of real
-        matrices of size 2n-by-2n via the map the sends each entry `z = a +
-        bi` to the block matrix ``[[a,b],[-b,a]]``.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
-
-        EXAMPLES::
-
-            sage: F = QuadraticField(-1, 'I')
-            sage: x1 = F(4 - 2*i)
-            sage: x2 = F(1 + 2*i)
-            sage: x3 = F(-i)
-            sage: x4 = F(6)
-            sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
-            sage: ComplexMatrixEJA.real_embed(M)
-            [ 4 -2| 1  2]
-            [ 2  4|-2  1]
-            [-----+-----]
-            [ 0 -1| 6  0]
-            [ 1  0| 0  6]
-
-        TESTS:
-
-        Embedding is a homomorphism (isomorphism, in fact)::
-
-            sage: set_random_seed()
-            sage: n = ZZ.random_element(3)
-            sage: F = QuadraticField(-1, 'I')
-            sage: X = random_matrix(F, n)
-            sage: Y = random_matrix(F, n)
-            sage: Xe = ComplexMatrixEJA.real_embed(X)
-            sage: Ye = ComplexMatrixEJA.real_embed(Y)
-            sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
-            sage: Xe*Ye == XYe
-            True
-
-        """
-        super(ComplexMatrixEJA,cls).real_embed(M)
-        n = M.nrows()
-
-        # We don't need any adjoined elements...
-        field = M.base_ring().base_ring()
-
-        blocks = []
-        for z in M.list():
-            a = z.real()
-            b = z.imag()
-            blocks.append(matrix(field, 2, [ [ a, b],
-                                             [-b, a] ]))
-
-        return matrix.block(field, n, blocks)
-
-
-    @classmethod
-    def real_unembed(cls,M):
-        """
-        The inverse of _embed_complex_matrix().
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
-
-        EXAMPLES::
-
-            sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
-            ....:                 [-2,  1,  -4,  3],
-            ....:                 [ 9,  10, 11, 12],
-            ....:                 [-10, 9, -12, 11] ])
-            sage: ComplexMatrixEJA.real_unembed(A)
-            [  2*I + 1   4*I + 3]
-            [ 10*I + 9 12*I + 11]
-
-        TESTS:
-
-        Unembedding is the inverse of embedding::
-
-            sage: set_random_seed()
-            sage: F = QuadraticField(-1, 'I')
-            sage: M = random_matrix(F, 3)
-            sage: Me = ComplexMatrixEJA.real_embed(M)
-            sage: ComplexMatrixEJA.real_unembed(Me) == M
-            True
-
-        """
-        super(ComplexMatrixEJA,cls).real_unembed(M)
-        n = ZZ(M.nrows())
-        d = cls.dimension_over_reals()
-        F = cls.complex_extension(M.base_ring())
-        i = F.gen()
-
-        # Go top-left to bottom-right (reading order), converting every
-        # 2-by-2 block we see to a single complex element.
-        elements = []
-        for k in range(n/d):
-            for j in range(n/d):
-                submat = M[d*k:d*k+d,d*j:d*j+d]
-                if submat[0,0] != submat[1,1]:
-                    raise ValueError('bad on-diagonal submatrix')
-                if submat[0,1] != -submat[1,0]:
-                    raise ValueError('bad off-diagonal submatrix')
-                z = submat[0,0] + submat[0,1]*i
-                elements.append(z)
-
-        return matrix(F, n/d, elements)
-
-
-class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
+class ComplexHermitianEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA):
     """
     The rank-n simple EJA consisting of complex Hermitian n-by-n
     matrices over the real numbers, the usual symmetric Jordan product,
     """
     The rank-n simple EJA consisting of complex Hermitian n-by-n
     matrices over the real numbers, the usual symmetric Jordan product,
@@ -2033,28 +2173,32 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
 
     EXAMPLES:
 
 
     EXAMPLES:
 
-    In theory, our "field" can be any subfield of the reals::
+    In theory, our "field" can be any subfield of the reals, but we
+    can't use inexact real fields at the moment because SageMath
+    doesn't know how to convert their elements into complex numbers,
+    or even into algebraic reals::
 
 
-        sage: ComplexHermitianEJA(2, field=RDF)
-        Euclidean Jordan algebra of dimension 4 over Real Double Field
-        sage: ComplexHermitianEJA(2, field=RR)
-        Euclidean Jordan algebra of dimension 4 over Real Field with
-        53 bits of precision
+        sage: QQbar(RDF(1))
+        Traceback (most recent call last):
+        ...
+        TypeError: Illegal initializer for algebraic number
+        sage: AA(RR(1))
+        Traceback (most recent call last):
+        ...
+        TypeError: Illegal initializer for algebraic number
 
     TESTS:
 
     The dimension of this algebra is `n^2`::
 
 
     TESTS:
 
     The dimension of this algebra is `n^2`::
 
-        sage: set_random_seed()
-        sage: n_max = ComplexHermitianEJA._max_random_instance_size()
-        sage: n = ZZ.random_element(1, n_max)
+        sage: d = ComplexHermitianEJA._max_random_instance_dimension()
+        sage: n = ComplexHermitianEJA._max_random_instance_size(d)
         sage: J = ComplexHermitianEJA(n)
         sage: J.dimension() == n^2
         True
 
     The Jordan multiplication is what we think it is::
 
         sage: J = ComplexHermitianEJA(n)
         sage: J.dimension() == n^2
         True
 
     The Jordan multiplication is what we think it is::
 
-        sage: set_random_seed()
         sage: J = ComplexHermitianEJA.random_instance()
         sage: x,y = J.random_elements(2)
         sage: actual = (x*y).to_matrix()
         sage: J = ComplexHermitianEJA.random_instance()
         sage: x,y = J.random_elements(2)
         sage: actual = (x*y).to_matrix()
@@ -2077,247 +2221,36 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """
+    def __init__(self, n, field=AA, **kwargs):
+        from mjo.hurwitz import ComplexMatrixAlgebra
+        A = ComplexMatrixAlgebra(n, scalars=field)
+        super().__init__(A, **kwargs)
 
 
-    @classmethod
-    def _denormalized_basis(cls, n):
-        """
-        Returns a basis for the space of complex Hermitian n-by-n matrices.
-
-        Why do we embed these? Basically, because all of numerical linear
-        algebra assumes that you're working with vectors consisting of `n`
-        entries from a field and scalars from the same field. There's no way
-        to tell SageMath that (for example) the vectors contain complex
-        numbers, while the scalar field is real.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
-
-        TESTS::
-
-            sage: set_random_seed()
-            sage: n = ZZ.random_element(1,5)
-            sage: B = ComplexHermitianEJA._denormalized_basis(n)
-            sage: all( M.is_symmetric() for M in  B)
-            True
-
-        """
-        field = ZZ
-        R = PolynomialRing(field, 'z')
-        z = R.gen()
-        F = field.extension(z**2 + 1, 'I')
-        I = F.gen(1)
-
-        # This is like the symmetric case, but we need to be careful:
-        #
-        #   * We want conjugate-symmetry, not just symmetry.
-        #   * The diagonal will (as a result) be real.
-        #
-        S = []
-        Eij = matrix.zero(F,n)
-        for i in range(n):
-            for j in range(i+1):
-                # "build" E_ij
-                Eij[i,j] = 1
-                if i == j:
-                    Sij = cls.real_embed(Eij)
-                    S.append(Sij)
-                else:
-                    # The second one has a minus because it's conjugated.
-                    Eij[j,i] = 1 # Eij = Eij + Eij.transpose()
-                    Sij_real = cls.real_embed(Eij)
-                    S.append(Sij_real)
-                    # Eij = I*Eij - I*Eij.transpose()
-                    Eij[i,j] = I
-                    Eij[j,i] = -I
-                    Sij_imag = cls.real_embed(Eij)
-                    S.append(Sij_imag)
-                    Eij[j,i] = 0
-                # "erase" E_ij
-                Eij[i,j] = 0
-
-        # Since we embedded these, we can drop back to the "field" that we
-        # started with instead of the complex extension "F".
-        return tuple( s.change_ring(field) for s in S )
-
-
-    def __init__(self, n, **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)
-        # 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)
-        self.one.set_cache(self(idV))
+        from mjo.eja.eja_cache import complex_hermitian_eja_coeffs
+        a = complex_hermitian_eja_coeffs(self)
+        if a is not None:
+            self.rational_algebra()._charpoly_coefficients.set_cache(a)
 
     @staticmethod
 
     @staticmethod
-    def _max_random_instance_size():
-        return 3 # Dimension 9
+    def _max_random_instance_size(max_dimension):
+        # Obtained by solving d = n^2.
+        # The ZZ-int-ZZ thing is just "floor."
+        return ZZ(int(ZZ(max_dimension).sqrt()))
 
     @classmethod
 
     @classmethod
-    def random_instance(cls, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this type of algebra.
         """
         """
         Return a random instance of this type of algebra.
         """
-        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
         return cls(n, **kwargs)
 
-class QuaternionMatrixEJA(MatrixEJA):
 
 
-    # A manual dictionary-cache for the quaternion_extension() method,
-    # since apparently @classmethods can't also be @cached_methods.
-    _quaternion_extension = {}
-
-    @classmethod
-    def quaternion_extension(cls,field):
-        r"""
-        The quaternion field that we embed/unembed, as an extension
-        of the given ``field``.
-        """
-        if field in cls._quaternion_extension:
-            return cls._quaternion_extension[field]
-
-        Q = QuaternionAlgebra(field,-1,-1)
-
-        cls._quaternion_extension[field] = Q
-        return Q
-
-    @staticmethod
-    def dimension_over_reals():
-        return 4
-
-    @classmethod
-    def real_embed(cls,M):
-        """
-        Embed the n-by-n quaternion matrix ``M`` into the space of real
-        matrices of size 4n-by-4n by first sending each quaternion entry `z
-        = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
-        c+di],[-c + di, a-bi]]`, and then embedding those into a real
-        matrix.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
-
-        EXAMPLES::
-
-            sage: Q = QuaternionAlgebra(QQ,-1,-1)
-            sage: i,j,k = Q.gens()
-            sage: x = 1 + 2*i + 3*j + 4*k
-            sage: M = matrix(Q, 1, [[x]])
-            sage: QuaternionMatrixEJA.real_embed(M)
-            [ 1  2  3  4]
-            [-2  1 -4  3]
-            [-3  4  1 -2]
-            [-4 -3  2  1]
-
-        Embedding is a homomorphism (isomorphism, in fact)::
-
-            sage: set_random_seed()
-            sage: n = ZZ.random_element(2)
-            sage: Q = QuaternionAlgebra(QQ,-1,-1)
-            sage: X = random_matrix(Q, n)
-            sage: Y = random_matrix(Q, n)
-            sage: Xe = QuaternionMatrixEJA.real_embed(X)
-            sage: Ye = QuaternionMatrixEJA.real_embed(Y)
-            sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
-            sage: Xe*Ye == XYe
-            True
-
-        """
-        super(QuaternionMatrixEJA,cls).real_embed(M)
-        quaternions = M.base_ring()
-        n = M.nrows()
-
-        F = QuadraticField(-1, 'I')
-        i = F.gen()
-
-        blocks = []
-        for z in M.list():
-            t = z.coefficient_tuple()
-            a = t[0]
-            b = t[1]
-            c = t[2]
-            d = t[3]
-            cplxM = matrix(F, 2, [[ a + b*i, c + d*i],
-                                 [-c + d*i, a - b*i]])
-            realM = ComplexMatrixEJA.real_embed(cplxM)
-            blocks.append(realM)
-
-        # We should have real entries by now, so use the realest field
-        # we've got for the return value.
-        return matrix.block(quaternions.base_ring(), n, blocks)
-
-
-
-    @classmethod
-    def real_unembed(cls,M):
-        """
-        The inverse of _embed_quaternion_matrix().
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
-
-        EXAMPLES::
-
-            sage: M = matrix(QQ, [[ 1,  2,  3,  4],
-            ....:                 [-2,  1, -4,  3],
-            ....:                 [-3,  4,  1, -2],
-            ....:                 [-4, -3,  2,  1]])
-            sage: QuaternionMatrixEJA.real_unembed(M)
-            [1 + 2*i + 3*j + 4*k]
-
-        TESTS:
-
-        Unembedding is the inverse of embedding::
-
-            sage: set_random_seed()
-            sage: Q = QuaternionAlgebra(QQ, -1, -1)
-            sage: M = random_matrix(Q, 3)
-            sage: Me = QuaternionMatrixEJA.real_embed(M)
-            sage: QuaternionMatrixEJA.real_unembed(Me) == M
-            True
-
-        """
-        super(QuaternionMatrixEJA,cls).real_unembed(M)
-        n = ZZ(M.nrows())
-        d = cls.dimension_over_reals()
-
-        # Use the base ring of the matrix to ensure that its entries can be
-        # multiplied by elements of the quaternion algebra.
-        Q = cls.quaternion_extension(M.base_ring())
-        i,j,k = Q.gens()
-
-        # Go top-left to bottom-right (reading order), converting every
-        # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
-        # quaternion block.
-        elements = []
-        for l in range(n/d):
-            for m in range(n/d):
-                submat = ComplexMatrixEJA.real_unembed(
-                    M[d*l:d*l+d,d*m:d*m+d] )
-                if submat[0,0] != submat[1,1].conjugate():
-                    raise ValueError('bad on-diagonal submatrix')
-                if submat[0,1] != -submat[1,0].conjugate():
-                    raise ValueError('bad off-diagonal submatrix')
-                z  = submat[0,0].real()
-                z += submat[0,0].imag()*i
-                z += submat[0,1].real()*j
-                z += submat[0,1].imag()*k
-                elements.append(z)
-
-        return matrix(Q, n/d, elements)
-
-
-class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
+class QuaternionHermitianEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA):
     r"""
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the
     r"""
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the
@@ -2332,9 +2265,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
 
@@ -2342,16 +2275,14 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
 
     The dimension of this algebra is `2*n^2 - n`::
 
 
     The dimension of this algebra is `2*n^2 - n`::
 
-        sage: set_random_seed()
-        sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
-        sage: n = ZZ.random_element(1, n_max)
+        sage: d = QuaternionHermitianEJA._max_random_instance_dimension()
+        sage: n = QuaternionHermitianEJA._max_random_instance_size(d)
         sage: J = QuaternionHermitianEJA(n)
         sage: J.dimension() == 2*(n^2) - n
         True
 
     The Jordan multiplication is what we think it is::
 
         sage: J = QuaternionHermitianEJA(n)
         sage: J.dimension() == 2*(n^2) - n
         True
 
     The Jordan multiplication is what we think it is::
 
-        sage: set_random_seed()
         sage: J = QuaternionHermitianEJA.random_instance()
         sage: x,y = J.random_elements(2)
         sage: actual = (x*y).to_matrix()
         sage: J = QuaternionHermitianEJA.random_instance()
         sage: x,y = J.random_elements(2)
         sage: actual = (x*y).to_matrix()
@@ -2374,120 +2305,203 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """
-    @classmethod
-    def _denormalized_basis(cls, n):
-        """
-        Returns a basis for the space of quaternion Hermitian n-by-n matrices.
+    def __init__(self, n, field=AA, **kwargs):
+        from mjo.hurwitz import QuaternionMatrixAlgebra
+        A = QuaternionMatrixAlgebra(n, scalars=field)
+        super().__init__(A, **kwargs)
 
 
-        Why do we embed these? Basically, because all of numerical
-        linear algebra assumes that you're working with vectors consisting
-        of `n` entries from a field and scalars from the same field. There's
-        no way to tell SageMath that (for example) the vectors contain
-        complex numbers, while the scalar field is real.
+        from mjo.eja.eja_cache import quaternion_hermitian_eja_coeffs
+        a = quaternion_hermitian_eja_coeffs(self)
+        if a is not None:
+            self.rational_algebra()._charpoly_coefficients.set_cache(a)
 
 
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
 
 
-        TESTS::
 
 
-            sage: set_random_seed()
-            sage: n = ZZ.random_element(1,5)
-            sage: B = QuaternionHermitianEJA._denormalized_basis(n)
-            sage: all( M.is_symmetric() for M in B )
-            True
+    @staticmethod
+    def _max_random_instance_size(max_dimension):
+        r"""
+        The maximum rank of a random QuaternionHermitianEJA.
+        """
+        # Obtained by solving d = 2n^2 - n.
+        # The ZZ-int-ZZ thing is just "floor."
+        return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/4 + 1/4))
 
 
+    @classmethod
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
+        """
+        Return a random instance of this type of algebra.
         """
         """
-        field = ZZ
-        Q = QuaternionAlgebra(QQ,-1,-1)
-        I,J,K = Q.gens()
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
+        return cls(n, **kwargs)
 
 
-        # This is like the symmetric case, but we need to be careful:
-        #
-        #   * We want conjugate-symmetry, not just symmetry.
-        #   * The diagonal will (as a result) be real.
-        #
-        S = []
-        Eij = matrix.zero(Q,n)
-        for i in range(n):
-            for j in range(i+1):
-                # "build" E_ij
-                Eij[i,j] = 1
-                if i == j:
-                    Sij = cls.real_embed(Eij)
-                    S.append(Sij)
-                else:
-                    # The second, third, and fourth ones have a minus
-                    # because they're conjugated.
-                    # Eij = Eij + Eij.transpose()
-                    Eij[j,i] = 1
-                    Sij_real = cls.real_embed(Eij)
-                    S.append(Sij_real)
-                    # Eij = I*(Eij - Eij.transpose())
-                    Eij[i,j] = I
-                    Eij[j,i] = -I
-                    Sij_I = cls.real_embed(Eij)
-                    S.append(Sij_I)
-                    # Eij = J*(Eij - Eij.transpose())
-                    Eij[i,j] = J
-                    Eij[j,i] = -J
-                    Sij_J = cls.real_embed(Eij)
-                    S.append(Sij_J)
-                    # Eij = K*(Eij - Eij.transpose())
-                    Eij[i,j] = K
-                    Eij[j,i] = -K
-                    Sij_K = cls.real_embed(Eij)
-                    S.append(Sij_K)
-                    Eij[j,i] = 0
-                # "erase" E_ij
-                Eij[i,j] = 0
-
-        # Since we embedded these, we can drop back to the "field" that we
-        # started with instead of the quaternion algebra "Q".
-        return tuple( s.change_ring(field) for s in S )
-
-
-    def __init__(self, n, **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
+class OctonionHermitianEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA):
+    r"""
+    SETUP::
 
 
-        super(QuaternionHermitianEJA, self).__init__(self._denormalized_basis(n),
-                                                     self.jordan_product,
-                                                     self.trace_inner_product,
-                                                     **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)
-        self.one.set_cache(self(idV))
+        sage: from mjo.eja.eja_algebra import (EJA,
+        ....:                                  OctonionHermitianEJA)
+        sage: from mjo.hurwitz import Octonions, OctonionMatrixAlgebra
 
 
+    EXAMPLES:
 
 
+    The 3-by-3 algebra satisfies the axioms of an EJA::
+
+        sage: OctonionHermitianEJA(3,                    # long time
+        ....:                      field=QQ,             # long time
+        ....:                      orthonormalize=False, # long time
+        ....:                      check_axioms=True)    # long time
+        Euclidean Jordan algebra of dimension 27 over Rational Field
+
+    After a change-of-basis, the 2-by-2 algebra has the same
+    multiplication table as the ten-dimensional Jordan spin algebra::
+
+        sage: A = OctonionMatrixAlgebra(2,Octonions(QQ),QQ)
+        sage: b = OctonionHermitianEJA._denormalized_basis(A)
+        sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],)
+        sage: jp = OctonionHermitianEJA.jordan_product
+        sage: ip = OctonionHermitianEJA.trace_inner_product
+        sage: J = EJA(basis,
+        ....:                          jp,
+        ....:                          ip,
+        ....:                          field=QQ,
+        ....:                          orthonormalize=False)
+        sage: J.multiplication_table()
+        +----++----+----+----+----+----+----+----+----+----+----+
+        | *  || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
+        +====++====+====+====+====+====+====+====+====+====+====+
+        | b0 || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
+        +----++----+----+----+----+----+----+----+----+----+----+
+        | b1 || b1 | b0 | 0  | 0  | 0  | 0  | 0  | 0  | 0  | 0  |
+        +----++----+----+----+----+----+----+----+----+----+----+
+        | b2 || b2 | 0  | b0 | 0  | 0  | 0  | 0  | 0  | 0  | 0  |
+        +----++----+----+----+----+----+----+----+----+----+----+
+        | b3 || b3 | 0  | 0  | b0 | 0  | 0  | 0  | 0  | 0  | 0  |
+        +----++----+----+----+----+----+----+----+----+----+----+
+        | b4 || b4 | 0  | 0  | 0  | b0 | 0  | 0  | 0  | 0  | 0  |
+        +----++----+----+----+----+----+----+----+----+----+----+
+        | b5 || b5 | 0  | 0  | 0  | 0  | b0 | 0  | 0  | 0  | 0  |
+        +----++----+----+----+----+----+----+----+----+----+----+
+        | b6 || b6 | 0  | 0  | 0  | 0  | 0  | b0 | 0  | 0  | 0  |
+        +----++----+----+----+----+----+----+----+----+----+----+
+        | b7 || b7 | 0  | 0  | 0  | 0  | 0  | 0  | b0 | 0  | 0  |
+        +----++----+----+----+----+----+----+----+----+----+----+
+        | b8 || b8 | 0  | 0  | 0  | 0  | 0  | 0  | 0  | b0 | 0  |
+        +----++----+----+----+----+----+----+----+----+----+----+
+        | b9 || b9 | 0  | 0  | 0  | 0  | 0  | 0  | 0  | 0  | b0 |
+        +----++----+----+----+----+----+----+----+----+----+----+
+
+    TESTS:
+
+    We can actually construct the 27-dimensional Albert algebra,
+    and we get the right unit element if we recompute it::
+
+        sage: J = OctonionHermitianEJA(3,                    # long time
+        ....:                          field=QQ,             # long time
+        ....:                          orthonormalize=False) # long time
+        sage: J.one.clear_cache()                            # long time
+        sage: J.one()                                        # long time
+        b0 + b9 + b26
+        sage: J.one().to_matrix()                            # long time
+        +----+----+----+
+        | e0 | 0  | 0  |
+        +----+----+----+
+        | 0  | e0 | 0  |
+        +----+----+----+
+        | 0  | 0  | e0 |
+        +----+----+----+
+
+    The 2-by-2 algebra is isomorphic to the ten-dimensional Jordan
+    spin algebra, but just to be sure, we recompute its rank::
+
+        sage: J = OctonionHermitianEJA(2,                    # long time
+        ....:                          field=QQ,             # long time
+        ....:                          orthonormalize=False) # long time
+        sage: J.rank.clear_cache()                           # long time
+        sage: J.rank()                                       # long time
+        2
+
+    """
     @staticmethod
     @staticmethod
-    def _max_random_instance_size():
+    def _max_random_instance_size(max_dimension):
         r"""
         r"""
-        The maximum rank of a random QuaternionHermitianEJA.
-        """
-        return 2 # Dimension 6
+        The maximum rank of a random OctonionHermitianEJA.
+        """
+        # There's certainly a formula for this, but with only four
+        # cases to worry about, I'm not that motivated to derive it.
+        if max_dimension >= 27:
+            return 3
+        elif max_dimension >= 10:
+            return 2
+        elif max_dimension >= 1:
+            return 1
+        else:
+            return 0
 
     @classmethod
 
     @classmethod
-    def random_instance(cls, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this type of algebra.
         """
         """
         Return a random instance of this type of algebra.
         """
-        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
         return cls(n, **kwargs)
 
+    def __init__(self, n, field=AA, **kwargs):
+        if n > 3:
+            # Otherwise we don't get an EJA.
+            raise ValueError("n cannot exceed 3")
 
 
-class HadamardEJA(ConcreteEJA):
+        # 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
+
+        from mjo.hurwitz import OctonionMatrixAlgebra
+        A = OctonionMatrixAlgebra(n, scalars=field)
+        super().__init__(A, **kwargs)
+
+        from mjo.eja.eja_cache import octonion_hermitian_eja_coeffs
+        a = octonion_hermitian_eja_coeffs(self)
+        if a is not None:
+            self.rational_algebra()._charpoly_coefficients.set_cache(a)
+
+
+class AlbertEJA(OctonionHermitianEJA):
+    r"""
+    The Albert algebra is the algebra of three-by-three Hermitian
+    matrices whose entries are octonions.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import AlbertEJA
+
+    EXAMPLES::
+
+        sage: AlbertEJA(field=QQ, orthonormalize=False)
+        Euclidean Jordan algebra of dimension 27 over Rational Field
+        sage: AlbertEJA() # long time
+        Euclidean Jordan algebra of dimension 27 over Algebraic Real Field
+
+    """
+    def __init__(self, *args, **kwargs):
+        super().__init__(3, *args, **kwargs)
+
+
+class HadamardEJA(RationalBasisEJA, ConcreteEJA):
     """
     """
-    Return the Euclidean Jordan Algebra corresponding to the set
-    `R^n` under the Hadamard product.
+    Return the Euclidean Jordan algebra on `R^n` with the Hadamard
+    (pointwise real-number multiplication) Jordan product and the
+    usual inner-product.
 
 
-    Note: this is nothing more than the Cartesian product of ``n``
-    copies of the spin algebra. Once Cartesian product algebras
-    are implemented, this can go.
+    This is nothing more than the Cartesian product of ``n`` copies of
+    the one-dimensional Jordan spin algebra, and is the most common
+    example of a non-simple Euclidean Jordan algebra.
 
     SETUP::
 
 
     SETUP::
 
@@ -2498,19 +2512,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:
 
@@ -2518,16 +2532,16 @@ class HadamardEJA(ConcreteEJA):
 
         sage: HadamardEJA(3, prefix='r').gens()
         (r0, r1, r2)
 
         sage: HadamardEJA(3, prefix='r').gens()
         (r0, r1, r2)
-
     """
     """
-    def __init__(self, n, **kwargs):
+    def __init__(self, n, field=AA, **kwargs):
+        MS = MatrixSpace(field, n, 1)
+
         if n == 0:
             jordan_product = lambda x,y: x
             inner_product = lambda x,y: x
         else:
             def jordan_product(x,y):
         if n == 0:
             jordan_product = lambda x,y: x
             inner_product = lambda x,y: x
         else:
             def jordan_product(x,y):
-                P = x.parent()
-                return P( xi*yi for (xi,yi) in zip(x,y) )
+                return MS( xi*yi for (xi,yi) in zip(x,y) )
 
             def inner_product(x,y):
                 return (x.T*y)[0,0]
 
             def inner_product(x,y):
                 return (x.T*y)[0,0]
@@ -2541,36 +2555,47 @@ 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() )
+        column_basis = tuple( MS(b) for b in FreeModule(field, n).basis() )
         super().__init__(column_basis,
                          jordan_product,
                          inner_product,
         super().__init__(column_basis,
                          jordan_product,
                          inner_product,
+                         field=field,
+                         matrix_space=MS,
                          associative=True,
                          **kwargs)
         self.rank.set_cache(n)
 
                          associative=True,
                          **kwargs)
         self.rank.set_cache(n)
 
-        if n == 0:
-            self.one.set_cache( self.zero() )
-        else:
-            self.one.set_cache( sum(self.gens()) )
+        self.one.set_cache( self.sum(self.gens()) )
 
     @staticmethod
 
     @staticmethod
-    def _max_random_instance_size():
+    def _max_random_instance_dimension():
         r"""
         r"""
-        The maximum dimension of a random HadamardEJA.
+        There's no reason to go higher than five here. That's
+        enough to get the point across.
         """
         return 5
 
         """
         return 5
 
+    @staticmethod
+    def _max_random_instance_size(max_dimension):
+        r"""
+        The maximum size (=dimension) of a random HadamardEJA.
+        """
+        return max_dimension
+
     @classmethod
     @classmethod
-    def random_instance(cls, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this type of algebra.
         """
         """
         Return a random instance of this type of algebra.
         """
-        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
 
         return cls(n, **kwargs)
 
 
-class BilinearFormEJA(ConcreteEJA):
+class BilinearFormEJA(RationalBasisEJA, ConcreteEJA):
     r"""
     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
     with the half-trace inner product and jordan product ``x*y =
     r"""
     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
     with the half-trace inner product and jordan product ``x*y =
@@ -2630,7 +2655,6 @@ class BilinearFormEJA(ConcreteEJA):
     matrix.  We opt not to orthonormalize the basis, because if we
     did, we would have to normalize the `s_{i}` in a similar manner::
 
     matrix.  We opt not to orthonormalize the basis, because if we
     did, we would have to normalize the `s_{i}` in a similar manner::
 
-        sage: set_random_seed()
         sage: n = ZZ.random_element(5)
         sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
         sage: B11 = matrix.identity(QQ,1)
         sage: n = ZZ.random_element(5)
         sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
         sage: B11 = matrix.identity(QQ,1)
@@ -2652,7 +2676,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...
@@ -2666,49 +2690,71 @@ class BilinearFormEJA(ConcreteEJA):
         # verify things, we'll skip the rest of the checks.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
         # verify things, we'll skip the rest of the checks.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
+        n = B.nrows()
+        MS = MatrixSpace(field, n, 1)
+
         def inner_product(x,y):
             return (y.T*B*x)[0,0]
 
         def jordan_product(x,y):
         def inner_product(x,y):
             return (y.T*B*x)[0,0]
 
         def jordan_product(x,y):
-            P = x.parent()
             x0 = x[0,0]
             xbar = x[1:,0]
             y0 = y[0,0]
             ybar = y[1:,0]
             z0 = inner_product(y,x)
             zbar = y0*xbar + x0*ybar
             x0 = x[0,0]
             xbar = x[1:,0]
             y0 = y[0,0]
             ybar = y[1:,0]
             z0 = inner_product(y,x)
             zbar = y0*xbar + x0*ybar
-            return P([z0] + zbar.list())
+            return MS([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( MS(b) 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,
+                         matrix_space=MS,
+                         associative=associative,
+                         **kwargs)
 
         # The rank of this algebra is two, unless we're in a
         # one-dimensional ambient space (because the rank is bounded
         # by the ambient dimension).
         self.rank.set_cache(min(n,2))
 
         # The rank of this algebra is two, unless we're in a
         # one-dimensional ambient space (because the rank is bounded
         # by the ambient dimension).
         self.rank.set_cache(min(n,2))
-
         if n == 0:
             self.one.set_cache( self.zero() )
         else:
             self.one.set_cache( self.monomial(0) )
 
     @staticmethod
         if n == 0:
             self.one.set_cache( self.zero() )
         else:
             self.one.set_cache( self.monomial(0) )
 
     @staticmethod
-    def _max_random_instance_size():
+    def _max_random_instance_dimension():
         r"""
         r"""
-        The maximum dimension of a random BilinearFormEJA.
+        There's no reason to go higher than five here. That's
+        enough to get the point across.
         """
         return 5
 
         """
         return 5
 
+    @staticmethod
+    def _max_random_instance_size(max_dimension):
+        r"""
+        The maximum size (=dimension) of a random BilinearFormEJA.
+        """
+        return max_dimension
+
     @classmethod
     @classmethod
-    def random_instance(cls, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this algebra.
         """
         """
         Return a random instance of this algebra.
         """
-        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
+
         if n.is_zero():
             B = matrix.identity(ZZ, n)
             return cls(B, **kwargs)
         if n.is_zero():
             B = matrix.identity(ZZ, n)
             return cls(B, **kwargs)
@@ -2719,6 +2765,7 @@ class BilinearFormEJA(ConcreteEJA):
         alpha = ZZ.zero()
         while alpha.is_zero():
             alpha = ZZ.random_element().abs()
         alpha = ZZ.zero()
         while alpha.is_zero():
             alpha = ZZ.random_element().abs()
+
         B22 = M.transpose()*M + alpha*I
 
         from sage.matrix.special import block_matrix
         B22 = M.transpose()*M + alpha*I
 
         from sage.matrix.special import block_matrix
@@ -2744,20 +2791,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::
@@ -2769,7 +2816,6 @@ class JordanSpinEJA(BilinearFormEJA):
 
         Ensure that we have the usual inner product on `R^n`::
 
 
         Ensure that we have the usual inner product on `R^n`::
 
-            sage: set_random_seed()
             sage: J = JordanSpinEJA.random_instance()
             sage: x,y = J.random_elements(2)
             sage: actual = x.inner_product(y)
             sage: J = JordanSpinEJA.random_instance()
             sage: x,y = J.random_elements(2)
             sage: actual = x.inner_product(y)
@@ -2778,7 +2824,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)
@@ -2789,27 +2835,24 @@ 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)
-
-    @staticmethod
-    def _max_random_instance_size():
-        r"""
-        The maximum dimension of a random JordanSpinEJA.
-        """
-        return 5
+        super().__init__(B, *args, **kwargs)
 
     @classmethod
 
     @classmethod
-    def random_instance(cls, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this type of algebra.
 
         Needed here to override the implementation for ``BilinearFormEJA``.
         """
         """
         Return a random instance of this type of algebra.
 
         Needed here to override the implementation for ``BilinearFormEJA``.
         """
-        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
 
         return cls(n, **kwargs)
 
 
-class TrivialEJA(ConcreteEJA):
+class TrivialEJA(RationalBasisEJA, ConcreteEJA):
     """
     The trivial Euclidean Jordan algebra consisting of only a zero element.
 
     """
     The trivial Euclidean Jordan algebra consisting of only a zero element.
 
@@ -2838,33 +2881,40 @@ class TrivialEJA(ConcreteEJA):
         0
 
     """
         0
 
     """
-    def __init__(self, **kwargs):
+    def __init__(self, field=AA, **kwargs):
         jordan_product = lambda x,y: x
         jordan_product = lambda x,y: x
-        inner_product = lambda x,y: 0
+        inner_product = lambda x,y: field.zero()
         basis = ()
         basis = ()
+        MS = MatrixSpace(field,0)
 
         # New defaults for keyword arguments
         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
 
         # New defaults for keyword arguments
         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,
+                         field=field,
+                         matrix_space=MS,
+                         **kwargs)
+
         # The rank is zero using my definition, namely the dimension of the
         # largest subalgebra generated by any element.
         self.rank.set_cache(0)
         self.one.set_cache( self.zero() )
 
     @classmethod
         # The rank is zero using my definition, namely the dimension of the
         # largest subalgebra generated by any element.
         self.rank.set_cache(0)
         self.one.set_cache( self.zero() )
 
     @classmethod
-    def random_instance(cls, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         # We don't take a "size" argument so the superclass method is
         # We don't take a "size" argument so the superclass method is
-        # inappropriate for us.
+        # inappropriate for us. The ``max_dimension`` argument is
+        # included so that if this method is called generically with a
+        # ``max_dimension=<whatever>`` argument, we don't try to pass
+        # it on to the algebra constructor.
         return cls(**kwargs)
 
 
         return cls(**kwargs)
 
 
-class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
-                          FiniteDimensionalEJA):
+class CartesianProductEJA(EJA):
     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
@@ -2876,6 +2926,7 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
 
         sage: from mjo.eja.eja_algebra import (random_eja,
         ....:                                  CartesianProductEJA,
 
         sage: from mjo.eja.eja_algebra import (random_eja,
         ....:                                  CartesianProductEJA,
+        ....:                                  ComplexHermitianEJA,
         ....:                                  HadamardEJA,
         ....:                                  JordanSpinEJA,
         ....:                                  RealSymmetricEJA)
         ....:                                  HadamardEJA,
         ....:                                  JordanSpinEJA,
         ....:                                  RealSymmetricEJA)
@@ -2885,7 +2936,6 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
     The Jordan product is inherited from our factors and implemented by
     our CombinatorialFreeModule Cartesian product superclass::
 
     The Jordan product is inherited from our factors and implemented by
     our CombinatorialFreeModule Cartesian product superclass::
 
-        sage: set_random_seed()
         sage: J1 = HadamardEJA(2)
         sage: J2 = RealSymmetricEJA(2)
         sage: J = cartesian_product([J1,J2])
         sage: J1 = HadamardEJA(2)
         sage: J2 = RealSymmetricEJA(2)
         sage: J = cartesian_product([J1,J2])
@@ -2960,6 +3010,55 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
         sage: CP2.is_associative()
         False
 
         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 |
+        +----++----+----+----+
+
+    The "matrix space" of a Cartesian product always consists of
+    ordered pairs (or triples, or...) whose components are the
+    matrix spaces of its factors::
+
+            sage: J1 = HadamardEJA(2)
+            sage: J2 = ComplexHermitianEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: J.matrix_space()
+            The Cartesian product of (Full MatrixSpace of 2 by 1 dense
+            matrices over Algebraic Real Field, Module of 2 by 2 matrices
+            with entries in Algebraic Field over the scalar ring Algebraic
+            Real Field)
+            sage: J.one().to_matrix()[0]
+            [1]
+            [1]
+            sage: J.one().to_matrix()[1]
+            +---+---+
+            | 1 | 0 |
+            +---+---+
+            | 0 | 1 |
+            +---+---+
+
     TESTS:
 
     All factors must share the same base field::
     TESTS:
 
     All factors must share the same base field::
@@ -2973,7 +3072,6 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
 
     The cached unit element is the same one that would be computed::
 
 
     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: J1 = random_eja()              # long time
         sage: J2 = random_eja()              # long time
         sage: J = cartesian_product([J1,J2]) # long time
@@ -2982,89 +3080,155 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
         sage: expected = J.one()             # long time
         sage: actual == expected             # long time
         True
         sage: expected = J.one()             # long time
         sage: actual == expected             # long time
         True
-
     """
     """
-    Element = FiniteDimensionalEJAElement
+    Element = CartesianProductEJAElement
+    def __init__(self, factors, **kwargs):
+        m = len(factors)
+        if m == 0:
+            return TrivialEJA()
 
 
+        self._sets = factors
 
 
-    def __init__(self, algebras, **kwargs):
-        CombinatorialFreeModule_CartesianProduct.__init__(self,
-                                                          algebras,
-                                                          **kwargs)
-        field = algebras[0].base_ring()
-        if not all( J.base_ring() == field for J in algebras ):
+        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")
 
-        associative = all( m.is_associative() for m in algebras )
-
-        # The definition of matrix_space() and self.basis() relies
-        # only on the stuff in the CFM_CartesianProduct class, which
-        # we've already initialized.
-        Js = self.cartesian_factors()
-        m = len(Js)
-        MS = self.matrix_space()
-        basis = tuple(
-            MS(tuple( self.cartesian_projection(i)(b).to_matrix()
-                      for i in range(m) ))
-            for b in self.basis()
-        )
+        # Figure out the category to use.
+        associative = all( f.is_associative() for f in factors )
+        category = EuclideanJordanAlgebras(field)
+        if associative: category = category.Associative()
+        category = category.join([category, category.CartesianProducts()])
+
+        # Compute my matrix space.  We don't simply use the
+        # ``cartesian_product()`` functor here because it acts
+        # differently on SageMath MatrixSpaces and our custom
+        # MatrixAlgebras, which are CombinatorialFreeModules. We
+        # always want the result to be represented (and indexed) as an
+        # ordered tuple. This category isn't perfect, but is good
+        # enough for what we need to do.
+        MS_cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
+        MS_cat = MS_cat.Unital().CartesianProducts()
+        MS_factors = tuple( J.matrix_space() for J in factors )
+        from sage.sets.cartesian_product import CartesianProduct
+        self._matrix_space = CartesianProduct(MS_factors, MS_cat)
+
+        self._matrix_basis = []
+        zero = self._matrix_space.zero()
+        for i in range(m):
+            for b in factors[i].matrix_basis():
+                z = list(zero)
+                z[i] = b
+                self._matrix_basis.append(z)
+
+        self._matrix_basis = tuple( self._matrix_space(b)
+                                    for b in self._matrix_basis )
+        n = len(self._matrix_basis)
+
+        # We already have what we need for the super-superclass constructor.
+        CombinatorialFreeModule.__init__(self,
+                                         field,
+                                         range(n),
+                                         prefix="b",
+                                         category=category,
+                                         bracket=False)
 
 
-        # Define jordan/inner products that operate on that matrix_basis.
-        def jordan_product(x,y):
-            return MS(tuple(
-                (Js[i](x[i])*Js[i](y[i])).to_matrix() for i in range(m)
-            ))
+        # Now create the vector space for the algebra, which will have
+        # its own set of non-ambient coordinates (in terms of the
+        # supplied basis).
+        degree = sum( f._matrix_span.ambient_vector_space().degree()
+                      for f in factors )
+        V = VectorSpace(field, degree)
+        vector_basis = tuple( V(_all2list(b)) for b in self._matrix_basis )
+
+        # Save the span of our matrix basis (when written out as long
+        # vectors) because otherwise we'll have to reconstruct it
+        # every time we want to coerce a matrix into the algebra.
+        self._matrix_span = V.span_of_basis( vector_basis, check=False)
+
+        # Since we don't (re)orthonormalize the basis, the FDEJA
+        # constructor is going to set self._deortho_matrix to the
+        # identity matrix. Here we set it to the correct value using
+        # the deortho matrices from our factors.
+        self._deortho_matrix = matrix.block_diagonal(
+            [J._deortho_matrix for J in factors]
+        )
 
 
-        def inner_product(x, y):
-            return sum(
-                Js[i](x[i]).inner_product(Js[i](y[i])) for i in range(m)
-            )
+        self._inner_product_matrix = matrix.block_diagonal(
+            [J._inner_product_matrix for J in factors]
+        )
+        self._inner_product_matrix._cache = {'hermitian': True}
+        self._inner_product_matrix.set_immutable()
 
 
-        # There's no need to check the field since it already came
-        # from an EJA. Likewise the axioms are guaranteed to be
-        # satisfied, unless the guy writing this class sucks.
-        #
-        # If you want the basis to be orthonormalized, orthonormalize
-        # the factors.
-        FiniteDimensionalEJA.__init__(self,
-                                      basis,
-                                      jordan_product,
-                                      inner_product,
-                                      field=field,
-                                      orthonormalize=False,
-                                      associative=associative,
-                                      cartesian_product=True,
-                                      check_field=False,
-                                      check_axioms=False)
-
-        ones = tuple(J.one() for J in algebras)
-        self.one.set_cache(self._cartesian_product_of_elements(ones))
-        self.rank.set_cache(sum(J.rank() for J in algebras))
+        # Building the multiplication table is a bit more tricky
+        # because we have to embed the entries of the factors'
+        # multiplication tables into the product EJA.
+        zed = self.zero()
+        self._multiplication_table = [ [zed for j in range(i+1)]
+                                       for i in range(n) ]
 
 
-    def matrix_space(self):
+        # Keep track of an offset that tallies the dimensions of all
+        # previous factors. If the second factor is dim=2 and if the
+        # first one is dim=3, then we want to skip the first 3x3 block
+        # when copying the multiplication table for the second factor.
+        offset = 0
+        for f in range(m):
+            phi_f = self.cartesian_embedding(f)
+            factor_dim = factors[f].dimension()
+            for i in range(factor_dim):
+                for j in range(i+1):
+                    f_ij = factors[f]._multiplication_table[i][j]
+                    e = phi_f(f_ij)
+                    self._multiplication_table[offset+i][offset+j] = e
+            offset += factor_dim
+
+        self.rank.set_cache(sum(J.rank() for J in factors))
+        ones = tuple(J.one().to_matrix() for J in factors)
+        self.one.set_cache(self(ones))
+
+    def _sets_keys(self):
         r"""
         r"""
-        Return the space that our matrix basis lives in as a Cartesian
-        product.
 
         SETUP::
 
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
             ....:                                  RealSymmetricEJA)
 
             ....:                                  RealSymmetricEJA)
 
-        EXAMPLES::
+        TESTS:
 
 
-            sage: J1 = HadamardEJA(1)
-            sage: J2 = RealSymmetricEJA(2)
+        The superclass uses ``_sets_keys()`` to implement its
+        ``cartesian_factors()`` method::
+
+            sage: J1 = RealSymmetricEJA(2,
+            ....:                       field=QQ,
+            ....:                       orthonormalize=False,
+            ....:                       prefix="a")
+            sage: J2 = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
             sage: J = cartesian_product([J1,J2])
             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)
+            sage: x = sum(i*J.gens()[i] for i in range(len(J.gens())))
+            sage: x.cartesian_factors()
+            (a1 + 2*a2, 3*b0 + 4*b1 + 5*b2 + 6*b3)
 
         """
 
         """
+        # Copy/pasted from CombinatorialFreeModule_CartesianProduct,
+        # but returning a tuple instead of a list.
+        return tuple(range(len(self.cartesian_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
         from sage.categories.cartesian_product import cartesian_product
-        return cartesian_product( [J.matrix_space()
-                                   for J in self.cartesian_factors()] )
+        return cartesian_product.symbol.join("%s" % factor
+                                             for factor in self._sets)
+
 
     @cached_method
     def cartesian_projection(self, i):
 
     @cached_method
     def cartesian_projection(self, i):
@@ -3126,7 +3290,6 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
 
         The answer never changes::
 
 
         The answer never changes::
 
-            sage: set_random_seed()
             sage: J1 = random_eja()
             sage: J2 = random_eja()
             sage: J = cartesian_product([J1,J2])
             sage: J1 = random_eja()
             sage: J2 = random_eja()
             sage: J = cartesian_product([J1,J2])
@@ -3136,10 +3299,13 @@ 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)
-        return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix())
+        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 EJAOperator(self,Ji,Pi.matrix())
 
     @cached_method
     def cartesian_embedding(self, i):
 
     @cached_method
     def cartesian_embedding(self, i):
@@ -3213,7 +3379,6 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
 
         The answer never changes::
 
 
         The answer never changes::
 
-            sage: set_random_seed()
             sage: J1 = random_eja()
             sage: J2 = random_eja()
             sage: J = cartesian_product([J1,J2])
             sage: J1 = random_eja()
             sage: J2 = random_eja()
             sage: J = cartesian_product([J1,J2])
@@ -3226,7 +3391,6 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
         produce the identity map, and mismatching them should produce
         the zero map::
 
         produce the identity map, and mismatching them should produce
         the zero map::
 
-            sage: set_random_seed()
             sage: J1 = random_eja()
             sage: J2 = random_eja()
             sage: J = cartesian_product([J1,J2])
             sage: J1 = random_eja()
             sage: J2 = random_eja()
             sage: J = cartesian_product([J1,J2])
@@ -3244,14 +3408,44 @@ 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)
-        return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix())
+        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 EJAOperator(Ji,self,Ei.matrix())
+
+
+    def subalgebra(self, basis, **kwargs):
+        r"""
+        Create a subalgebra of this algebra from the given basis.
+
+        Only overridden to allow us to use a special Cartesian product
+        subalgebra class.
 
 
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            ....:                                  QuaternionHermitianEJA)
+
+        EXAMPLES:
 
 
+        Subalgebras of Cartesian product EJAs have a different class
+        than those of non-Cartesian-product EJAs::
+
+            sage: J1 = HadamardEJA(2,field=QQ,orthonormalize=False)
+            sage: J2 = QuaternionHermitianEJA(0,field=QQ,orthonormalize=False)
+            sage: J = cartesian_product([J1,J2])
+            sage: K1 = J1.subalgebra((J1.one(),), orthonormalize=False)
+            sage: K = J.subalgebra((J.one(),), orthonormalize=False)
+            sage: K1.__class__ is K.__class__
+            False
+
+        """
+        from mjo.eja.eja_subalgebra import CartesianProductEJASubalgebra
+        return CartesianProductEJASubalgebra(self, basis, **kwargs)
 
 
-FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA
+EJA.CartesianProduct = CartesianProductEJA
 
 class RationalBasisCartesianProductEJA(CartesianProductEJA,
                                        RationalBasisEJA):
 
 class RationalBasisCartesianProductEJA(CartesianProductEJA,
                                        RationalBasisEJA):
@@ -3261,7 +3455,9 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
 
     SETUP::
 
 
     SETUP::
 
-        sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+        sage: from mjo.eja.eja_algebra import (EJA,
+        ....:                                  HadamardEJA,
+        ....:                                  JordanSpinEJA,
         ....:                                  RealSymmetricEJA)
 
     EXAMPLES:
         ....:                                  RealSymmetricEJA)
 
     EXAMPLES:
@@ -3278,17 +3474,308 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
         sage: J.rank()
         5
 
         sage: J.rank()
         5
 
+    TESTS:
+
+    The ``cartesian_product()`` function only uses the first factor to
+    decide where the result will live; thus we have to be careful to
+    check that all factors do indeed have a ``rational_algebra()`` method
+    before we construct an algebra that claims to have a rational basis::
+
+        sage: J1 = HadamardEJA(2)
+        sage: jp = lambda X,Y: X*Y
+        sage: ip = lambda X,Y: X[0,0]*Y[0,0]
+        sage: b1 = matrix(QQ, [[1]])
+        sage: J2 = EJA((b1,), jp, ip)
+        sage: cartesian_product([J2,J1]) # factor one not RationalBasisEJA
+        Euclidean Jordan algebra of dimension 1 over Algebraic Real
+        Field (+) Euclidean Jordan algebra of dimension 2 over Algebraic
+        Real Field
+        sage: cartesian_product([J1,J2]) # factor one is RationalBasisEJA
+        Traceback (most recent call last):
+        ...
+        ValueError: factor not a RationalBasisEJA
+
     """
     def __init__(self, algebras, **kwargs):
     """
     def __init__(self, algebras, **kwargs):
+        if not all( hasattr(r, "rational_algebra") for r in algebras ):
+            raise ValueError("factor not a RationalBasisEJA")
+
         CartesianProductEJA.__init__(self, algebras, **kwargs)
 
         CartesianProductEJA.__init__(self, algebras, **kwargs)
 
-        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
-            ])
+    @cached_method
+    def rational_algebra(self):
+        if self.base_ring() is QQ:
+            return self
+
+        return cartesian_product([
+            r.rational_algebra() for r in self.cartesian_factors()
+        ])
 
 
 RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
 
 
 
 RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
 
-random_eja = ConcreteEJA.random_instance
+def random_eja(max_dimension=None, *args, **kwargs):
+    r"""
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import random_eja
+
+    TESTS::
+
+        sage: n = ZZ.random_element(1,5)
+        sage: J = random_eja(max_dimension=n, field=QQ, orthonormalize=False)
+        sage: J.dimension() <= n
+        True
+
+    """
+    # Use the ConcreteEJA default as the total upper bound (regardless
+    # of any whether or not any individual factors set a lower limit).
+    if max_dimension is None:
+        max_dimension = ConcreteEJA._max_random_instance_dimension()
+    J1 = ConcreteEJA.random_instance(max_dimension, *args, **kwargs)
+
+
+    # Roll the dice to see if we attempt a Cartesian product.
+    dice_roll = ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1)
+    new_max_dimension = max_dimension - J1.dimension()
+    if new_max_dimension == 0 or dice_roll != 0:
+        # If it's already as big as we're willing to tolerate, just
+        # return it and don't worry about Cartesian products.
+        return J1
+    else:
+        # Use random_eja() again so we can get more than two factors
+        # if the sub-call also Decides on a cartesian product.
+        J2 = random_eja(new_max_dimension, *args, **kwargs)
+        return cartesian_product([J1,J2])
+
+
+class ComplexSkewSymmetricEJA(RationalBasisEJA, ConcreteEJA):
+    r"""
+    The skew-symmetric EJA of order `2n` described in Faraut and
+    Koranyi's Exercise III.1.b. It has dimension `2n^2 - n`.
+
+    It is (not obviously) isomorphic to the QuaternionHermitianEJA of
+    order `n`, as can be inferred by comparing rank/dimension or
+    explicitly from their "characteristic polynomial of" functions,
+    which just so happen to align nicely.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import (ComplexSkewSymmetricEJA,
+        ....:                                  QuaternionHermitianEJA)
+        sage: from mjo.eja.eja_operator import EJAOperator
+
+    EXAMPLES:
+
+    This EJA is isomorphic to the quaternions::
+
+        sage: J = ComplexSkewSymmetricEJA(2, field=QQ, orthonormalize=False)
+        sage: K = QuaternionHermitianEJA(2, field=QQ, orthonormalize=False)
+        sage: jordan_isom_matrix = matrix.diagonal(QQ,[-1,1,1,1,1,-1])
+        sage: phi = EJAOperator(J,K,jordan_isom_matrix)
+        sage: all( phi(x*y) == phi(x)*phi(y)
+        ....:      for x in J.gens()
+        ....:      for y in J.gens() )
+        True
+        sage: x,y = J.random_elements(2)
+        sage: phi(x*y) == phi(x)*phi(y)
+        True
+
+    TESTS:
+
+    Random elements should satisfy the same conditions that the basis
+    elements do::
+
+        sage: K = ComplexSkewSymmetricEJA.random_instance(field=QQ,
+        ....:                                             orthonormalize=False)
+        sage: x,y = K.random_elements(2)
+        sage: z = x*y
+        sage: x = x.to_matrix()
+        sage: y = y.to_matrix()
+        sage: z = z.to_matrix()
+        sage: all( e.is_skew_symmetric() for e in (x,y,z) )
+        True
+        sage: J = -K.one().to_matrix()
+        sage: all( e*J == J*e.conjugate() for e in (x,y,z) )
+        True
+
+    The power law in Faraut & Koranyi's II.7.a is satisfied.
+    We're in a subalgebra of theirs, but powers are still
+    defined the same::
+
+        sage: K = ComplexSkewSymmetricEJA.random_instance(field=QQ,
+        ....:                                             orthonormalize=False)
+        sage: x = K.random_element()
+        sage: k = ZZ.random_element(5)
+        sage: actual = x^k
+        sage: J = -K.one().to_matrix()
+        sage: expected = K(-J*(J*x.to_matrix())^k)
+        sage: actual == expected
+        True
+
+    """
+    @staticmethod
+    def _max_random_instance_size(max_dimension):
+        # Obtained by solving d = 2n^2 - n, which comes from noticing
+        # that, in 2x2 block form, any element of this algebra has a
+        # free skew-symmetric top-left block, a Hermitian top-right
+        # block, and two bottom blocks that are determined by the top.
+        # The ZZ-int-ZZ thing is just "floor."
+        return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/4 + 1/4))
+
+    @classmethod
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
+        """
+        Return a random instance of this type of algebra.
+        """
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
+        return cls(n, **kwargs)
+
+    @staticmethod
+    def _denormalized_basis(A):
+        """
+        SETUP::
+
+            sage: from mjo.hurwitz import ComplexMatrixAlgebra
+            sage: from mjo.eja.eja_algebra import ComplexSkewSymmetricEJA
+
+        TESTS:
+
+        The basis elements are all skew-Hermitian::
+
+            sage: d_max = ComplexSkewSymmetricEJA._max_random_instance_dimension()
+            sage: n_max = ComplexSkewSymmetricEJA._max_random_instance_size(d_max)
+            sage: n = ZZ.random_element(n_max + 1)
+            sage: A = ComplexMatrixAlgebra(2*n, scalars=QQ)
+            sage: B = ComplexSkewSymmetricEJA._denormalized_basis(A)
+            sage: all( M.is_skew_symmetric() for M in  B)
+            True
+
+        The basis elements ``b`` all satisfy ``b*J == J*b.conjugate()``,
+        as in the definition of the algebra::
+
+            sage: d_max = ComplexSkewSymmetricEJA._max_random_instance_dimension()
+            sage: n_max = ComplexSkewSymmetricEJA._max_random_instance_size(d_max)
+            sage: n = ZZ.random_element(n_max + 1)
+            sage: A = ComplexMatrixAlgebra(2*n, scalars=QQ)
+            sage: I_n = matrix.identity(ZZ, n)
+            sage: J = matrix.block(ZZ, 2, 2, (0, I_n, -I_n, 0), subdivide=False)
+            sage: J = A.from_list(J.rows())
+            sage: B = ComplexSkewSymmetricEJA._denormalized_basis(A)
+            sage: all( b*J == J*b.conjugate()  for b in B )
+            True
+
+        """
+        es = A.entry_algebra_gens()
+        gen = lambda A,m: A.monomial(m)
+
+        basis = []
+
+        # The size of the blocks. We're going to treat these thing as
+        # 2x2 block matrices,
+        #
+        #   [  x1        x2      ]
+        #   [ -x2-conj   x1-conj ]
+        #
+        # where x1 is skew-symmetric and x2 is Hermitian.
+        #
+        m = A.nrows()/2
+
+        # We only loop through the top half of the matrix, because the
+        # bottom can be constructed from the top.
+        for i in range(m):
+            # First do the top-left block, which is skew-symmetric.
+            # We can compute the bottom-right block in the process.
+            for j in range(i+1):
+                if i != j:
+                    # Skew-symmetry implies zeros for (i == j).
+                    for e in es:
+                        # Top-left block's entry.
+                        E_ij  = gen(A, (i,j,e))
+                        E_ij -= gen(A, (j,i,e))
+
+                        # Bottom-right block's entry.
+                        F_ij  = gen(A, (i+m,j+m,e)).conjugate()
+                        F_ij -= gen(A, (j+m,i+m,e)).conjugate()
+
+                        basis.append(E_ij + F_ij)
+
+            # Now do the top-right block, which is Hermitian, and compute
+            # the bottom-left block along the way.
+            for j in range(m,i+m+1):
+                if (i+m) == j:
+                    # Hermitian matrices have real diagonal entries.
+                    # Top-right block's entry.
+                    E_ii = gen(A, (i,j,es[0]))
+
+                    # Bottom-left block's entry. Don't conjugate
+                    # 'cause it's real.
+                    E_ii -= gen(A, (i+m,j-m,es[0]))
+                    basis.append(E_ii)
+                else:
+                    for e in es:
+                        # Top-right block's entry. BEWARE! We're not
+                        # reflecting across the main diagonal as in
+                        # (i,j)~(j,i). We're only reflecting across
+                        # the diagonal for the top-right block.
+                        E_ij  = gen(A, (i,j,e))
+
+                        # Shift it back to non-offset coords, transpose,
+                        # conjugate, and put it back:
+                        #
+                        # (i,j) -> (i,j-m) -> (j-m, i) -> (j-m, i+m)
+                        E_ij += gen(A, (j-m,i+m,e)).conjugate()
+
+                        # Bottom-left's block's below-diagonal entry.
+                        # Just shift the top-right coords down m and
+                        # left m.
+                        F_ij  = -gen(A, (i+m,j-m,e)).conjugate()
+                        F_ij += -gen(A, (j,i,e)) # double-conjugate cancels
+
+                        basis.append(E_ij + F_ij)
+
+        return tuple( basis )
+
+    @staticmethod
+    @cached_method
+    def _J_matrix(matrix_space):
+        n = matrix_space.nrows() // 2
+        F = matrix_space.base_ring()
+        I_n = matrix.identity(F, n)
+        J = matrix.block(F, 2, 2, (0, I_n, -I_n, 0), subdivide=False)
+        return matrix_space.from_list(J.rows())
+
+    def J_matrix(self):
+        return ComplexSkewSymmetricEJA._J_matrix(self.matrix_space())
+
+    def __init__(self, n, field=AA, **kwargs):
+        # New code; always check the axioms.
+        #if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
+
+        from mjo.hurwitz import ComplexMatrixAlgebra
+        A = ComplexMatrixAlgebra(2*n, scalars=field)
+        J = ComplexSkewSymmetricEJA._J_matrix(A)
+
+        def jordan_product(X,Y):
+            return (X*J*Y + Y*J*X)/2
+
+        def inner_product(X,Y):
+            return (X.conjugate_transpose()*Y).trace().real()
+
+        super().__init__(self._denormalized_basis(A),
+                         jordan_product,
+                         inner_product,
+                         field=field,
+                         matrix_space=A,
+                         **kwargs)
+
+        # This algebra is conjectured (by me) to be isomorphic to
+        # the quaternion Hermitian EJA of size n, and the rank
+        # would follow from that.
+        #self.rank.set_cache(n)
+        self.one.set_cache( self(-J) )