]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: add the AlbertEJA class.
[sage.d.git] / mjo / eja / eja_algebra.py
index 60a23c8b148f32f44cbbacbd71140c0ed5fc408f..ca3df5fea55949225a51f8255a2c13a33c37d9f2 100644 (file)
@@ -32,22 +32,25 @@ for these simple algebras:
   * :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 two other example constructions,
 
+  * :class:`JordanSpinEJA`
   * :class:`HadamardEJA`
+  * :class:`AlbertEJA`
   * :class:`TrivialEJA`
 
 The Jordan spin algebra is a bilinear form algebra where the bilinear
 form is the identity. The Hadamard EJA is simply a Cartesian product
-of one-dimensional spin algebras. And last but not least, the trivial
-EJA is exactly what you think. Cartesian products of these are also
-supported using the usual ``cartesian_product()`` function; as a
-result, we support (up to isomorphism) all Euclidean Jordan algebras
-that don't involve octonions.
+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.
 
 SETUP::
 
@@ -59,13 +62,10 @@ EXAMPLES::
     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.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
@@ -142,19 +142,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                  cartesian_product=False,
                  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)
-        if n > 0:
-            if hasattr(basis[0], 'cartesian_factors'):
-                self._matrix_basis_is_cartesian = True
 
         if check_field:
             if not field.is_subring(RR):
@@ -163,21 +153,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                 # we've specified a real embedding.
                 raise ValueError("scalar field is not real")
 
-        # If the basis given to us wasn't over the field that it's
-        # supposed to be over, fix that. Or, you know, crash.
-        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 )
-
-
         if check_axioms:
             # Check commutativity of the Jordan and inner-products.
             # This has to be done before we build the multiplication
@@ -213,7 +188,10 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             # Element subalgebras can take advantage of this.
             category = category.Associative()
         if cartesian_product:
-            category = category.CartesianProducts()
+            # Use join() here because otherwise we only get the
+            # "Cartesian product of..." and not the things themselves.
+            category = category.join([category,
+                                      category.CartesianProducts()])
 
         # Call the superclass constructor so that we can use its from_vector()
         # method to build our multiplication table.
@@ -228,7 +206,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
-        # we see in things like x = 1*e1 + 2*e2.
+        # we see in things like x = 1*b1 + 2*b2.
         vector_basis = basis
 
         degree = 0
@@ -355,16 +333,16 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: set_random_seed()
             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)
-            ....:     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
 
         """
@@ -473,9 +451,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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()) )
 
@@ -555,9 +533,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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)*z - x*(y*z)
 
                     if diff.norm() > epsilon:
@@ -586,9 +564,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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)
 
                     if diff.abs() > epsilon:
@@ -636,7 +614,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             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
 
         TESTS:
 
@@ -911,15 +889,15 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: J = JordanSpinEJA(4)
             sage: J.multiplication_table()
             +----++----+----+----+----+
-            | *  || e0 | e1 | e2 | e3 |
+            | *  || b0 | b1 | b2 | b3 |
             +====++====+====+====+====+
-            | e0 || e0 | e1 | e2 | e3 |
+            | b0 || b0 | b1 | b2 | b3 |
             +----++----+----+----+----+
-            | e1 || e1 | e0 | 0  | 0  |
+            | b1 || b1 | b0 | 0  | 0  |
             +----++----+----+----+----+
-            | e2 || e2 | 0  | e0 | 0  |
+            | b2 || b2 | 0  | b0 | 0  |
             +----++----+----+----+----+
-            | e3 || e3 | 0  | 0  | e0 |
+            | b3 || b3 | 0  | 0  | b0 |
             +----++----+----+----+----+
 
         """
@@ -929,7 +907,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         # And to each subsequent row, prepend an entry that belongs to
         # the left-side "header column."
-        M += [ [self.gens()[i]] + [ self.gens()[i]*self.gens()[j]
+        M += [ [self.monomial(i)] + [ self.monomial(i)*self.monomial(j)
                                     for j in range(n) ]
                for i in range(n) ]
 
@@ -973,7 +951,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: J = RealSymmetricEJA(2)
             sage: J.basis()
-            Finite family {0: e0, 1: e1, 2: e2}
+            Finite family {0: b0, 1: b1, 2: b2}
             sage: J.matrix_basis()
             (
             [1 0]  [                  0 0.7071067811865475?]  [0 0]
@@ -984,7 +962,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: J = JordanSpinEJA(2)
             sage: J.basis()
-            Finite family {0: e0, 1: e1}
+            Finite family {0: b0, 1: b1}
             sage: J.matrix_basis()
             (
             [1]  [0]
@@ -1066,20 +1044,20 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: J = HadamardEJA(5)
             sage: J.one()
-            e0 + e1 + e2 + e3 + e4
+            b0 + b1 + b2 + b3 + b4
 
         The unit element in the Hadamard EJA is inherited in the
         subalgebras generated by its elements::
 
             sage: J = HadamardEJA(5)
             sage: J.one()
-            e0 + e1 + e2 + e3 + e4
+            b0 + b1 + b2 + b3 + b4
             sage: x = sum(J.gens())
             sage: A = x.subalgebra_generated_by(orthonormalize=False)
             sage: A.one()
-            f0
+            c0
             sage: A.one().superalgebra_element()
-            e0 + e1 + e2 + e3 + e4
+            b0 + b1 + b2 + b3 + b4
 
         TESTS:
 
@@ -1431,7 +1409,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         def L_x_i_j(i,j):
             # From a result in my book, these are the entries of the
             # basis representation of L_x.
-            return sum( vars[k]*self.gens()[k].operator().matrix()[i,j]
+            return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
                         for k in range(n) )
 
         L_x = matrix(F, n, n, L_x_i_j)
@@ -1566,7 +1544,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
 class RationalBasisEJA(FiniteDimensionalEJA):
     r"""
-    New class for algebras whose supplied basis elements have all rational entries.
+    Algebras whose supplied basis elements have all rational entries.
 
     SETUP::
 
@@ -1679,7 +1657,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 )
 
-class ConcreteEJA(RationalBasisEJA):
+class ConcreteEJA(FiniteDimensionalEJA):
     r"""
     A class for the Euclidean Jordan algebras that we know by name.
 
@@ -1748,6 +1726,19 @@ class ConcreteEJA(RationalBasisEJA):
 
 
 class MatrixEJA:
+    @staticmethod
+    def jordan_product(X,Y):
+        return (X*Y + Y*X)/2
+
+    @staticmethod
+    def trace_inner_product(X,Y):
+        r"""
+        A trace inner-product for matrices that aren't embedded in the
+        reals. It takes MATRICES as arguments, not EJA elements.
+        """
+        return (X*Y).trace().real()
+
+class RealEmbeddedMatrixEJA(MatrixEJA):
     @staticmethod
     def dimension_over_reals():
         r"""
@@ -1793,9 +1784,6 @@ class MatrixEJA:
             raise ValueError("the matrix 'M' must be a real embedding")
         return M
 
-    @staticmethod
-    def jordan_product(X,Y):
-        return (X*Y + Y*X)/2
 
     @classmethod
     def trace_inner_product(cls,X,Y):
@@ -1804,29 +1792,11 @@ class MatrixEJA:
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
-            ....:                                  ComplexHermitianEJA,
+            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
             ....:                                  QuaternionHermitianEJA)
 
         EXAMPLES::
 
-        This gives the same answer as it would if we computed the trace
-        from the unembedded (original) matrices::
-
-            sage: set_random_seed()
-            sage: J = RealSymmetricEJA.random_instance()
-            sage: x,y = J.random_elements(2)
-            sage: Xe = x.to_matrix()
-            sage: Ye = y.to_matrix()
-            sage: X = J.real_unembed(Xe)
-            sage: Y = J.real_unembed(Ye)
-            sage: expected = (X*Y).trace()
-            sage: actual = J.trace_inner_product(Xe,Ye)
-            sage: actual == expected
-            True
-
-        ::
-
             sage: set_random_seed()
             sage: J = ComplexHermitianEJA.random_instance()
             sage: x,y = J.random_elements(2)
@@ -1854,27 +1824,15 @@ class MatrixEJA:
             True
 
         """
-        Xu = cls.real_unembed(X)
-        Yu = cls.real_unembed(Y)
-        tr = (Xu*Yu).trace()
-
-        try:
-            # Works in QQ, AA, RDF, et cetera.
-            return tr.real()
-        except AttributeError:
-            # A quaternion doesn't have a real() method, but does
-            # have coefficient_tuple() method that returns the
-            # coefficients of 1, i, j, and k -- in that order.
-            return tr.coefficient_tuple()[0]
-
-
-class RealMatrixEJA(MatrixEJA):
-    @staticmethod
-    def dimension_over_reals():
-        return 1
-
+        # This does in fact compute the real part of the trace.
+        # If we compute the trace of e.g. a complex matrix M,
+        # then we do so by adding up its diagonal entries --
+        # call them z_1 through z_n. The real embedding of z_1
+        # will be a 2-by-2 REAL matrix [a, b; -b, a] whose trace
+        # as a REAL matrix will be 2*a = 2*Re(z_1). And so forth.
+        return (X*Y).trace()/cls.dimension_over_reals()
 
-class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
+class RealSymmetricEJA(ConcreteEJA, RationalBasisEJA, MatrixEJA):
     """
     The rank-n simple EJA consisting of real symmetric n-by-n
     matrices, the usual symmetric Jordan product, and the trace inner
@@ -1887,13 +1845,13 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
     EXAMPLES::
 
         sage: J = RealSymmetricEJA(2)
-        sage: e0, e1, e2 = J.gens()
-        sage: e0*e0
-        e0
-        sage: e1*e1
-        1/2*e0 + 1/2*e2
-        sage: e2*e2
-        e2
+        sage: b0, b1, b2 = J.gens()
+        sage: b0*b0
+        b0
+        sage: b1*b1
+        1/2*b0 + 1/2*b2
+        sage: b2*b2
+        b2
 
     In theory, our "field" can be any subfield of the reals::
 
@@ -1940,7 +1898,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
 
     """
     @classmethod
-    def _denormalized_basis(cls, n):
+    def _denormalized_basis(cls, n, field):
         """
         Return a basis for the space of real symmetric n-by-n matrices.
 
@@ -1952,7 +1910,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
-            sage: B = RealSymmetricEJA._denormalized_basis(n)
+            sage: B = RealSymmetricEJA._denormalized_basis(n,ZZ)
             sage: all( M.is_symmetric() for M in  B)
             True
 
@@ -1962,7 +1920,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
         S = []
         for i in range(n):
             for j in range(i+1):
-                Eij = matrix(ZZ, n, lambda k,l: k==i and l==j)
+                Eij = matrix(field, n, lambda k,l: k==i and l==j)
                 if i == j:
                     Sij = Eij
                 else:
@@ -1983,7 +1941,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
         return cls(n, **kwargs)
 
-    def __init__(self, n, **kwargs):
+    def __init__(self, n, field=AA, **kwargs):
         # We know this is a valid EJA, but will double-check
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
@@ -1992,9 +1950,10 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
         if n <= 1:
             associative = True
 
-        super().__init__(self._denormalized_basis(n),
+        super().__init__(self._denormalized_basis(n,field),
                          self.jordan_product,
                          self.trace_inner_product,
+                         field=field,
                          associative=associative,
                          **kwargs)
 
@@ -2002,12 +1961,12 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         self.rank.set_cache(n)
-        idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
+        idV = self.matrix_space().one()
         self.one.set_cache(self(idV))
 
 
 
-class ComplexMatrixEJA(MatrixEJA):
+class ComplexMatrixEJA(RealEmbeddedMatrixEJA):
     # A manual dictionary-cache for the complex_extension() method,
     # since apparently @classmethods can't also be @cached_methods.
     _complex_extension = {}
@@ -2155,7 +2114,7 @@ class ComplexMatrixEJA(MatrixEJA):
         return matrix(F, n/d, elements)
 
 
-class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
+class ComplexHermitianEJA(ConcreteEJA, RationalBasisEJA, ComplexMatrixEJA):
     """
     The rank-n simple EJA consisting of complex Hermitian n-by-n
     matrices over the real numbers, the usual symmetric Jordan product,
@@ -2214,7 +2173,7 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
     """
 
     @classmethod
-    def _denormalized_basis(cls, n):
+    def _denormalized_basis(cls, n, field):
         """
         Returns a basis for the space of complex Hermitian n-by-n matrices.
 
@@ -2232,15 +2191,14 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
-            sage: B = ComplexHermitianEJA._denormalized_basis(n)
+            sage: B = ComplexHermitianEJA._denormalized_basis(n,ZZ)
             sage: all( M.is_symmetric() for M in  B)
             True
 
         """
-        field = ZZ
-        R = PolynomialRing(field, 'z')
+        R = PolynomialRing(ZZ, 'z')
         z = R.gen()
-        F = field.extension(z**2 + 1, 'I')
+        F = ZZ.extension(z**2 + 1, 'I')
         I = F.gen(1)
 
         # This is like the symmetric case, but we need to be careful:
@@ -2271,12 +2229,12 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
                 # "erase" E_ij
                 Eij[i,j] = 0
 
-        # Since we embedded these, we can drop back to the "field" that we
-        # started with instead of the complex extension "F".
+        # Since we embedded the entries, we can drop back to the
+        # desired real "field" instead of the extension "F".
         return tuple( s.change_ring(field) for s in S )
 
 
-    def __init__(self, n, **kwargs):
+    def __init__(self, n, field=AA, **kwargs):
         # We know this is a valid EJA, but will double-check
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
@@ -2285,9 +2243,10 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
         if n <= 1:
             associative = True
 
-        super().__init__(self._denormalized_basis(n),
+        super().__init__(self._denormalized_basis(n,field),
                          self.jordan_product,
                          self.trace_inner_product,
+                         field=field,
                          associative=associative,
                          **kwargs)
         # TODO: this could be factored out somehow, but is left here
@@ -2309,7 +2268,7 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
         return cls(n, **kwargs)
 
-class QuaternionMatrixEJA(MatrixEJA):
+class QuaternionMatrixEJA(RealEmbeddedMatrixEJA):
 
     # A manual dictionary-cache for the quaternion_extension() method,
     # since apparently @classmethods can't also be @cached_methods.
@@ -2457,7 +2416,9 @@ class QuaternionMatrixEJA(MatrixEJA):
         return matrix(Q, n/d, elements)
 
 
-class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
+class QuaternionHermitianEJA(ConcreteEJA,
+                             RationalBasisEJA,
+                             QuaternionMatrixEJA):
     r"""
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the
@@ -2515,7 +2476,7 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
 
     """
     @classmethod
-    def _denormalized_basis(cls, n):
+    def _denormalized_basis(cls, n, field):
         """
         Returns a basis for the space of quaternion Hermitian n-by-n matrices.
 
@@ -2533,12 +2494,11 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
-            sage: B = QuaternionHermitianEJA._denormalized_basis(n)
+            sage: B = QuaternionHermitianEJA._denormalized_basis(n,ZZ)
             sage: all( M.is_symmetric() for M in B )
             True
 
         """
-        field = ZZ
         Q = QuaternionAlgebra(QQ,-1,-1)
         I,J,K = Q.gens()
 
@@ -2582,12 +2542,12 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
                 # "erase" E_ij
                 Eij[i,j] = 0
 
-        # Since we embedded these, we can drop back to the "field" that we
-        # started with instead of the quaternion algebra "Q".
+        # Since we embedded the entries, we can drop back to the
+        # desired real "field" instead of the quaternion algebra "Q".
         return tuple( s.change_ring(field) for s in S )
 
 
-    def __init__(self, n, **kwargs):
+    def __init__(self, n, field=AA, **kwargs):
         # We know this is a valid EJA, but will double-check
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
@@ -2596,9 +2556,10 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
         if n <= 1:
             associative = True
 
-        super().__init__(self._denormalized_basis(n),
+        super().__init__(self._denormalized_basis(n,field),
                          self.jordan_product,
                          self.trace_inner_product,
+                         field=field,
                          associative=associative,
                          **kwargs)
 
@@ -2625,8 +2586,215 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
         return cls(n, **kwargs)
 
+class OctonionHermitianEJA(ConcreteEJA, MatrixEJA):
+    r"""
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
+        ....:                                  OctonionHermitianEJA)
+
+    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: b = OctonionHermitianEJA._denormalized_basis(2,QQ)
+        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 = FiniteDimensionalEJA(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
 
-class HadamardEJA(ConcreteEJA):
+    """
+    @staticmethod
+    def _max_random_instance_size():
+        r"""
+        The maximum rank of a random QuaternionHermitianEJA.
+        """
+        return 1 # Dimension 1
+
+    @classmethod
+    def random_instance(cls, **kwargs):
+        """
+        Return a random instance of this type of algebra.
+        """
+        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        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")
+
+        # 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(n,field),
+                         self.jordan_product,
+                         self.trace_inner_product,
+                         field=field,
+                         **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 = self.matrix_space().one()
+        self.one.set_cache(self(idV))
+
+
+    @classmethod
+    def _denormalized_basis(cls, n, field):
+        """
+        Returns a basis for the space of octonion Hermitian n-by-n
+        matrices.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
+
+        EXAMPLES::
+
+            sage: B = OctonionHermitianEJA._denormalized_basis(3,QQ)
+            sage: all( M.is_hermitian() for M in B )
+            True
+            sage: len(B)
+            27
+
+        """
+        from mjo.octonions import OctonionMatrixAlgebra
+        MS = OctonionMatrixAlgebra(n, scalars=field)
+        es = MS.entry_algebra().gens()
+
+        basis = []
+        for i in range(n):
+            for j in range(i+1):
+                if i == j:
+                    E_ii = MS.monomial( (i,j,es[0]) )
+                    basis.append(E_ii)
+                else:
+                    for e in es:
+                        E_ij  = MS.monomial( (i,j,e)             )
+                        ec = e.conjugate()
+                        # If the conjugate has a negative sign in front
+                        # of it, (j,i,ec) won't be a monomial!
+                        if (j,i,ec) in MS.indices():
+                            E_ij += MS.monomial( (j,i,ec) )
+                        else:
+                            E_ij -= MS.monomial( (j,i,-ec) )
+                        basis.append(E_ij)
+
+        return tuple( basis )
+
+    @staticmethod
+    def trace_inner_product(X,Y):
+        r"""
+        The octonions don't know that the reals are embedded in them,
+        so we have to take the e0 component ourselves.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
+
+        TESTS::
+
+            sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2
+
+        """
+        return (X*Y).trace().real().coefficient(0)
+
+
+class AlbertEJA(OctonionHermitianEJA):
+    r"""
+    The Albert algebra is the algebra of three-by-three Hermitian
+    matrices whose entries are octonions.
+
+    SETUP::
+
+        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(ConcreteEJA, RationalBasisEJA):
     """
     Return the Euclidean Jordan Algebra corresponding to the set
     `R^n` under the Hadamard product.
@@ -2644,19 +2812,19 @@ class HadamardEJA(ConcreteEJA):
     This multiplication table can be verified by hand::
 
         sage: J = HadamardEJA(3)
-        sage: e0,e1,e2 = J.gens()
-        sage: e0*e0
-        e0
-        sage: e0*e1
+        sage: b0,b1,b2 = J.gens()
+        sage: b0*b0
+        b0
+        sage: b0*b1
         0
-        sage: e0*e2
+        sage: b0*b2
         0
-        sage: e1*e1
-        e1
-        sage: e1*e2
+        sage: b1*b1
+        b1
+        sage: b1*b2
         0
-        sage: e2*e2
-        e2
+        sage: b2*b2
+        b2
 
     TESTS:
 
@@ -2666,7 +2834,7 @@ class HadamardEJA(ConcreteEJA):
         (r0, r1, r2)
 
     """
-    def __init__(self, n, **kwargs):
+    def __init__(self, n, field=AA, **kwargs):
         if n == 0:
             jordan_product = lambda x,y: x
             inner_product = lambda x,y: x
@@ -2687,10 +2855,12 @@ class HadamardEJA(ConcreteEJA):
         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
-        column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )
+        column_basis = tuple( b.column()
+                              for b in FreeModule(field, n).basis() )
         super().__init__(column_basis,
                          jordan_product,
                          inner_product,
+                         field=field,
                          associative=True,
                          **kwargs)
         self.rank.set_cache(n)
@@ -2716,7 +2886,7 @@ class HadamardEJA(ConcreteEJA):
         return cls(n, **kwargs)
 
 
-class BilinearFormEJA(ConcreteEJA):
+class BilinearFormEJA(ConcreteEJA, RationalBasisEJA):
     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 =
@@ -2798,7 +2968,7 @@ class BilinearFormEJA(ConcreteEJA):
         True
 
     """
-    def __init__(self, B, **kwargs):
+    def __init__(self, B, field=AA, **kwargs):
         # The matrix "B" is supplied by the user in most cases,
         # so it makes sense to check whether or not its positive-
         # definite unless we are specifically asked not to...
@@ -2826,7 +2996,8 @@ class BilinearFormEJA(ConcreteEJA):
             return P([z0] + zbar.list())
 
         n = B.nrows()
-        column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )
+        column_basis = tuple( b.column()
+                              for b in FreeModule(field, n).basis() )
 
         # TODO: I haven't actually checked this, but it seems legit.
         associative = False
@@ -2836,6 +3007,7 @@ class BilinearFormEJA(ConcreteEJA):
         super().__init__(column_basis,
                          jordan_product,
                          inner_product,
+                         field=field,
                          associative=associative,
                          **kwargs)
 
@@ -2897,20 +3069,20 @@ class JordanSpinEJA(BilinearFormEJA):
     This multiplication table can be verified by hand::
 
         sage: J = JordanSpinEJA(4)
-        sage: e0,e1,e2,e3 = J.gens()
-        sage: e0*e0
-        e0
-        sage: e0*e1
-        e1
-        sage: e0*e2
-        e2
-        sage: e0*e3
-        e3
-        sage: e1*e2
+        sage: b0,b1,b2,b3 = J.gens()
+        sage: b0*b0
+        b0
+        sage: b0*b1
+        b1
+        sage: b0*b2
+        b2
+        sage: b0*b3
+        b3
+        sage: b1*b2
         0
-        sage: e1*e3
+        sage: b1*b3
         0
-        sage: e2*e3
+        sage: b2*b3
         0
 
     We can change the generator prefix::
@@ -2931,7 +3103,7 @@ class JordanSpinEJA(BilinearFormEJA):
             True
 
     """
-    def __init__(self, n, **kwargs):
+    def __init__(self, n, *args, **kwargs):
         # This is a special case of the BilinearFormEJA with the
         # identity matrix as its bilinear form.
         B = matrix.identity(ZZ, n)
@@ -2942,7 +3114,7 @@ class JordanSpinEJA(BilinearFormEJA):
 
         # But also don't pass check_field=False here, because the user
         # can pass in a field!
-        super().__init__(B, **kwargs)
+        super().__init__(B, *args, **kwargs)
 
     @staticmethod
     def _max_random_instance_size():
@@ -2962,7 +3134,7 @@ class JordanSpinEJA(BilinearFormEJA):
         return cls(n, **kwargs)
 
 
-class TrivialEJA(ConcreteEJA):
+class TrivialEJA(ConcreteEJA, RationalBasisEJA):
     """
     The trivial Euclidean Jordan algebra consisting of only a zero element.
 
@@ -3018,8 +3190,7 @@ class TrivialEJA(ConcreteEJA):
         return cls(**kwargs)
 
 
-class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
-                          FiniteDimensionalEJA):
+class CartesianProductEJA(FiniteDimensionalEJA):
     r"""
     The external (orthogonal) direct sum of two or more Euclidean
     Jordan algebras. Every Euclidean Jordan algebra decomposes into an
@@ -3115,6 +3286,33 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
         sage: CP2.is_associative()
         False
 
+    Cartesian products of Cartesian products work::
+
+        sage: J1 = JordanSpinEJA(1)
+        sage: J2 = JordanSpinEJA(1)
+        sage: J3 = JordanSpinEJA(1)
+        sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
+        sage: J.multiplication_table()
+        +----++----+----+----+
+        | *  || b0 | b1 | b2 |
+        +====++====+====+====+
+        | b0 || b0 | 0  | 0  |
+        +----++----+----+----+
+        | b1 || 0  | b1 | 0  |
+        +----++----+----+----+
+        | b2 || 0  | 0  | b2 |
+        +----++----+----+----+
+        sage: HadamardEJA(3).multiplication_table()
+        +----++----+----+----+
+        | *  || b0 | b1 | b2 |
+        +====++====+====+====+
+        | b0 || b0 | 0  | 0  |
+        +----++----+----+----+
+        | b1 || 0  | b1 | 0  |
+        +----++----+----+----+
+        | b2 || 0  | 0  | b2 |
+        +----++----+----+----+
+
     TESTS:
 
     All factors must share the same base field::
@@ -3142,37 +3340,41 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
     Element = FiniteDimensionalEJAElement
 
 
-    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 ):
+    def __init__(self, factors, **kwargs):
+        m = len(factors)
+        if m == 0:
+            return TrivialEJA()
+
+        self._sets = factors
+
+        field = factors[0].base_ring()
+        if not all( J.base_ring() == field for J in factors ):
             raise ValueError("all factors must share the same base field")
 
-        associative = all( m.is_associative() for m in algebras )
+        associative = all( f.is_associative() for f in factors )
 
-        # 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()
-        )
+        basis = []
+        zero = MS.zero()
+        for i in range(m):
+            for b in factors[i].matrix_basis():
+                z = list(zero)
+                z[i] = b
+                basis.append(z)
+
+        basis = tuple( MS(b) for b in basis )
 
         # Define jordan/inner products that operate on that matrix_basis.
         def jordan_product(x,y):
             return MS(tuple(
-                (Js[i](x[i])*Js[i](y[i])).to_matrix() for i in range(m)
+                (factors[i](x[i])*factors[i](y[i])).to_matrix()
+                for i in range(m)
             ))
 
         def inner_product(x, y):
             return sum(
-                Js[i](x[i]).inner_product(Js[i](y[i])) for i in range(m)
+                factors[i](x[i]).inner_product(factors[i](y[i]))
+                for i in range(m)
             )
 
         # There's no need to check the field since it already came
@@ -3192,96 +3394,42 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
                                       check_field=False,
                                       check_axioms=False)
 
-        ones = tuple(J.one() for J in algebras)
-        self.one.set_cache(self._cartesian_product_of_elements(ones))
-        self.rank.set_cache(sum(J.rank() for J in algebras))
-
-    def _monomial_to_generator(self, mon):
-        r"""
-        Convert a monomial index into a generator index.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import random_eja
-
-        TESTS::
-
-            sage: J1 = random_eja(field=QQ, orthonormalize=False)
-            sage: J2 = random_eja(field=QQ, orthonormalize=False)
-            sage: J = cartesian_product([J1,J2])
-            sage: all( J.monomial(m)
-            ....:      ==
-            ....:      J.gens()[J._monomial_to_generator(m)]
-            ....:      for m in J.basis().keys() )
-
-        """
-        # The superclass method indexes into a matrix, so we have to
-        # turn the tuples i and j into integers. This is easy enough
-        # given that the first coordinate of i and j corresponds to
-        # the factor, and the second coordinate corresponds to the
-        # index of the generator within that factor.
-        try:
-            factor = mon[0]
-        except TypeError: # 'int' object is not subscriptable
-            return mon
-        idx_in_factor = self._monomial_to_generator(mon[1])
+        ones = tuple(J.one().to_matrix() for J in factors)
+        self.one.set_cache(self(ones))
+        self.rank.set_cache(sum(J.rank() for J in factors))
 
-        offset = sum( f.dimension()
-                      for f in self.cartesian_factors()[:factor] )
-        return offset + idx_in_factor
+    def cartesian_factors(self):
+        # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
+        return self._sets
 
-    def product_on_basis(self, i, j):
+    def cartesian_factor(self, i):
         r"""
-        Return the product of the monomials indexed by ``i`` and ``j``.
-
-        This overrides the superclass method because here, both ``i``
-        and ``j`` will be ordered pairs.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import (HadamardEJA,
-            ....:                                  JordanSpinEJA,
-            ....:                                  QuaternionHermitianEJA,
-            ....:                                  RealSymmetricEJA,)
-
-        EXAMPLES::
-
-            sage: J1 = JordanSpinEJA(2, field=QQ)
-            sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
-            sage: J3 = HadamardEJA(1, field=QQ)
-            sage: K1 = cartesian_product([J1,J2])
-            sage: K2 = cartesian_product([K1,J3])
-            sage: list(K2.basis())
-            [e(0, (0, 0)), e(0, (0, 1)), e(0, (1, 0)), e(0, (1, 1)),
-            e(0, (1, 2)), e(1, 0)]
-            sage: sage: g = K2.gens()
-            sage: (g[0] + 2*g[3]) * (g[1] - 4*g[2])
-            e(0, (0, 1)) - 4*e(0, (1, 1))
-
-        TESTS::
-
-            sage: J1 = RealSymmetricEJA(1,field=QQ)
-            sage: J2 = QuaternionHermitianEJA(1,field=QQ)
-            sage: J = cartesian_product([J1,J2])
-            sage: x = sum(J.gens())
-            sage: x == J.one()
-            True
-            sage: x*x == x
-            True
-
+        Return the ``i``th factor of this algebra.
         """
-        l = self._monomial_to_generator(i)
-        m = self._monomial_to_generator(j)
-        return FiniteDimensionalEJA.product_on_basis(self, l, m)
+        return self._sets[i]
+
+    def _repr_(self):
+        # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
+        from sage.categories.cartesian_product import cartesian_product
+        return cartesian_product.symbol.join("%s" % factor
+                                             for factor in self._sets)
 
     def matrix_space(self):
         r"""
         Return the space that our matrix basis lives in as a Cartesian
         product.
 
+        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.
+
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
+            ....:                                  HadamardEJA,
+            ....:                                  OctonionHermitianEJA,
             ....:                                  RealSymmetricEJA)
 
         EXAMPLES::
@@ -3294,10 +3442,44 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
             matrices over Algebraic Real Field, Full MatrixSpace of 2
             by 2 dense matrices over Algebraic Real Field)
 
+        ::
+
+            sage: J1 = ComplexHermitianEJA(1)
+            sage: J2 = ComplexHermitianEJA(1)
+            sage: J = cartesian_product([J1,J2])
+            sage: J.one().to_matrix()[0]
+            [1 0]
+            [0 1]
+            sage: J.one().to_matrix()[1]
+            [1 0]
+            [0 1]
+
+        ::
+
+            sage: J1 = OctonionHermitianEJA(1)
+            sage: J2 = OctonionHermitianEJA(1)
+            sage: J = cartesian_product([J1,J2])
+            sage: J.one().to_matrix()[0]
+            +----+
+            | e0 |
+            +----+
+            sage: J.one().to_matrix()[1]
+            +----+
+            | e0 |
+            +----+
+
         """
-        from sage.categories.cartesian_product import cartesian_product
-        return cartesian_product( [J.matrix_space()
-                                   for J in self.cartesian_factors()] )
+        scalars = self.cartesian_factor(0).base_ring()
+
+        # This category isn't perfect, but is good enough for what we
+        # need to do.
+        cat = MagmaticAlgebras(scalars).FiniteDimensional().WithBasis()
+        cat = cat.Unital().CartesianProducts()
+        factors = tuple( J.matrix_space() for J in self.cartesian_factors() )
+
+        from sage.sets.cartesian_product import CartesianProduct
+        return CartesianProduct(factors, cat)
+
 
     @cached_method
     def cartesian_projection(self, i):
@@ -3369,9 +3551,12 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
             True
 
         """
-        Ji = self.cartesian_factors()[i]
-        # Requires the fix on Trac 31421/31422 to work!
-        Pi = super().cartesian_projection(i)
+        offset = sum( self.cartesian_factor(k).dimension()
+                      for k in range(i) )
+        Ji = self.cartesian_factor(i)
+        Pi = self._module_morphism(lambda j: Ji.monomial(j - offset),
+                                   codomain=Ji)
+
         return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix())
 
     @cached_method
@@ -3477,9 +3662,11 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
             True
 
         """
-        Ji = self.cartesian_factors()[i]
-        # Requires the fix on Trac 31421/31422 to work!
-        Ei = super().cartesian_embedding(i)
+        offset = sum( self.cartesian_factor(k).dimension()
+                      for k in range(i) )
+        Ji = self.cartesian_factor(i)
+        Ei = Ji._module_morphism(lambda j: self.monomial(j + offset),
+                                 codomain=self)
         return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix())
 
 
@@ -3494,7 +3681,9 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
 
     SETUP::
 
-        sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+        sage: from mjo.eja.eja_algebra import (HadamardEJA,
+        ....:                                  JordanSpinEJA,
+        ....:                                  OctonionHermitianEJA,
         ....:                                  RealSymmetricEJA)
 
     EXAMPLES:
@@ -3511,30 +3700,45 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
         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` member
+    before we try to access it::
+
+        sage: J1 = OctonionHermitianEJA(1) # no rational basis
+        sage: J2 = HadamardEJA(2)
+        sage: cartesian_product([J1,J2])
+        Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
+        (+) Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
+        sage: cartesian_product([J2,J1])
+        Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
+        (+) Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
+
     """
     def __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
-            ])
+            if all( hasattr(r, "_rational_algebra") for r in algebras ):
+                self._rational_algebra = cartesian_product([
+                    r._rational_algebra for r in algebras
+                ])
 
 
 RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
 
-random_eja = ConcreteEJA.random_instance
-
-# def random_eja(*args, **kwargs):
-#     J1 = ConcreteEJA.random_instance(*args, **kwargs)
-
-#     # This might make Cartesian products appear roughly as often as
-#     # any other ConcreteEJA.
-#     if ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1) == 0:
-#         # Use random_eja() again so we can get more than two factors.
-#         J2 = random_eja(*args, **kwargs)
-#         J = cartesian_product([J1,J2])
-#         return J
-#     else:
-#         return J1
+def random_eja(*args, **kwargs):
+    J1 = ConcreteEJA.random_instance(*args, **kwargs)
+
+    # This might make Cartesian products appear roughly as often as
+    # any other ConcreteEJA.
+    if ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1) == 0:
+        # Use random_eja() again so we can get more than two factors.
+        J2 = random_eja(*args, **kwargs)
+        J = cartesian_product([J1,J2])
+        return J
+    else:
+        return J1