]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: factor out a class for real-embedded matrices.
[sage.d.git] / mjo / eja / eja_algebra.py
index 40b4bcac6583e75f6877260f633d752b6693d8d7..5bf597565b2ae292032c3f0a532763d7889cd2a8 100644 (file)
@@ -64,8 +64,7 @@ from itertools import repeat
 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
 from sage.categories.magmatic_algebras import MagmaticAlgebras
 from sage.categories.sets_cat import cartesian_product
-from sage.combinat.free_module import (CombinatorialFreeModule,
-                                       CombinatorialFreeModule_CartesianProduct)
+from sage.combinat.free_module import CombinatorialFreeModule
 from sage.matrix.constructor import matrix
 from sage.matrix.matrix_space import MatrixSpace
 from sage.misc.cachefunc import cached_method
@@ -142,19 +141,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 +152,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                 # we've specified a real embedding.
                 raise ValueError("scalar field is not real")
 
-        # If the basis given to us wasn't over the field that it's
-        # supposed to be over, fix that. Or, you know, crash.
-        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 +187,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 +205,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 +332,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
 
         """
@@ -460,9 +437,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         this algebra was constructed with ``check_axioms=False`` and
         passed an invalid multiplication table.
         """
-        return all( self.product_on_basis(i,j) == self.product_on_basis(i,j)
-                    for i in range(self.dimension())
-                    for j in range(self.dimension()) )
+        return all( x*y == y*x for x in self.gens() for y in self.gens() )
 
     def _is_jordanian(self):
         r"""
@@ -475,9 +450,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()) )
 
@@ -557,9 +532,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:
@@ -588,9 +563,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:
@@ -638,7 +613,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:
 
@@ -913,15 +888,15 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: J = JordanSpinEJA(4)
             sage: J.multiplication_table()
             +----++----+----+----+----+
-            | *  || e0 | e1 | e2 | e3 |
+            | *  || b0 | b1 | b2 | b3 |
             +====++====+====+====+====+
-            | e0 || e0 | e1 | e2 | e3 |
+            | b0 || b0 | b1 | b2 | b3 |
             +----++----+----+----+----+
-            | e1 || e1 | e0 | 0  | 0  |
+            | b1 || b1 | b0 | 0  | 0  |
             +----++----+----+----+----+
-            | e2 || e2 | 0  | e0 | 0  |
+            | b2 || b2 | 0  | b0 | 0  |
             +----++----+----+----+----+
-            | e3 || e3 | 0  | 0  | e0 |
+            | b3 || b3 | 0  | 0  | b0 |
             +----++----+----+----+----+
 
         """
@@ -931,7 +906,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.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) ]
 
@@ -975,7 +950,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: J = RealSymmetricEJA(2)
             sage: J.basis()
-            Finite family {0: e0, 1: e1, 2: e2}
+            Finite family {0: b0, 1: b1, 2: b2}
             sage: J.matrix_basis()
             (
             [1 0]  [                  0 0.7071067811865475?]  [0 0]
@@ -986,7 +961,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: J = JordanSpinEJA(2)
             sage: J.basis()
-            Finite family {0: e0, 1: e1}
+            Finite family {0: b0, 1: b1}
             sage: J.matrix_basis()
             (
             [1]  [0]
@@ -1068,20 +1043,20 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: J = HadamardEJA(5)
             sage: J.one()
-            e0 + e1 + e2 + e3 + e4
+            b0 + b1 + b2 + b3 + b4
 
         The unit element in the Hadamard EJA is inherited in the
         subalgebras generated by its elements::
 
             sage: J = HadamardEJA(5)
             sage: J.one()
-            e0 + e1 + e2 + e3 + e4
+            b0 + b1 + b2 + b3 + b4
             sage: x = sum(J.gens())
             sage: A = x.subalgebra_generated_by(orthonormalize=False)
             sage: A.one()
-            f0
+            c0
             sage: A.one().superalgebra_element()
-            e0 + e1 + e2 + e3 + e4
+            b0 + b1 + b2 + b3 + b4
 
         TESTS:
 
@@ -1433,7 +1408,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)
@@ -1750,6 +1725,21 @@ class ConcreteEJA(RationalBasisEJA):
 
 
 class MatrixEJA:
+    @staticmethod
+    def jordan_product(X,Y):
+        return (X*Y + Y*X)/2
+
+    @staticmethod
+    def trace_inner_product(X,Y):
+        r"""
+        A trace inner-product for matrices that aren't embedded in the
+        reals.
+        """
+        # We take the norm (absolute value) because Octonions() isn't
+        # smart enough yet to coerce its one() into the base field.
+        return (X*Y).trace().abs()
+
+class RealEmbeddedMatrixEJA(MatrixEJA):
     @staticmethod
     def dimension_over_reals():
         r"""
@@ -1795,9 +1785,6 @@ class MatrixEJA:
             raise ValueError("the matrix 'M' must be a real embedding")
         return M
 
-    @staticmethod
-    def jordan_product(X,Y):
-        return (X*Y + Y*X)/2
 
     @classmethod
     def trace_inner_product(cls,X,Y):
@@ -1806,29 +1793,11 @@ class MatrixEJA:
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
-            ....:                                  ComplexHermitianEJA,
+            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
             ....:                                  QuaternionHermitianEJA)
 
         EXAMPLES::
 
-        This gives the same answer as it would if we computed the trace
-        from the unembedded (original) matrices::
-
-            sage: set_random_seed()
-            sage: J = RealSymmetricEJA.random_instance()
-            sage: x,y = J.random_elements(2)
-            sage: Xe = x.to_matrix()
-            sage: Ye = y.to_matrix()
-            sage: X = J.real_unembed(Xe)
-            sage: Y = J.real_unembed(Ye)
-            sage: expected = (X*Y).trace()
-            sage: actual = J.trace_inner_product(Xe,Ye)
-            sage: actual == expected
-            True
-
-        ::
-
             sage: set_random_seed()
             sage: J = ComplexHermitianEJA.random_instance()
             sage: x,y = J.random_elements(2)
@@ -1856,27 +1825,15 @@ class MatrixEJA:
             True
 
         """
-        Xu = cls.real_unembed(X)
-        Yu = cls.real_unembed(Y)
-        tr = (Xu*Yu).trace()
-
-        try:
-            # Works in QQ, AA, RDF, et cetera.
-            return tr.real()
-        except AttributeError:
-            # A quaternion doesn't have a real() method, but does
-            # have coefficient_tuple() method that returns the
-            # coefficients of 1, i, j, and k -- in that order.
-            return tr.coefficient_tuple()[0]
-
+        # This does in fact compute the real part of the trace.
+        # If we compute the trace of e.g. a complex matrix M,
+        # then we do so by adding up its diagonal entries --
+        # call them z_1 through z_n. The real embedding of z_1
+        # will be a 2-by-2 REAL matrix [a, b; -b, a] whose trace
+        # as a REAL matrix will be 2*a = 2*Re(z_1). And so forth.
+        return (X*Y).trace()/cls.dimension_over_reals()
 
-class RealMatrixEJA(MatrixEJA):
-    @staticmethod
-    def dimension_over_reals():
-        return 1
-
-
-class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
+class RealSymmetricEJA(ConcreteEJA, MatrixEJA):
     """
     The rank-n simple EJA consisting of real symmetric n-by-n
     matrices, the usual symmetric Jordan product, and the trace inner
@@ -1889,13 +1846,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::
 
@@ -1942,7 +1899,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
 
     """
     @classmethod
-    def _denormalized_basis(cls, n):
+    def _denormalized_basis(cls, n, field):
         """
         Return a basis for the space of real symmetric n-by-n matrices.
 
@@ -1954,7 +1911,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
-            sage: B = RealSymmetricEJA._denormalized_basis(n)
+            sage: B = RealSymmetricEJA._denormalized_basis(n,ZZ)
             sage: all( M.is_symmetric() for M in  B)
             True
 
@@ -1964,7 +1921,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
         S = []
         for i in range(n):
             for j in range(i+1):
-                Eij = matrix(ZZ, n, lambda k,l: k==i and l==j)
+                Eij = matrix(field, n, lambda k,l: k==i and l==j)
                 if i == j:
                     Sij = Eij
                 else:
@@ -1985,7 +1942,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
@@ -1994,22 +1951,23 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
         if n <= 1:
             associative = True
 
-        super(RealSymmetricEJA, self).__init__(self._denormalized_basis(n),
-                                               self.jordan_product,
-                                               self.trace_inner_product,
-                                               associative=associative,
-                                               **kwargs)
+        super().__init__(self._denormalized_basis(n,field),
+                         self.jordan_product,
+                         self.trace_inner_product,
+                         field=field,
+                         associative=associative,
+                         **kwargs)
 
         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         self.rank.set_cache(n)
-        idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
+        idV = self.matrix_space().one()
         self.one.set_cache(self(idV))
 
 
 
-class ComplexMatrixEJA(MatrixEJA):
+class ComplexMatrixEJA(RealEmbeddedMatrixEJA):
     # A manual dictionary-cache for the complex_extension() method,
     # since apparently @classmethods can't also be @cached_methods.
     _complex_extension = {}
@@ -2088,7 +2046,7 @@ class ComplexMatrixEJA(MatrixEJA):
             True
 
         """
-        super(ComplexMatrixEJA,cls).real_embed(M)
+        super().real_embed(M)
         n = M.nrows()
 
         # We don't need any adjoined elements...
@@ -2135,7 +2093,7 @@ class ComplexMatrixEJA(MatrixEJA):
             True
 
         """
-        super(ComplexMatrixEJA,cls).real_unembed(M)
+        super().real_unembed(M)
         n = ZZ(M.nrows())
         d = cls.dimension_over_reals()
         F = cls.complex_extension(M.base_ring())
@@ -2216,7 +2174,7 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
     """
 
     @classmethod
-    def _denormalized_basis(cls, n):
+    def _denormalized_basis(cls, n, field):
         """
         Returns a basis for the space of complex Hermitian n-by-n matrices.
 
@@ -2234,15 +2192,14 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
-            sage: B = ComplexHermitianEJA._denormalized_basis(n)
+            sage: B = ComplexHermitianEJA._denormalized_basis(n,ZZ)
             sage: all( M.is_symmetric() for M in  B)
             True
 
         """
-        field = ZZ
-        R = PolynomialRing(field, 'z')
+        R = PolynomialRing(ZZ, 'z')
         z = R.gen()
-        F = field.extension(z**2 + 1, 'I')
+        F = ZZ.extension(z**2 + 1, 'I')
         I = F.gen(1)
 
         # This is like the symmetric case, but we need to be careful:
@@ -2273,12 +2230,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
@@ -2287,11 +2244,12 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
         if n <= 1:
             associative = True
 
-        super(ComplexHermitianEJA, self).__init__(self._denormalized_basis(n),
-                                                  self.jordan_product,
-                                                  self.trace_inner_product,
-                                                  associative=associative,
-                                                  **kwargs)
+        super().__init__(self._denormalized_basis(n,field),
+                         self.jordan_product,
+                         self.trace_inner_product,
+                         field=field,
+                         associative=associative,
+                         **kwargs)
         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
@@ -2311,7 +2269,7 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
         return cls(n, **kwargs)
 
-class QuaternionMatrixEJA(MatrixEJA):
+class QuaternionMatrixEJA(RealEmbeddedMatrixEJA):
 
     # A manual dictionary-cache for the quaternion_extension() method,
     # since apparently @classmethods can't also be @cached_methods.
@@ -2374,7 +2332,7 @@ class QuaternionMatrixEJA(MatrixEJA):
             True
 
         """
-        super(QuaternionMatrixEJA,cls).real_embed(M)
+        super().real_embed(M)
         quaternions = M.base_ring()
         n = M.nrows()
 
@@ -2429,7 +2387,7 @@ class QuaternionMatrixEJA(MatrixEJA):
             True
 
         """
-        super(QuaternionMatrixEJA,cls).real_unembed(M)
+        super().real_unembed(M)
         n = ZZ(M.nrows())
         d = cls.dimension_over_reals()
 
@@ -2517,7 +2475,7 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
 
     """
     @classmethod
-    def _denormalized_basis(cls, n):
+    def _denormalized_basis(cls, n, field):
         """
         Returns a basis for the space of quaternion Hermitian n-by-n matrices.
 
@@ -2535,12 +2493,11 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
-            sage: B = QuaternionHermitianEJA._denormalized_basis(n)
+            sage: B = QuaternionHermitianEJA._denormalized_basis(n,ZZ)
             sage: all( M.is_symmetric() for M in B )
             True
 
         """
-        field = ZZ
         Q = QuaternionAlgebra(QQ,-1,-1)
         I,J,K = Q.gens()
 
@@ -2584,12 +2541,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
@@ -2598,11 +2555,12 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
         if n <= 1:
             associative = True
 
-        super(QuaternionHermitianEJA, self).__init__(self._denormalized_basis(n),
-                                                     self.jordan_product,
-                                                     self.trace_inner_product,
-                                                     associative=associative,
-                                                     **kwargs)
+        super().__init__(self._denormalized_basis(n,field),
+                         self.jordan_product,
+                         self.trace_inner_product,
+                         field=field,
+                         associative=associative,
+                         **kwargs)
 
         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
@@ -2646,19 +2604,19 @@ class HadamardEJA(ConcreteEJA):
     This multiplication table can be verified by hand::
 
         sage: J = HadamardEJA(3)
-        sage: e0,e1,e2 = J.gens()
-        sage: e0*e0
-        e0
-        sage: e0*e1
+        sage: b0,b1,b2 = J.gens()
+        sage: b0*b0
+        b0
+        sage: b0*b1
         0
-        sage: e0*e2
+        sage: b0*b2
         0
-        sage: e1*e1
-        e1
-        sage: e1*e2
+        sage: b1*b1
+        b1
+        sage: b1*b2
         0
-        sage: e2*e2
-        e2
+        sage: b2*b2
+        b2
 
     TESTS:
 
@@ -2668,7 +2626,7 @@ class HadamardEJA(ConcreteEJA):
         (r0, r1, r2)
 
     """
-    def __init__(self, n, **kwargs):
+    def __init__(self, n, field=AA, **kwargs):
         if n == 0:
             jordan_product = lambda x,y: x
             inner_product = lambda x,y: x
@@ -2689,10 +2647,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)
@@ -2800,7 +2760,7 @@ class BilinearFormEJA(ConcreteEJA):
         True
 
     """
-    def __init__(self, B, **kwargs):
+    def __init__(self, B, field=AA, **kwargs):
         # The matrix "B" is supplied by the user in most cases,
         # so it makes sense to check whether or not its positive-
         # definite unless we are specifically asked not to...
@@ -2828,18 +2788,20 @@ class BilinearFormEJA(ConcreteEJA):
             return P([z0] + zbar.list())
 
         n = B.nrows()
-        column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )
+        column_basis = tuple( b.column()
+                              for b in FreeModule(field, n).basis() )
 
         # TODO: I haven't actually checked this, but it seems legit.
         associative = False
         if n <= 2:
             associative = True
 
-        super(BilinearFormEJA, self).__init__(column_basis,
-                                              jordan_product,
-                                              inner_product,
-                                              associative=associative,
-                                              **kwargs)
+        super().__init__(column_basis,
+                         jordan_product,
+                         inner_product,
+                         field=field,
+                         associative=associative,
+                         **kwargs)
 
         # The rank of this algebra is two, unless we're in a
         # one-dimensional ambient space (because the rank is bounded
@@ -2899,20 +2861,20 @@ class JordanSpinEJA(BilinearFormEJA):
     This multiplication table can be verified by hand::
 
         sage: J = JordanSpinEJA(4)
-        sage: e0,e1,e2,e3 = J.gens()
-        sage: e0*e0
-        e0
-        sage: e0*e1
-        e1
-        sage: e0*e2
-        e2
-        sage: e0*e3
-        e3
-        sage: e1*e2
+        sage: b0,b1,b2,b3 = J.gens()
+        sage: b0*b0
+        b0
+        sage: b0*b1
+        b1
+        sage: b0*b2
+        b2
+        sage: b0*b3
+        b3
+        sage: b1*b2
         0
-        sage: e1*e3
+        sage: b1*b3
         0
-        sage: e2*e3
+        sage: b2*b3
         0
 
     We can change the generator prefix::
@@ -2933,7 +2895,7 @@ class JordanSpinEJA(BilinearFormEJA):
             True
 
     """
-    def __init__(self, n, **kwargs):
+    def __init__(self, n, *args, **kwargs):
         # This is a special case of the BilinearFormEJA with the
         # identity matrix as its bilinear form.
         B = matrix.identity(ZZ, n)
@@ -2944,7 +2906,7 @@ class JordanSpinEJA(BilinearFormEJA):
 
         # But also don't pass check_field=False here, because the user
         # can pass in a field!
-        super(JordanSpinEJA, self).__init__(B, **kwargs)
+        super().__init__(B, *args, **kwargs)
 
     @staticmethod
     def _max_random_instance_size():
@@ -3002,11 +2964,11 @@ class TrivialEJA(ConcreteEJA):
         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
-        super(TrivialEJA, self).__init__(basis,
-                                         jordan_product,
-                                         inner_product,
-                                         associative=True,
-                                         **kwargs)
+        super().__init__(basis,
+                         jordan_product,
+                         inner_product,
+                         associative=True,
+                         **kwargs)
 
         # The rank is zero using my definition, namely the dimension of the
         # largest subalgebra generated by any element.
@@ -3020,8 +2982,7 @@ class TrivialEJA(ConcreteEJA):
         return cls(**kwargs)
 
 
-class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
-                          FiniteDimensionalEJA):
+class CartesianProductEJA(FiniteDimensionalEJA):
     r"""
     The external (orthogonal) direct sum of two or more Euclidean
     Jordan algebras. Every Euclidean Jordan algebra decomposes into an
@@ -3117,6 +3078,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::
@@ -3144,37 +3132,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
@@ -3194,9 +3186,25 @@ 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))
+        ones = tuple(J.one().to_matrix() for J in factors)
+        self.one.set_cache(self(ones))
+        self.rank.set_cache(sum(J.rank() for J in factors))
+
+    def cartesian_factors(self):
+        # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
+        return self._sets
+
+    def cartesian_factor(self, i):
+        r"""
+        Return the ``i``th factor of this algebra.
+        """
+        return self._sets[i]
+
+    def _repr_(self):
+        # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
+        from sage.categories.cartesian_product import cartesian_product
+        return cartesian_product.symbol.join("%s" % factor
+                                             for factor in self._sets)
 
     def matrix_space(self):
         r"""
@@ -3293,9 +3301,12 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
             True
 
         """
-        Ji = self.cartesian_factors()[i]
-        # Requires the fix on Trac 31421/31422 to work!
-        Pi = super().cartesian_projection(i)
+        offset = sum( self.cartesian_factor(k).dimension()
+                      for k in range(i) )
+        Ji = self.cartesian_factor(i)
+        Pi = self._module_morphism(lambda j: Ji.monomial(j - offset),
+                                   codomain=Ji)
+
         return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix())
 
     @cached_method
@@ -3401,9 +3412,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())
 
 
@@ -3448,4 +3461,15 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
 
 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