eja: cache the charpoly coefficients for the AlbertEJA.

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

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index 5d96a53f402f4f3343cc396049b4311d13fa3ba1..587d8e339463adaafb0a390ec164e76a6e5ca0c4 100644 (file)
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -32,22 +32,25 @@ for these simple algebras:
   * :class:`RealSymmetricEJA`
   * :class:`ComplexHermitianEJA`
   * :class:`QuaternionHermitianEJA`
   * :class:`RealSymmetricEJA`
   * :class:`ComplexHermitianEJA`
   * :class:`QuaternionHermitianEJA`

+  * :class:`OctonionHermitianEJA`

 
 

-Missing from this list is the algebra of three-by-three octononion
-Hermitian matrices, as there is (as of yet) no implementation of the
-octonions in SageMath. In addition to these, we provide two other
-example constructions,
+In addition to these, we provide two other example constructions,

 
 

+  * :class:`JordanSpinEJA`

   * :class:`HadamardEJA`
   * :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
   * :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::
 
 
 SETUP::
 


@@ -59,13 +62,10 @@ EXAMPLES::
     Euclidean Jordan algebra of dimension...
 """
 
     Euclidean Jordan algebra of dimension...
 """
 

-from itertools import repeat
-

 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
 from sage.categories.magmatic_algebras import MagmaticAlgebras
 from sage.categories.sets_cat import cartesian_product
 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
 from sage.categories.magmatic_algebras import MagmaticAlgebras
 from sage.categories.sets_cat import cartesian_product

-from sage.combinat.free_module import (CombinatorialFreeModule,
-                                       CombinatorialFreeModule_CartesianProduct)
+from sage.combinat.free_module import CombinatorialFreeModule

 from sage.matrix.constructor import matrix
 from sage.matrix.matrix_space import MatrixSpace
 from sage.misc.cachefunc import cached_method
 from sage.matrix.constructor import matrix
 from sage.matrix.matrix_space import MatrixSpace
 from sage.misc.cachefunc import cached_method


@@ -104,6 +104,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         product. This will be applied to ``basis`` to compute an
         inner-product table (basically a matrix) for this algebra.
 
         product. This will be applied to ``basis`` to compute an
         inner-product table (basically a matrix) for this algebra.
 

+      - ``matrix_space`` -- the space that your matrix basis lives in,
+        or ``None`` (the default). So long as your basis does not have
+        length zero you can omit this. But in trivial algebras, it is
+        required.
+

       - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
         field for the algebra.
 
       - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
         field for the algebra.
 


@@ -128,7 +133,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         sage: basis = tuple(b.superalgebra_element() for b in A.basis())
         sage: J.subalgebra(basis, orthonormalize=False).is_associative()
         True
         sage: basis = tuple(b.superalgebra_element() for b in A.basis())
         sage: J.subalgebra(basis, orthonormalize=False).is_associative()
         True

-

     """
     Element = FiniteDimensionalEJAElement
 
     """
     Element = FiniteDimensionalEJAElement
 


@@ -137,24 +141,15 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                  jordan_product,
                  inner_product,
                  field=AA,
                  jordan_product,
                  inner_product,
                  field=AA,

+                 matrix_space=None,

                  orthonormalize=True,
                  associative=None,
                  cartesian_product=False,
                  check_field=True,
                  check_axioms=True,
                  orthonormalize=True,
                  associative=None,
                  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 +158,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


@@ -197,6 +177,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         category = MagmaticAlgebras(field).FiniteDimensional()
         category = category.WithBasis().Unital().Commutative()
 
         category = MagmaticAlgebras(field).FiniteDimensional()
         category = category.WithBasis().Unital().Commutative()
 

+        if n <= 1:
+            # All zero- and one-dimensional algebras are just the real
+            # numbers with (some positive multiples of) the usual
+            # multiplication as its Jordan and inner-product.
+            associative = True

         if associative is None:
             # We should figure it out. As with check_axioms, we have to do
             # this without the help of the _jordan_product_is_associative()
         if associative is None:
             # We should figure it out. As with check_axioms, we have to do
             # this without the help of the _jordan_product_is_associative()


@@ -213,7 +198,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 +216,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


@@ -253,8 +241,14 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             basis = tuple(gram_schmidt(basis, inner_product))
 
         # Save the (possibly orthonormalized) matrix basis for
             basis = tuple(gram_schmidt(basis, inner_product))
 
         # Save the (possibly orthonormalized) matrix basis for

-        # later...
+        # later, as well as the space that its elements live in.
+        # In most cases we can deduce the matrix space, but when
+        # n == 0 (that is, there are no basis elements) we cannot.

         self._matrix_basis = basis
         self._matrix_basis = basis

+        if matrix_space is None:
+            self._matrix_space = self._matrix_basis[0].parent()
+        else:
+            self._matrix_space = matrix_space

 
         # Now create the vector space for the algebra, which will have
         # its own set of non-ambient coordinates (in terms of the
 
         # Now create the vector space for the algebra, which will have
         # its own set of non-ambient coordinates (in terms of the


@@ -355,16 +349,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 +467,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 +549,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 +580,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 +630,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 +905,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 +923,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 +967,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 +978,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]


@@ -1038,16 +1032,16 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
             sage: J.matrix_space()
 
             sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
             sage: J.matrix_space()

-            Full MatrixSpace of 4 by 4 dense matrices over Rational Field
+            Module of 2 by 2 matrices with entries in Algebraic Field over
+            the scalar ring Rational Field

             sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
             sage: J.matrix_space()
             sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
             sage: J.matrix_space()

-            Full MatrixSpace of 4 by 4 dense matrices over Rational Field
+            Module of 1 by 1 matrices with entries in Quaternion
+            Algebra (-1, -1) with base ring Rational Field over
+            the scalar ring Rational Field

 
         """
 
         """

-        if self.is_trivial():
-            return MatrixSpace(self.base_ring(), 0)
-        else:
-            return self.matrix_basis()[0].parent()
+        return self._matrix_space

 
 
     @cached_method
 
 
     @cached_method


@@ -1066,20 +1060,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 +1425,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)


@@ -1566,7 +1560,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
 class RationalBasisEJA(FiniteDimensionalEJA):
     r"""
 
 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::
 
 
     SETUP::
 


@@ -1597,7 +1591,11 @@ class RationalBasisEJA(FiniteDimensionalEJA):
         if check_field:
             # Abuse the check_field parameter to check that the entries of
             # out basis (in ambient coordinates) are in the field QQ.
         if check_field:
             # Abuse the check_field parameter to check that the entries of
             # out basis (in ambient coordinates) are in the field QQ.

-            if not all( all(b_i in QQ for b_i in b.list()) for b in basis ):
+            # Use _all2list to get the vector coordinates of octonion
+            # entries and not the octonions themselves (which are not
+            # rational).
+            if not all( all(b_i in QQ for b_i in _all2list(b))
+                        for b in basis ):

                 raise TypeError("basis not rational")
 
         super().__init__(basis,
                 raise TypeError("basis not rational")
 
         super().__init__(basis,


@@ -1620,6 +1618,7 @@ class RationalBasisEJA(FiniteDimensionalEJA):
                                        jordan_product,
                                        inner_product,
                                        field=QQ,
                                        jordan_product,
                                        inner_product,
                                        field=QQ,

+                                       matrix_space=self.matrix_space(),

                                        associative=self.is_associative(),
                                        orthonormalize=False,
                                        check_field=False,
                                        associative=self.is_associative(),
                                        orthonormalize=False,
                                        check_field=False,


@@ -1679,7 +1678,7 @@ class RationalBasisEJA(FiniteDimensionalEJA):
         subs_dict = { X[i]: BX[i] for i in range(len(X)) }
         return tuple( a_i.subs(subs_dict) for a_i in a )
 
         subs_dict = { X[i]: BX[i] for i in range(len(X)) }
         return tuple( a_i.subs(subs_dict) for a_i in a )
 

-class ConcreteEJA(RationalBasisEJA):
+class ConcreteEJA(FiniteDimensionalEJA):

     r"""
     A class for the Euclidean Jordan algebras that we know by name.
 
     r"""
     A class for the Euclidean Jordan algebras that we know by name.
 


@@ -1747,134 +1746,164 @@ class ConcreteEJA(RationalBasisEJA):
         return eja_class.random_instance(*args, **kwargs)
 
 
         return eja_class.random_instance(*args, **kwargs)
 
 

-class MatrixEJA:
+class MatrixEJA(FiniteDimensionalEJA):

     @staticmethod
     @staticmethod

-    def dimension_over_reals():
-        r"""
-        The dimension of this matrix's base ring over the reals.
+    def _denormalized_basis(A):
+        """
+        Returns a basis for the space of complex Hermitian n-by-n matrices.

 
 

-        The reals are dimension one over themselves, obviously; that's
-        just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
-        have dimension two. Finally, the quaternions have dimension
-        four over the reals.
+        Why do we embed these? Basically, because all of numerical linear
+        algebra assumes that you're working with vectors consisting of `n`
+        entries from a field and scalars from the same field. There's no way
+        to tell SageMath that (for example) the vectors contain complex
+        numbers, while the scalar field is real.

 
 

-        This is used to determine the size of the matrix returned from
-        :meth:`real_embed`, among other things.
-        """
-        raise NotImplementedError
+        SETUP::

 
 

-    @classmethod
-    def real_embed(cls,M):
-        """
-        Embed the matrix ``M`` into a space of real matrices.
+            sage: from mjo.hurwitz import (ComplexMatrixAlgebra,
+            ....:                          QuaternionMatrixAlgebra,
+            ....:                          OctonionMatrixAlgebra)
+            sage: from mjo.eja.eja_algebra import MatrixEJA
+
+        TESTS::

 
 

-        The matrix ``M`` can have entries in any field at the moment:
-        the real numbers, complex numbers, or quaternions. And although
-        they are not a field, we can probably support octonions at some
-        point, too. This function returns a real matrix that "acts like"
-        the original with respect to matrix multiplication; i.e.
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: A = MatrixSpace(QQ, n)
+            sage: B = MatrixEJA._denormalized_basis(A)
+            sage: all( M.is_hermitian() for M in  B)
+            True

 
 

-          real_embed(M*N) = real_embed(M)*real_embed(N)
+        ::

 
 

-        """
-        if M.ncols() != M.nrows():
-            raise ValueError("the matrix 'M' must be square")
-        return M
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: A = ComplexMatrixAlgebra(n, scalars=QQ)
+            sage: B = MatrixEJA._denormalized_basis(A)
+            sage: all( M.is_hermitian() for M in  B)
+            True

 
 

+        ::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: A = QuaternionMatrixAlgebra(n, scalars=QQ)
+            sage: B = MatrixEJA._denormalized_basis(A)
+            sage: all( M.is_hermitian() for M in B )
+            True
+
+        ::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: A = OctonionMatrixAlgebra(n, scalars=QQ)
+            sage: B = MatrixEJA._denormalized_basis(A)
+            sage: all( M.is_hermitian() for M in B )
+            True

 
 

-    @classmethod
-    def real_unembed(cls,M):
-        """
-        The inverse of :meth:`real_embed`.

         """
         """

-        if M.ncols() != M.nrows():
-            raise ValueError("the matrix 'M' must be square")
-        if not ZZ(M.nrows()).mod(cls.dimension_over_reals()).is_zero():
-            raise ValueError("the matrix 'M' must be a real embedding")
-        return M
+        # These work for real MatrixSpace, whose monomials only have
+        # two coordinates (because the last one would always be "1").
+        es = A.base_ring().gens()
+        gen = lambda A,m: A.monomial(m[:2])
+
+        if hasattr(A, 'entry_algebra_gens'):
+            # We've got a MatrixAlgebra, and its monomials will have
+            # three coordinates.
+            es = A.entry_algebra_gens()
+            gen = lambda A,m: A.monomial(m)
+
+        basis = []
+        for i in range(A.nrows()):
+            for j in range(i+1):
+                if i == j:
+                    E_ii = gen(A, (i,j,es[0]))
+                    basis.append(E_ii)
+                else:
+                    for e in es:
+                        E_ij  = gen(A, (i,j,e))
+                        E_ij += E_ij.conjugate_transpose()
+                        basis.append(E_ij)
+
+        return tuple( basis )

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

-    @classmethod
-    def trace_inner_product(cls,X,Y):
+    @staticmethod
+    def trace_inner_product(X,Y):

         r"""
         r"""

-        Compute the trace inner-product of two real-embeddings.
+        A trace inner-product for matrices that aren't embedded in the
+        reals. It takes MATRICES as arguments, not EJA elements.

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

-            ....:                                  QuaternionHermitianEJA)
+            ....:                                  QuaternionHermitianEJA,
+            ....:                                  OctonionHermitianEJA)

 
         EXAMPLES::
 
 
         EXAMPLES::
 

-        This gives the same answer as it would if we computed the trace
-        from the unembedded (original) matrices::
+            sage: J = RealSymmetricEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2

 
 

-            sage: set_random_seed()
-            sage: J = RealSymmetricEJA.random_instance()
-            sage: x,y = J.random_elements(2)
-            sage: Xe = x.to_matrix()
-            sage: Ye = y.to_matrix()
-            sage: X = J.real_unembed(Xe)
-            sage: Y = J.real_unembed(Ye)
-            sage: expected = (X*Y).trace()
-            sage: actual = J.trace_inner_product(Xe,Ye)
-            sage: actual == expected
-            True
+        ::
+
+            sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2

 
         ::
 
 
         ::
 

-            sage: set_random_seed()
-            sage: J = ComplexHermitianEJA.random_instance()
-            sage: x,y = J.random_elements(2)
-            sage: Xe = x.to_matrix()
-            sage: Ye = y.to_matrix()
-            sage: X = J.real_unembed(Xe)
-            sage: Y = J.real_unembed(Ye)
-            sage: expected = (X*Y).trace().real()
-            sage: actual = J.trace_inner_product(Xe,Ye)
-            sage: actual == expected
-            True
+            sage: J = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2

 
         ::
 
 
         ::
 

-            sage: set_random_seed()
-            sage: J = QuaternionHermitianEJA.random_instance()
-            sage: x,y = J.random_elements(2)
-            sage: Xe = x.to_matrix()
-            sage: Ye = y.to_matrix()
-            sage: X = J.real_unembed(Xe)
-            sage: Y = J.real_unembed(Ye)
-            sage: expected = (X*Y).trace().coefficient_tuple()[0]
-            sage: actual = J.trace_inner_product(Xe,Ye)
-            sage: actual == expected
-            True
+            sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2

 
         """
 
         """

-        Xu = cls.real_unembed(X)
-        Yu = cls.real_unembed(Y)
-        tr = (Xu*Yu).trace()
-
-        try:
-            # Works in QQ, AA, RDF, et cetera.
-            return tr.real()
-        except AttributeError:
-            # A quaternion doesn't have a real() method, but does
-            # have coefficient_tuple() method that returns the
-            # coefficients of 1, i, j, and k -- in that order.
+        tr = (X*Y).trace()
+        if hasattr(tr, 'coefficient'):
+            # Works for octonions, and has to come first because they
+            # also have a "real()" method that doesn't return an
+            # element of the scalar ring.
+            return tr.coefficient(0)
+        elif hasattr(tr, 'coefficient_tuple'):
+            # Works for quaternions.

             return tr.coefficient_tuple()[0]
 
             return tr.coefficient_tuple()[0]
 

+        # Works for real and complex numbers.
+        return tr.real()

 
 

-class RealMatrixEJA(MatrixEJA):
-    @staticmethod
-    def dimension_over_reals():
-        return 1

 
 

+    def __init__(self, matrix_space, **kwargs):
+        # We know this is a valid EJA, but will double-check
+        # if the user passes check_axioms=True.
+        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
+
+
+        super().__init__(self._denormalized_basis(matrix_space),
+                         self.jordan_product,
+                         self.trace_inner_product,
+                         field=matrix_space.base_ring(),
+                         matrix_space=matrix_space,
+                         **kwargs)
+
+        self.rank.set_cache(matrix_space.nrows())
+        self.one.set_cache( self(matrix_space.one()) )

 
 

-class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
+class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):

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


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


@@ -1939,38 +1968,6 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """

-    @classmethod
-    def _denormalized_basis(cls, n):
-        """
-        Return a basis for the space of real symmetric n-by-n matrices.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import RealSymmetricEJA
-
-        TESTS::
-
-            sage: set_random_seed()
-            sage: n = ZZ.random_element(1,5)
-            sage: B = RealSymmetricEJA._denormalized_basis(n)
-            sage: all( M.is_symmetric() for M in  B)
-            True
-
-        """
-        # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
-        # coordinates.
-        S = []
-        for i in range(n):
-            for j in range(i+1):
-                Eij = matrix(ZZ, n, lambda k,l: k==i and l==j)
-                if i == j:
-                    Sij = Eij
-                else:
-                    Sij = Eij + Eij.transpose()
-                S.append(Sij)
-        return tuple(S)
-
-

     @staticmethod
     def _max_random_instance_size():
         return 4 # Dimension 10
     @staticmethod
     def _max_random_instance_size():
         return 4 # Dimension 10


@@ -1983,179 +1980,17 @@ 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
 

-        associative = False
-        if n <= 1:
-            associative = True
-
-        super().__init__(self._denormalized_basis(n),
-                         self.jordan_product,
-                         self.trace_inner_product,
-                         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)
-        self.one.set_cache(self(idV))
-
-
-
-class ComplexMatrixEJA(MatrixEJA):
-    # A manual dictionary-cache for the complex_extension() method,
-    # since apparently @classmethods can't also be @cached_methods.
-    _complex_extension = {}
-
-    @classmethod
-    def complex_extension(cls,field):
-        r"""
-        The complex field that we embed/unembed, as an extension
-        of the given ``field``.
-        """
-        if field in cls._complex_extension:
-            return cls._complex_extension[field]
-
-        # Sage doesn't know how to adjoin the complex "i" (the root of
-        # x^2 + 1) to a field in a general way. Here, we just enumerate
-        # all of the cases that I have cared to support so far.
-        if field is AA:
-            # Sage doesn't know how to embed AA into QQbar, i.e. how
-            # to adjoin sqrt(-1) to AA.
-            F = QQbar
-        elif not field.is_exact():
-            # RDF or RR
-            F = field.complex_field()
-        else:
-            # Works for QQ and... maybe some other fields.
-            R = PolynomialRing(field, 'z')
-            z = R.gen()
-            F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
-
-        cls._complex_extension[field] = F
-        return F
-
-    @staticmethod
-    def dimension_over_reals():
-        return 2
-
-    @classmethod
-    def real_embed(cls,M):
-        """
-        Embed the n-by-n complex matrix ``M`` into the space of real
-        matrices of size 2n-by-2n via the map the sends each entry `z = a +
-        bi` to the block matrix ``[[a,b],[-b,a]]``.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
-
-        EXAMPLES::
-
-            sage: F = QuadraticField(-1, 'I')
-            sage: x1 = F(4 - 2*i)
-            sage: x2 = F(1 + 2*i)
-            sage: x3 = F(-i)
-            sage: x4 = F(6)
-            sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
-            sage: ComplexMatrixEJA.real_embed(M)
-            [ 4 -2| 1  2]
-            [ 2  4|-2  1]
-            [-----+-----]
-            [ 0 -1| 6  0]
-            [ 1  0| 0  6]
-
-        TESTS:
-
-        Embedding is a homomorphism (isomorphism, in fact)::
-
-            sage: set_random_seed()
-            sage: n = ZZ.random_element(3)
-            sage: F = QuadraticField(-1, 'I')
-            sage: X = random_matrix(F, n)
-            sage: Y = random_matrix(F, n)
-            sage: Xe = ComplexMatrixEJA.real_embed(X)
-            sage: Ye = ComplexMatrixEJA.real_embed(Y)
-            sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
-            sage: Xe*Ye == XYe
-            True
-
-        """
-        super().real_embed(M)
-        n = M.nrows()
-
-        # We don't need any adjoined elements...
-        field = M.base_ring().base_ring()
-
-        blocks = []
-        for z in M.list():
-            a = z.real()
-            b = z.imag()
-            blocks.append(matrix(field, 2, [ [ a, b],
-                                             [-b, a] ]))
-
-        return matrix.block(field, n, blocks)
-
-
-    @classmethod
-    def real_unembed(cls,M):
-        """
-        The inverse of _embed_complex_matrix().
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
-
-        EXAMPLES::
-
-            sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
-            ....:                 [-2,  1,  -4,  3],
-            ....:                 [ 9,  10, 11, 12],
-            ....:                 [-10, 9, -12, 11] ])
-            sage: ComplexMatrixEJA.real_unembed(A)
-            [  2*I + 1   4*I + 3]
-            [ 10*I + 9 12*I + 11]
-
-        TESTS:
+        A = MatrixSpace(field, n)
+        super().__init__(A, **kwargs)

 
 

-        Unembedding is the inverse of embedding::

 
 

-            sage: set_random_seed()
-            sage: F = QuadraticField(-1, 'I')
-            sage: M = random_matrix(F, 3)
-            sage: Me = ComplexMatrixEJA.real_embed(M)
-            sage: ComplexMatrixEJA.real_unembed(Me) == M
-            True
-
-        """
-        super().real_unembed(M)
-        n = ZZ(M.nrows())
-        d = cls.dimension_over_reals()
-        F = cls.complex_extension(M.base_ring())
-        i = F.gen()

 
 

-        # Go top-left to bottom-right (reading order), converting every
-        # 2-by-2 block we see to a single complex element.
-        elements = []
-        for k in range(n/d):
-            for j in range(n/d):
-                submat = M[d*k:d*k+d,d*j:d*j+d]
-                if submat[0,0] != submat[1,1]:
-                    raise ValueError('bad on-diagonal submatrix')
-                if submat[0,1] != -submat[1,0]:
-                    raise ValueError('bad off-diagonal submatrix')
-                z = submat[0,0] + submat[0,1]*i
-                elements.append(z)
-
-        return matrix(F, n/d, elements)
-
-
-class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
+class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):

     """
     The rank-n simple EJA consisting of complex Hermitian n-by-n
     matrices over the real numbers, the usual symmetric Jordan product,
     """
     The rank-n simple EJA consisting of complex Hermitian n-by-n
     matrices over the real numbers, the usual symmetric Jordan product,


@@ -2168,13 +2003,28 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
 
     EXAMPLES:
 
 
     EXAMPLES:
 

-    In theory, our "field" can be any subfield of the reals::
+    In theory, our "field" can be any subfield of the reals, but we
+    can't use inexact real fields at the moment because SageMath
+    doesn't know how to convert their elements into complex numbers,
+    or even into algebraic reals::

 
 

-        sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True)
-        Euclidean Jordan algebra of dimension 4 over Real Double Field
-        sage: ComplexHermitianEJA(2, field=RR, check_axioms=True)
-        Euclidean Jordan algebra of dimension 4 over Real Field with
-        53 bits of precision
+        sage: QQbar(RDF(1))
+        Traceback (most recent call last):
+        ...
+        TypeError: Illegal initializer for algebraic number
+        sage: AA(RR(1))
+        Traceback (most recent call last):
+        ...
+        TypeError: Illegal initializer for algebraic number
+
+    This causes the following error when we try to scale a matrix of
+    complex numbers by an inexact real number::
+
+        sage: ComplexHermitianEJA(2,field=RR)
+        Traceback (most recent call last):
+        ...
+        TypeError: Unable to coerce entries (=(1.00000000000000,
+        -0.000000000000000)) to coefficients in Algebraic Real Field

 
     TESTS:
 
 
     TESTS:
 


@@ -2210,92 +2060,16 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
 
         sage: ComplexHermitianEJA(0)
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
         sage: ComplexHermitianEJA(0)
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field

-

     """
     """

-
-    @classmethod
-    def _denormalized_basis(cls, n):
-        """
-        Returns a basis for the space of complex Hermitian n-by-n matrices.
-
-        Why do we embed these? Basically, because all of numerical linear
-        algebra assumes that you're working with vectors consisting of `n`
-        entries from a field and scalars from the same field. There's no way
-        to tell SageMath that (for example) the vectors contain complex
-        numbers, while the scalar field is real.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
-
-        TESTS::
-
-            sage: set_random_seed()
-            sage: n = ZZ.random_element(1,5)
-            sage: B = ComplexHermitianEJA._denormalized_basis(n)
-            sage: all( M.is_symmetric() for M in  B)
-            True
-
-        """
-        field = ZZ
-        R = PolynomialRing(field, 'z')
-        z = R.gen()
-        F = field.extension(z**2 + 1, 'I')
-        I = F.gen(1)
-
-        # This is like the symmetric case, but we need to be careful:
-        #
-        #   * We want conjugate-symmetry, not just symmetry.
-        #   * The diagonal will (as a result) be real.
-        #
-        S = []
-        Eij = matrix.zero(F,n)
-        for i in range(n):
-            for j in range(i+1):
-                # "build" E_ij
-                Eij[i,j] = 1
-                if i == j:
-                    Sij = cls.real_embed(Eij)
-                    S.append(Sij)
-                else:
-                    # The second one has a minus because it's conjugated.
-                    Eij[j,i] = 1 # Eij = Eij + Eij.transpose()
-                    Sij_real = cls.real_embed(Eij)
-                    S.append(Sij_real)
-                    # Eij = I*Eij - I*Eij.transpose()
-                    Eij[i,j] = I
-                    Eij[j,i] = -I
-                    Sij_imag = cls.real_embed(Eij)
-                    S.append(Sij_imag)
-                    Eij[j,i] = 0
-                # "erase" E_ij
-                Eij[i,j] = 0
-
-        # Since we embedded these, we can drop back to the "field" that we
-        # started with instead of the complex extension "F".
-        return tuple( s.change_ring(field) for s in S )
-
-
-    def __init__(self, n, **kwargs):
+    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
 

-        associative = False
-        if n <= 1:
-            associative = True
+        from mjo.hurwitz import ComplexMatrixAlgebra
+        A = ComplexMatrixAlgebra(n, scalars=field)
+        super().__init__(A, **kwargs)

 
 

-        super().__init__(self._denormalized_basis(n),
-                         self.jordan_product,
-                         self.trace_inner_product,
-                         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)
-        self.one.set_cache(self(idV))

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


@@ -2309,155 +2083,8 @@ 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):
-
-    # A manual dictionary-cache for the quaternion_extension() method,
-    # since apparently @classmethods can't also be @cached_methods.
-    _quaternion_extension = {}
-
-    @classmethod
-    def quaternion_extension(cls,field):
-        r"""
-        The quaternion field that we embed/unembed, as an extension
-        of the given ``field``.
-        """
-        if field in cls._quaternion_extension:
-            return cls._quaternion_extension[field]
-
-        Q = QuaternionAlgebra(field,-1,-1)
-
-        cls._quaternion_extension[field] = Q
-        return Q
-
-    @staticmethod
-    def dimension_over_reals():
-        return 4
-
-    @classmethod
-    def real_embed(cls,M):
-        """
-        Embed the n-by-n quaternion matrix ``M`` into the space of real
-        matrices of size 4n-by-4n by first sending each quaternion entry `z
-        = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
-        c+di],[-c + di, a-bi]]`, and then embedding those into a real
-        matrix.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
-
-        EXAMPLES::
-
-            sage: Q = QuaternionAlgebra(QQ,-1,-1)
-            sage: i,j,k = Q.gens()
-            sage: x = 1 + 2*i + 3*j + 4*k
-            sage: M = matrix(Q, 1, [[x]])
-            sage: QuaternionMatrixEJA.real_embed(M)
-            [ 1  2  3  4]
-            [-2  1 -4  3]
-            [-3  4  1 -2]
-            [-4 -3  2  1]
-
-        Embedding is a homomorphism (isomorphism, in fact)::
-
-            sage: set_random_seed()
-            sage: n = ZZ.random_element(2)
-            sage: Q = QuaternionAlgebra(QQ,-1,-1)
-            sage: X = random_matrix(Q, n)
-            sage: Y = random_matrix(Q, n)
-            sage: Xe = QuaternionMatrixEJA.real_embed(X)
-            sage: Ye = QuaternionMatrixEJA.real_embed(Y)
-            sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
-            sage: Xe*Ye == XYe
-            True
-
-        """
-        super().real_embed(M)
-        quaternions = M.base_ring()
-        n = M.nrows()
-
-        F = QuadraticField(-1, 'I')
-        i = F.gen()
-
-        blocks = []
-        for z in M.list():
-            t = z.coefficient_tuple()
-            a = t[0]
-            b = t[1]
-            c = t[2]
-            d = t[3]
-            cplxM = matrix(F, 2, [[ a + b*i, c + d*i],
-                                 [-c + d*i, a - b*i]])
-            realM = ComplexMatrixEJA.real_embed(cplxM)
-            blocks.append(realM)
-
-        # We should have real entries by now, so use the realest field
-        # we've got for the return value.
-        return matrix.block(quaternions.base_ring(), n, blocks)
-
-
-
-    @classmethod
-    def real_unembed(cls,M):
-        """
-        The inverse of _embed_quaternion_matrix().
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
-
-        EXAMPLES::
-
-            sage: M = matrix(QQ, [[ 1,  2,  3,  4],
-            ....:                 [-2,  1, -4,  3],
-            ....:                 [-3,  4,  1, -2],
-            ....:                 [-4, -3,  2,  1]])
-            sage: QuaternionMatrixEJA.real_unembed(M)
-            [1 + 2*i + 3*j + 4*k]
-
-        TESTS:
-
-        Unembedding is the inverse of embedding::
-
-            sage: set_random_seed()
-            sage: Q = QuaternionAlgebra(QQ, -1, -1)
-            sage: M = random_matrix(Q, 3)
-            sage: Me = QuaternionMatrixEJA.real_embed(M)
-            sage: QuaternionMatrixEJA.real_unembed(Me) == M
-            True

 
 

-        """
-        super().real_unembed(M)
-        n = ZZ(M.nrows())
-        d = cls.dimension_over_reals()
-
-        # Use the base ring of the matrix to ensure that its entries can be
-        # multiplied by elements of the quaternion algebra.
-        Q = cls.quaternion_extension(M.base_ring())
-        i,j,k = Q.gens()
-
-        # Go top-left to bottom-right (reading order), converting every
-        # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
-        # quaternion block.
-        elements = []
-        for l in range(n/d):
-            for m in range(n/d):
-                submat = ComplexMatrixEJA.real_unembed(
-                    M[d*l:d*l+d,d*m:d*m+d] )
-                if submat[0,0] != submat[1,1].conjugate():
-                    raise ValueError('bad on-diagonal submatrix')
-                if submat[0,1] != -submat[1,0].conjugate():
-                    raise ValueError('bad off-diagonal submatrix')
-                z  = submat[0,0].real()
-                z += submat[0,0].imag()*i
-                z += submat[0,1].real()*j
-                z += submat[0,1].imag()*k
-                elements.append(z)
-
-        return matrix(Q, n/d, elements)
-
-
-class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
+class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):

     r"""
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the
     r"""
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the


@@ -2514,108 +2141,124 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """

-    @classmethod
-    def _denormalized_basis(cls, n):
-        """
-        Returns a basis for the space of quaternion Hermitian n-by-n matrices.
-
-        Why do we embed these? Basically, because all of numerical
-        linear algebra assumes that you're working with vectors consisting
-        of `n` entries from a field and scalars from the same field. There's
-        no way to tell SageMath that (for example) the vectors contain
-        complex numbers, while the scalar field is real.
-
-        SETUP::
+    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

 
 

-            sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
+        from mjo.hurwitz import QuaternionMatrixAlgebra
+        A = QuaternionMatrixAlgebra(n, scalars=field)
+        super().__init__(A, **kwargs)

 
 

-        TESTS::

 
 

-            sage: set_random_seed()
-            sage: n = ZZ.random_element(1,5)
-            sage: B = QuaternionHermitianEJA._denormalized_basis(n)
-            sage: all( M.is_symmetric() for M in B )
-            True
+    @staticmethod
+    def _max_random_instance_size():
+        r"""
+        The maximum rank of a random QuaternionHermitianEJA.
+        """
+        return 2 # Dimension 6

 
 

+    @classmethod
+    def random_instance(cls, **kwargs):
+        """
+        Return a random instance of this type of algebra.

         """
         """

-        field = ZZ
-        Q = QuaternionAlgebra(QQ,-1,-1)
-        I,J,K = Q.gens()
+        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        return cls(n, **kwargs)

 
 

-        # This is like the symmetric case, but we need to be careful:
-        #
-        #   * We want conjugate-symmetry, not just symmetry.
-        #   * The diagonal will (as a result) be real.
-        #
-        S = []
-        Eij = matrix.zero(Q,n)
-        for i in range(n):
-            for j in range(i+1):
-                # "build" E_ij
-                Eij[i,j] = 1
-                if i == j:
-                    Sij = cls.real_embed(Eij)
-                    S.append(Sij)
-                else:
-                    # The second, third, and fourth ones have a minus
-                    # because they're conjugated.
-                    # Eij = Eij + Eij.transpose()
-                    Eij[j,i] = 1
-                    Sij_real = cls.real_embed(Eij)
-                    S.append(Sij_real)
-                    # Eij = I*(Eij - Eij.transpose())
-                    Eij[i,j] = I
-                    Eij[j,i] = -I
-                    Sij_I = cls.real_embed(Eij)
-                    S.append(Sij_I)
-                    # Eij = J*(Eij - Eij.transpose())
-                    Eij[i,j] = J
-                    Eij[j,i] = -J
-                    Sij_J = cls.real_embed(Eij)
-                    S.append(Sij_J)
-                    # Eij = K*(Eij - Eij.transpose())
-                    Eij[i,j] = K
-                    Eij[j,i] = -K
-                    Sij_K = cls.real_embed(Eij)
-                    S.append(Sij_K)
-                    Eij[j,i] = 0
-                # "erase" E_ij
-                Eij[i,j] = 0
-
-        # Since we embedded these, we can drop back to the "field" that we
-        # started with instead of the quaternion algebra "Q".
-        return tuple( s.change_ring(field) for s in S )
-
-
-    def __init__(self, n, **kwargs):
-        # We know this is a valid EJA, but will double-check
-        # if the user passes check_axioms=True.
-        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
+class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
+    r"""
+    SETUP::

 
 

-        associative = False
-        if n <= 1:
-            associative = True
+        sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
+        ....:                                  OctonionHermitianEJA)
+        sage: from mjo.hurwitz import Octonions, OctonionMatrixAlgebra

 
 

-        super().__init__(self._denormalized_basis(n),
-                         self.jordan_product,
-                         self.trace_inner_product,
-                         associative=associative,
-                         **kwargs)
+    EXAMPLES:

 
 

-        # TODO: this could be factored out somehow, but is left here
-        # because the MatrixEJA is not presently a subclass of the
-        # FDEJA class that defines rank() and one().
-        self.rank.set_cache(n)
-        idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
-        self.one.set_cache(self(idV))
+    The 3-by-3 algebra satisfies the axioms of an EJA::
+
+        sage: OctonionHermitianEJA(3,                    # long time
+        ....:                      field=QQ,             # long time
+        ....:                      orthonormalize=False, # long time
+        ....:                      check_axioms=True)    # long time
+        Euclidean Jordan algebra of dimension 27 over Rational Field
+
+    After a change-of-basis, the 2-by-2 algebra has the same
+    multiplication table as the ten-dimensional Jordan spin algebra::
+
+        sage: A = OctonionMatrixAlgebra(2,Octonions(QQ),QQ)
+        sage: b = OctonionHermitianEJA._denormalized_basis(A)
+        sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],)
+        sage: jp = OctonionHermitianEJA.jordan_product
+        sage: ip = OctonionHermitianEJA.trace_inner_product
+        sage: J = 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

 
 

+    """

     @staticmethod
     def _max_random_instance_size():
         r"""
         The maximum rank of a random QuaternionHermitianEJA.
         """
     @staticmethod
     def _max_random_instance_size():
         r"""
         The maximum rank of a random QuaternionHermitianEJA.
         """

-        return 2 # Dimension 6
+        return 1 # Dimension 1

 
     @classmethod
     def random_instance(cls, **kwargs):
 
     @classmethod
     def random_instance(cls, **kwargs):


@@ -2625,15 +2268,58 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
         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, 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
+
+        from mjo.hurwitz import OctonionMatrixAlgebra
+        A = OctonionMatrixAlgebra(n, scalars=field)
+        super().__init__(A, **kwargs)
+
+        if n == 3:
+            from mjo.eja.eja_cache import albert_eja_coeffs
+            a = albert_eja_coeffs(self.coordinate_polynomial_ring())
+            if self._rational_algebra is None:
+                self._charpoly_coefficients.set_cache(a)
+            else:
+                self._rational_algebra._charpoly_coefficients.set_cache(a)
+
+
+class AlbertEJA(OctonionHermitianEJA):
+    r"""
+    The Albert algebra is the algebra of three-by-three Hermitian
+    matrices whose entries are octonions.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import AlbertEJA
+
+    EXAMPLES::
+
+        sage: AlbertEJA(field=QQ, orthonormalize=False)
+        Euclidean Jordan algebra of dimension 27 over Rational Field
+        sage: AlbertEJA() # long time
+        Euclidean Jordan algebra of dimension 27 over Algebraic Real Field

 
 

-class HadamardEJA(ConcreteEJA):

     """
     """

-    Return the Euclidean Jordan Algebra corresponding to the set
-    `R^n` under the Hadamard product.
+    def __init__(self, *args, **kwargs):
+        super().__init__(3, *args, **kwargs)

 
 

-    Note: this is nothing more than the Cartesian product of ``n``
-    copies of the spin algebra. Once Cartesian product algebras
-    are implemented, this can go.
+
+class HadamardEJA(RationalBasisEJA, ConcreteEJA):
+    """
+    Return the Euclidean Jordan algebra on `R^n` with the Hadamard
+    (pointwise real-number multiplication) Jordan product and the
+    usual inner-product.
+
+    This is nothing more than the Cartesian product of ``n`` copies of
+    the one-dimensional Jordan spin algebra, and is the most common
+    example of a non-simple Euclidean Jordan algebra.

 
     SETUP::
 
 
     SETUP::
 


@@ -2644,19 +2330,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:
 


@@ -2664,16 +2350,16 @@ class HadamardEJA(ConcreteEJA):
 
         sage: HadamardEJA(3, prefix='r').gens()
         (r0, r1, r2)
 
         sage: HadamardEJA(3, prefix='r').gens()
         (r0, r1, r2)

-

     """
     """

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

         if n == 0:
             jordan_product = lambda x,y: x
             inner_product = lambda x,y: x
         else:
             def jordan_product(x,y):
         if n == 0:
             jordan_product = lambda x,y: x
             inner_product = lambda x,y: x
         else:
             def jordan_product(x,y):

-                P = x.parent()
-                return P( xi*yi for (xi,yi) in zip(x,y) )
+                return MS( xi*yi for (xi,yi) in zip(x,y) )

 
             def inner_product(x,y):
                 return (x.T*y)[0,0]
 
             def inner_product(x,y):
                 return (x.T*y)[0,0]


@@ -2687,18 +2373,17 @@ class HadamardEJA(ConcreteEJA):
         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 

-        column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )
+        column_basis = tuple( MS(b) for b in FreeModule(field, n).basis() )

         super().__init__(column_basis,
                          jordan_product,
                          inner_product,
         super().__init__(column_basis,
                          jordan_product,
                          inner_product,

+                         field=field,
+                         matrix_space=MS,

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

-        if n == 0:
-            self.one.set_cache( self.zero() )
-        else:
-            self.one.set_cache( sum(self.gens()) )
+        self.one.set_cache( self.sum(self.gens()) )

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


@@ -2716,7 +2401,7 @@ class HadamardEJA(ConcreteEJA):
         return cls(n, **kwargs)
 
 
         return cls(n, **kwargs)
 
 

-class BilinearFormEJA(ConcreteEJA):
+class BilinearFormEJA(RationalBasisEJA, ConcreteEJA):

     r"""
     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
     with the half-trace inner product and jordan product ``x*y =
     r"""
     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
     with the half-trace inner product and jordan product ``x*y =


@@ -2798,7 +2483,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...


@@ -2812,21 +2497,22 @@ class BilinearFormEJA(ConcreteEJA):
         # verify things, we'll skip the rest of the checks.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
         # verify things, we'll skip the rest of the checks.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 

+        n = B.nrows()
+        MS = MatrixSpace(field, n, 1)
+

         def inner_product(x,y):
             return (y.T*B*x)[0,0]
 
         def jordan_product(x,y):
         def inner_product(x,y):
             return (y.T*B*x)[0,0]
 
         def jordan_product(x,y):

-            P = x.parent()

             x0 = x[0,0]
             xbar = x[1:,0]
             y0 = y[0,0]
             ybar = y[1:,0]
             z0 = inner_product(y,x)
             zbar = y0*xbar + x0*ybar
             x0 = x[0,0]
             xbar = x[1:,0]
             y0 = y[0,0]
             ybar = y[1:,0]
             z0 = inner_product(y,x)
             zbar = y0*xbar + x0*ybar

-            return P([z0] + zbar.list())
+            return MS([z0] + zbar.list())

 
 

-        n = B.nrows()
-        column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )
+        column_basis = tuple( MS(b) for b in FreeModule(field, n).basis() )

 
         # TODO: I haven't actually checked this, but it seems legit.
         associative = False
 
         # TODO: I haven't actually checked this, but it seems legit.
         associative = False


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

+                         field=field,
+                         matrix_space=MS,

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


@@ -2843,7 +2531,6 @@ class BilinearFormEJA(ConcreteEJA):
         # one-dimensional ambient space (because the rank is bounded
         # by the ambient dimension).
         self.rank.set_cache(min(n,2))
         # one-dimensional ambient space (because the rank is bounded
         # by the ambient dimension).
         self.rank.set_cache(min(n,2))

-

         if n == 0:
             self.one.set_cache( self.zero() )
         else:
         if n == 0:
             self.one.set_cache( self.zero() )
         else:


@@ -2897,20 +2584,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 +2618,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 +2629,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():


@@ -2962,7 +2649,7 @@ class JordanSpinEJA(BilinearFormEJA):
         return cls(n, **kwargs)
 
 
         return cls(n, **kwargs)
 
 

-class TrivialEJA(ConcreteEJA):
+class TrivialEJA(RationalBasisEJA, ConcreteEJA):

     """
     The trivial Euclidean Jordan algebra consisting of only a zero element.
 
     """
     The trivial Euclidean Jordan algebra consisting of only a zero element.
 


@@ -2991,10 +2678,11 @@ class TrivialEJA(ConcreteEJA):
         0
 
     """
         0
 
     """

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

         jordan_product = lambda x,y: x
         jordan_product = lambda x,y: x

-        inner_product = lambda x,y: 0
+        inner_product = lambda x,y: field.zero()

         basis = ()
         basis = ()

+        MS = MatrixSpace(field,0)

 
         # New defaults for keyword arguments
         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
 
         # New defaults for keyword arguments
         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False


@@ -3004,6 +2692,8 @@ class TrivialEJA(ConcreteEJA):
                          jordan_product,
                          inner_product,
                          associative=True,
                          jordan_product,
                          inner_product,
                          associative=True,

+                         field=field,
+                         matrix_space=MS,

                          **kwargs)
 
         # The rank is zero using my definition, namely the dimension of the
                          **kwargs)
 
         # The rank is zero using my definition, namely the dimension of the


@@ -3018,8 +2708,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


@@ -3122,24 +2811,24 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
         sage: J3 = JordanSpinEJA(1)
         sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
         sage: J.multiplication_table()
         sage: J3 = JordanSpinEJA(1)
         sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
         sage: J.multiplication_table()

-        +--------------++---------+--------------+--------------+
-        | *            || e(0, 0) | e(1, (0, 0)) | e(1, (1, 0)) |
-        +==============++=========+==============+==============+
-        | e(0, 0)      || e(0, 0) | 0            | 0            |
-        +--------------++---------+--------------+--------------+
-        | e(1, (0, 0)) || 0       | e(1, (0, 0)) | 0            |
-        +--------------++---------+--------------+--------------+
-        | e(1, (1, 0)) || 0       | 0            | e(1, (1, 0)) |
-        +--------------++---------+--------------+--------------+
+        +----++----+----+----+
+        | *  || b0 | b1 | b2 |
+        +====++====+====+====+
+        | b0 || b0 | 0  | 0  |
+        +----++----+----+----+
+        | b1 || 0  | b1 | 0  |
+        +----++----+----+----+
+        | b2 || 0  | 0  | b2 |
+        +----++----+----+----+

         sage: HadamardEJA(3).multiplication_table()
         +----++----+----+----+
         sage: HadamardEJA(3).multiplication_table()
         +----++----+----+----+

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

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

-        | e0 || e0 | 0  | 0  |
+        | b0 || b0 | 0  | 0  |

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

-        | e1 || 0  | e1 | 0  |
+        | b1 || 0  | b1 | 0  |

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

-        | e2 || 0  | 0  | e2 |
+        | b2 || 0  | 0  | b2 |

         +----++----+----+----+
 
     TESTS:
         +----++----+----+----+
 
     TESTS:


@@ -3169,37 +2858,48 @@ 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 )
-
-        # 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()
-        )
+        associative = all( f.is_associative() for f in factors )
+
+        # Compute my matrix space. This category isn't perfect, but
+        # is good enough for what we need to do.
+        MS_cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
+        MS_cat = MS_cat.Unital().CartesianProducts()
+        MS_factors = tuple( J.matrix_space() for J in factors )
+        from sage.sets.cartesian_product import CartesianProduct
+        MS = CartesianProduct(MS_factors, MS_cat)
+
+        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


@@ -3213,106 +2913,49 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
                                       jordan_product,
                                       inner_product,
                                       field=field,
                                       jordan_product,
                                       inner_product,
                                       field=field,

+                                      matrix_space=MS,

                                       orthonormalize=False,
                                       associative=associative,
                                       cartesian_product=True,
                                       check_field=False,
                                       check_axioms=False)
 
                                       orthonormalize=False,
                                       associative=associative,
                                       cartesian_product=True,
                                       check_field=False,
                                       check_axioms=False)
 

-        ones = tuple(J.one() for J in algebras)
-        self.one.set_cache(self._cartesian_product_of_elements(ones))
-        self.rank.set_cache(sum(J.rank() for J in algebras))
-
-    def _monomial_to_generator(self, mon):
-        r"""
-        Convert a monomial index into a generator index.
-
-        This is needed in product algebras because the multiplication
-        table is in terms of the generator indices.
-
-        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() )
-            True
-
-        """
-        # This works recursively so that we can handle Cartesian
-        # products of Cartesian products.
-        try:
-            # monomial is an ordered pair
-            factor = mon[0]
-        except TypeError: # 'int' object is not subscriptable
-            # base case where the monomials are integers
-            return mon
-
-        idx_in_factor = self._monomial_to_generator(mon[1])
+        self.rank.set_cache(sum(J.rank() for J in factors))
+        ones = tuple(J.one().to_matrix() for J in factors)
+        self.one.set_cache(self(ones))

 
 

-        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"""
         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: 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.
 
 
     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::
 
         SETUP::
 

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

             ....:                                  RealSymmetricEJA)
 
         EXAMPLES::
             ....:                                  RealSymmetricEJA)
 
         EXAMPLES::


@@ -3325,10 +2968,37 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
             matrices over Algebraic Real Field, Full MatrixSpace of 2
             by 2 dense matrices over Algebraic Real Field)
 
             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 |
+            +---+
+            sage: J.one().to_matrix()[1]
+            +---+
+            | 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()] )
+        return super().matrix_space()
+

 
     @cached_method
     def cartesian_projection(self, i):
 
     @cached_method
     def cartesian_projection(self, i):


@@ -3400,9 +3070,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


@@ -3508,9 +3181,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())
 
 


@@ -3525,7 +3200,9 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
 
     SETUP::
 
 
     SETUP::
 

-        sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+        sage: from mjo.eja.eja_algebra import (HadamardEJA,
+        ....:                                  JordanSpinEJA,
+        ....:                                  OctonionHermitianEJA,

         ....:                                  RealSymmetricEJA)
 
     EXAMPLES:
         ....:                                  RealSymmetricEJA)
 
     EXAMPLES:


@@ -3542,30 +3219,45 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
         sage: J.rank()
         5
 
         sage: J.rank()
         5
 

+    TESTS:
+
+    The ``cartesian_product()`` function only uses the first factor to
+    decide where the result will live; thus we have to be careful to
+    check that all factors do indeed have a `_rational_algebra` 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:
     """
     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
 
 
 
 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