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

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

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index 031359fbfdab4ea8dc852f24ad7cfea8e54bfe4d..5bf597565b2ae292032c3f0a532763d7889cd2a8 100644 (file)
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -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.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


@@ -142,19 +141,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                  cartesian_product=False,
                  check_field=True,
                  check_axioms=True,
                  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)
         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):
 
         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")
 
                 # 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
         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:
             # 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.
 
         # 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
         # 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
         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: 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)
             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
 
         """


@@ -473,9 +450,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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()) )
 


@@ -555,9 +532,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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)*z - x*(y*z)
 
                     if diff.norm() > epsilon:
                     diff = (x*y)*z - x*(y*z)
 
                     if diff.norm() > epsilon:


@@ -586,9 +563,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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)
 
                     if diff.abs() > epsilon:
                     diff = (x*y).inner_product(z) - x.inner_product(y*z)
 
                     if diff.abs() > epsilon:


@@ -636,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]) )
             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:
 
 
         TESTS:
 


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

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

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

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

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

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

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

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

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

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

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


@@ -929,7 +906,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.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) ]
 
                                     for j in range(n) ]
                for i in range(n) ]
 


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

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

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


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

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

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


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

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

 
         The unit element in the Hadamard EJA is inherited in the
         subalgebras generated by its elements::
 
             sage: J = HadamardEJA(5)
             sage: J.one()
 
         The unit element in the Hadamard EJA is inherited in the
         subalgebras generated by its elements::
 
             sage: J = HadamardEJA(5)
             sage: J.one()

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

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

-            f0
+            c0

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

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

 
         TESTS:
 
 
         TESTS:
 


@@ -1431,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.
         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)
                         for k in range(n) )
 
         L_x = matrix(F, n, n, L_x_i_j)


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

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

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


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

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

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


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

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

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

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

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


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

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

 
 

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

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


@@ -1887,13 +1846,13 @@ 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::
 


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

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

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


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

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

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


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

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

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


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

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

         # We know this is a valid EJA, but will double-check
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
         # We know this is a valid EJA, but will double-check
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False


@@ -1992,9 +1951,10 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
         if n <= 1:
             associative = True
 
         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,
                          self.jordan_product,
                          self.trace_inner_product,

+                         field=field,

                          associative=associative,
                          **kwargs)
 
                          associative=associative,
                          **kwargs)
 


@@ -2002,12 +1962,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)
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         self.rank.set_cache(n)

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

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

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

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


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

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

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


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

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

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

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

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

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

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


@@ -2271,12 +2230,12 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
                 # "erase" E_ij
                 Eij[i,j] = 0
 
                 # "erase" E_ij
                 Eij[i,j] = 0
 

-        # Since we embedded these, we can drop back to the "field" that we
-        # started with instead of the complex extension "F".
+        # Since we embedded the entries, we can drop back to the
+        # desired real "field" instead of the extension "F".

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

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

         # We know this is a valid EJA, but will double-check
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
         # We know this is a valid EJA, but will double-check
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False


@@ -2285,9 +2244,10 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
         if n <= 1:
             associative = True
 
         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,
                          self.jordan_product,
                          self.trace_inner_product,

+                         field=field,

                          associative=associative,
                          **kwargs)
         # TODO: this could be factored out somehow, but is left here
                          associative=associative,
                          **kwargs)
         # TODO: this could be factored out somehow, but is left here


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

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

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


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

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

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


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

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

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

-        field = ZZ

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


@@ -2582,12 +2541,12 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
                 # "erase" E_ij
                 Eij[i,j] = 0
 
                 # "erase" E_ij
                 Eij[i,j] = 0
 

-        # Since we embedded these, we can drop back to the "field" that we
-        # started with instead of the quaternion algebra "Q".
+        # Since we embedded the entries, we can drop back to the
+        # desired real "field" instead of the quaternion algebra "Q".

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

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

         # We know this is a valid EJA, but will double-check
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
         # We know this is a valid EJA, but will double-check
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False


@@ -2596,9 +2555,10 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
         if n <= 1:
             associative = True
 
         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,
                          self.jordan_product,
                          self.trace_inner_product,

+                         field=field,

                          associative=associative,
                          **kwargs)
 
                          associative=associative,
                          **kwargs)
 


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

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

         0
         0

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

         0
         0

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

         0
         0

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

 
     TESTS:
 
 
     TESTS:
 


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

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

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


@@ -2687,10 +2647,12 @@ 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( b.column()
+                              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,

                          associative=True,
                          **kwargs)
         self.rank.set_cache(n)
                          associative=True,
                          **kwargs)
         self.rank.set_cache(n)


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

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

         # The matrix "B" is supplied by the user in most cases,
         # so it makes sense to check whether or not its positive-
         # definite unless we are specifically asked not to...
         # The matrix "B" is supplied by the user in most cases,
         # so it makes sense to check whether or not its positive-
         # definite unless we are specifically asked not to...


@@ -2826,7 +2788,8 @@ class BilinearFormEJA(ConcreteEJA):
             return P([z0] + zbar.list())
 
         n = B.nrows()
             return P([z0] + zbar.list())
 
         n = B.nrows()

-        column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )
+        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
 
         # TODO: I haven't actually checked this, but it seems legit.
         associative = False


@@ -2836,6 +2799,7 @@ class BilinearFormEJA(ConcreteEJA):
         super().__init__(column_basis,
                          jordan_product,
                          inner_product,
         super().__init__(column_basis,
                          jordan_product,
                          inner_product,

+                         field=field,

                          associative=associative,
                          **kwargs)
 
                          associative=associative,
                          **kwargs)
 


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

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

         0
         0

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

         0
         0

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

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


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

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

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


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

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

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


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

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

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


@@ -3115,6 +3078,33 @@ 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 |
+        +----++----+----+----+
+

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


@@ -3142,37 +3132,41 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
     Element = FiniteDimensionalEJAElement
 
 
     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")
 
             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()
         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(
 
         # 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(
             ))
 
         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
             )
 
         # There's no need to check the field since it already came


@@ -3192,87 +3186,25 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
                                       check_field=False,
                                       check_axioms=False)
 
                                       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()
+        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))

 
 

-        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])
+    def cartesian_factors(self):
+        # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
+        return self._sets

 
 

-        offset = sum( f.dimension()
-                      for f in self.cartesian_factors()[:factor] )
-        return offset + idx_in_factor
-
-    def product_on_basis(self, i, j):
+    def cartesian_factor(self, i):

         r"""
         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"""
 
     def matrix_space(self):
         r"""


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

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

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


@@ -3477,9 +3412,11 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
             True
 
         """
             True
 
         """

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

         return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix())
 
 
         return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix())
 
 


@@ -3524,17 +3461,15 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
 
 RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
 
 
 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