]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: don't allow creation of invalid RationalBasisEJAs.
[sage.d.git] / mjo / eja / eja_algebra.py
index 8b37a83602ef9b00b5a14c55785c2163b44df32b..af4080b0807d6ac6848849f23edd237657896413 100644 (file)
@@ -170,6 +170,17 @@ from mjo.eja.eja_element import FiniteDimensionalEJAElement
 from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
 from mjo.eja.eja_utils import _all2list, _mat2vec
 
+def EuclideanJordanAlgebras(field):
+    r"""
+    The category of Euclidean Jordan algebras over ``field``, which
+    must be a subfield of the real numbers. For now this is just a
+    convenient wrapper around all of the other category axioms that
+    apply to all EJAs.
+    """
+    category = MagmaticAlgebras(field).FiniteDimensional()
+    category = category.WithBasis().Unital().Commutative()
+    return category
+
 class FiniteDimensionalEJA(CombinatorialFreeModule):
     r"""
     A finite-dimensional Euclidean Jordan algebra.
@@ -228,6 +239,26 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
     """
     Element = FiniteDimensionalEJAElement
 
+    @staticmethod
+    def _check_input_field(field):
+        if not field.is_subring(RR):
+            # Note: this does return true for the real algebraic
+            # field, the rationals, and any quadratic field where
+            # we've specified a real embedding.
+            raise ValueError("scalar field is not real")
+
+    @staticmethod
+    def _check_input_axioms(basis, jordan_product, inner_product):
+        if not all( jordan_product(bi,bj) == jordan_product(bj,bi)
+                    for bi in basis
+                    for bj in basis ):
+            raise ValueError("Jordan product is not commutative")
+
+        if not all( inner_product(bi,bj) == inner_product(bj,bi)
+                    for bi in basis
+                    for bj in basis ):
+            raise ValueError("inner-product is not commutative")
+
     def __init__(self,
                  basis,
                  jordan_product,
@@ -236,7 +267,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                  matrix_space=None,
                  orthonormalize=True,
                  associative=None,
-                 cartesian_product=False,
                  check_field=True,
                  check_axioms=True,
                  prefix="b"):
@@ -244,30 +274,14 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         n = len(basis)
 
         if check_field:
-            if not field.is_subring(RR):
-                # Note: this does return true for the real algebraic
-                # field, the rationals, and any quadratic field where
-                # we've specified a real embedding.
-                raise ValueError("scalar field is not real")
+            self._check_input_field(field)
 
         if check_axioms:
             # Check commutativity of the Jordan and inner-products.
             # This has to be done before we build the multiplication
             # and inner-product tables/matrices, because we take
             # advantage of symmetry in the process.
-            if not all( jordan_product(bi,bj) == jordan_product(bj,bi)
-                        for bi in basis
-                        for bj in basis ):
-                raise ValueError("Jordan product is not commutative")
-
-            if not all( inner_product(bi,bj) == inner_product(bj,bi)
-                        for bi in basis
-                        for bj in basis ):
-                raise ValueError("inner-product is not commutative")
-
-
-        category = MagmaticAlgebras(field).FiniteDimensional()
-        category = category.WithBasis().Unital().Commutative()
+            self._check_input_axioms(basis, jordan_product, inner_product)
 
         if n <= 1:
             # All zero- and one-dimensional algebras are just the real
@@ -286,14 +300,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                                for bj in basis
                                for bk in basis)
 
+        category = EuclideanJordanAlgebras(field)
+
         if associative:
             # Element subalgebras can take advantage of this.
             category = category.Associative()
-        if cartesian_product:
-            # 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.
@@ -309,7 +320,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # 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*b1 + 2*b2.
-        vector_basis = basis
 
         degree = 0
         if n > 0:
@@ -319,9 +329,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # written out as "long vectors."
         V = VectorSpace(field, degree)
 
-        # The matrix that will hole the orthonormal -> unorthonormal
-        # coordinate transformation.
-        self._deortho_matrix = None
+        # The matrix that will hold the orthonormal -> unorthonormal
+        # coordinate transformation. Default to an identity matrix of
+        # the appropriate size to avoid special cases for None
+        # everywhere.
+        self._deortho_matrix = matrix.identity(field,n)
 
         if orthonormalize:
             # Save a copy of the un-orthonormalized basis for later.
@@ -346,23 +358,29 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # its own set of non-ambient coordinates (in terms of the
         # supplied basis).
         vector_basis = tuple( V(_all2list(b)) for b in basis )
-        W = V.span_of_basis( vector_basis, check=check_axioms)
+
+        # Save the span of our matrix basis (when written out as long
+        # vectors) because otherwise we'll have to reconstruct it
+        # every time we want to coerce a matrix into the algebra.
+        self._matrix_span = V.span_of_basis( vector_basis, check=check_axioms)
 
         if orthonormalize:
-            # Now "W" is the vector space of our algebra coordinates. The
-            # variables "X1", "X2",...  refer to the entries of vectors in
-            # W. Thus to convert back and forth between the orthonormal
-            # coordinates and the given ones, we need to stick the original
-            # basis in W.
+            # Now "self._matrix_span" is the vector space of our
+            # algebra coordinates. The variables "X1", "X2",...  refer
+            # to the entries of vectors in self._matrix_span. Thus to
+            # convert back and forth between the orthonormal
+            # coordinates and the given ones, we need to stick the
+            # original basis in self._matrix_span.
             U = V.span_of_basis( deortho_vector_basis, check=check_axioms)
-            self._deortho_matrix = matrix( U.coordinate_vector(q)
-                                           for q in vector_basis )
+            self._deortho_matrix = matrix.column( U.coordinate_vector(q)
+                                                  for q in vector_basis )
 
 
         # Now we actually compute the multiplication and inner-product
         # tables/matrices using the possibly-orthonormalized basis.
         self._inner_product_matrix = matrix.identity(field, n)
-        self._multiplication_table = [ [0 for j in range(i+1)]
+        zed = self.zero()
+        self._multiplication_table = [ [zed for j in range(i+1)]
                                        for i in range(n) ]
 
         # Note: the Jordan and inner-products are defined in terms
@@ -377,7 +395,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                 # The jordan product returns a matrixy answer, so we
                 # have to convert it to the algebra coordinates.
                 elt = jordan_product(q_i, q_j)
-                elt = W.coordinate_vector(V(_all2list(elt)))
+                elt = self._matrix_span.coordinate_vector(V(_all2list(elt)))
                 self._multiplication_table[i][j] = self.from_vector(elt)
 
                 if not orthonormalize:
@@ -684,8 +702,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
     def _element_constructor_(self, elt):
         """
-        Construct an element of this algebra from its vector or matrix
-        representation.
+        Construct an element of this algebra or a subalgebra from its
+        EJA element, vector, or matrix representation.
 
         This gets called only after the parent element _call_ method
         fails to find a coercion for the argument.
@@ -724,6 +742,16 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
             b1 + b5
 
+        Subalgebra elements are embedded into the superalgebra::
+
+            sage: J = JordanSpinEJA(3)
+            sage: J.one()
+            b0
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by()
+            sage: J(A.one())
+            b0
+
         TESTS:
 
         Ensure that we can convert any element back and forth
@@ -748,6 +776,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             Traceback (most recent call last):
             ...
             ValueError: not an element of this algebra
+
         """
         msg = "not an element of this algebra"
         if elt in self.base_ring():
@@ -757,13 +786,16 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             # that the integer 3 belongs to the space of 2-by-2 matrices.
             raise ValueError(msg)
 
-        try:
-            # Try to convert a vector into a column-matrix...
+        if hasattr(elt, 'superalgebra_element'):
+            # Handle subalgebra elements
+            if elt.parent().superalgebra() == self:
+                return elt.superalgebra_element()
+
+        if hasattr(elt, 'sparse_vector'):
+            # Convert a vector into a column-matrix. We check for
+            # "sparse_vector" and not "column" because matrices also
+            # have a "column" method.
             elt = elt.column()
-        except (AttributeError, TypeError):
-            # and ignore failure, because we weren't really expecting
-            # a vector as an argument anyway.
-            pass
 
         if elt not in self.matrix_space():
             raise ValueError(msg)
@@ -780,15 +812,10 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # is that we're already converting everything to long vectors,
         # and that strategy works for tuples as well.
         #
-        # We pass check=False because the matrix basis is "guaranteed"
-        # to be linearly independent... right? Ha ha.
-        elt = _all2list(elt)
-        V = VectorSpace(self.base_ring(), len(elt))
-        W = V.span_of_basis( (V(_all2list(s)) for s in self.matrix_basis()),
-                             check=False)
+        elt = self._matrix_span.ambient_vector_space()(_all2list(elt))
 
         try:
-            coords = W.coordinate_vector(V(elt))
+            coords = self._matrix_span.coordinate_vector(elt)
         except ArithmeticError:  # vector is not in free module
             raise ValueError(msg)
 
@@ -1394,7 +1421,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # corresponding to trivial spaces (e.g. it returns only the
         # eigenspace corresponding to lambda=1 if you take the
         # decomposition relative to the identity element).
-        trivial = self.subalgebra(())
+        trivial = self.subalgebra((), check_axioms=False)
         J0 = trivial                          # eigenvalue zero
         J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
         J1 = trivial                          # eigenvalue one
@@ -1716,6 +1743,15 @@ class RationalBasisEJA(FiniteDimensionalEJA):
                                        check_field=False,
                                        check_axioms=False)
 
+    def rational_algebra(self):
+        # Using None as a flag here (rather than just assigning "self"
+        # to self._rational_algebra by default) feels a little bit
+        # more sane to me in a garbage-collected environment.
+        if self._rational_algebra is None:
+            return self
+        else:
+            return self._rational_algebra
+
     @cached_method
     def _charpoly_coefficients(self):
         r"""
@@ -1740,25 +1776,15 @@ class RationalBasisEJA(FiniteDimensionalEJA):
             Algebraic Real Field
 
         """
-        if self._rational_algebra is None:
-            # There's no need to construct *another* algebra over the
-            # rationals if this one is already over the
-            # rationals. Likewise, if we never orthonormalized our
-            # basis, we might as well just use the given one.
+        if self.rational_algebra() is self:
+            # Bypass the hijinks if they won't benefit us.
             return super()._charpoly_coefficients()
 
         # Do the computation over the rationals. The answer will be
         # the same, because all we've done is a change of basis.
         # Then, change back from QQ to our real base ring
         a = ( a_i.change_ring(self.base_ring())
-              for a_i in self._rational_algebra._charpoly_coefficients() )
-
-        if self._deortho_matrix is None:
-            # This can happen if our base ring was, say, AA and we
-            # chose not to (or didn't need to) orthonormalize. It's
-            # still faster to do the computations over QQ even if
-            # the numbers in the boxes stay the same.
-            return tuple(a)
+              for a_i in self.rational_algebra()._charpoly_coefficients() )
 
         # Otherwise, convert the coordinate variables back to the
         # deorthonormalized ones.
@@ -1833,11 +1859,36 @@ class ConcreteEJA(FiniteDimensionalEJA):
     def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this type of algebra whose dimension
-        is less than or equal to ``max_dimension``. If the dimension bound
-        is omitted, then the ``_max_random_instance_dimension()`` is used
-        to get a suitable bound.
+        is less than or equal to the lesser of ``max_dimension`` and
+        the value returned by ``_max_random_instance_dimension()``. If
+        the dimension bound is omitted, then only the
+        ``_max_random_instance_dimension()`` is used as a bound.
 
         This method should be implemented in each subclass.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import ConcreteEJA
+
+        TESTS:
+
+        Both the class bound and the ``max_dimension`` argument are upper
+        bounds on the dimension of the algebra returned::
+
+            sage: from sage.misc.prandom import choice
+            sage: eja_class = choice(ConcreteEJA.__subclasses__())
+            sage: class_max_d = eja_class._max_random_instance_dimension()
+            sage: J = eja_class.random_instance(max_dimension=20,
+            ....:                               field=QQ,
+            ....:                               orthonormalize=False)
+            sage: J.dimension() <= class_max_d
+            True
+            sage: J = eja_class.random_instance(max_dimension=2,
+            ....:                               field=QQ,
+            ....:                               orthonormalize=False)
+            sage: J.dimension() <= 2
+            True
+
         """
         from sage.misc.prandom import choice
         eja_class = choice(cls.__subclasses__())
@@ -1845,7 +1896,7 @@ class ConcreteEJA(FiniteDimensionalEJA):
         # These all bubble up to the RationalBasisEJA superclass
         # constructor, so any (kw)args valid there are also valid
         # here.
-        return eja_class.random_instance(*args, **kwargs)
+        return eja_class.random_instance(max_dimension, *args, **kwargs)
 
 
 class MatrixEJA(FiniteDimensionalEJA):
@@ -2077,14 +2128,15 @@ class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/2 - 1/2))
 
     @classmethod
-    def random_instance(cls, max_dimension=None, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this type of algebra.
         """
-        if max_dimension is None:
-            max_dimension = cls._max_random_instance_dimension()
-        max_size = cls._max_random_instance_size(max_dimension) + 1
-        n = ZZ.random_element(max_size)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
     def __init__(self, n, field=AA, **kwargs):
@@ -2098,10 +2150,7 @@ class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         from mjo.eja.eja_cache import real_symmetric_eja_coeffs
         a = real_symmetric_eja_coeffs(self)
         if a is not None:
-            if self._rational_algebra is None:
-                self._charpoly_coefficients.set_cache(a)
-            else:
-                self._rational_algebra._charpoly_coefficients.set_cache(a)
+            self.rational_algebra()._charpoly_coefficients.set_cache(a)
 
 
 
@@ -2189,10 +2238,7 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         from mjo.eja.eja_cache import complex_hermitian_eja_coeffs
         a = complex_hermitian_eja_coeffs(self)
         if a is not None:
-            if self._rational_algebra is None:
-                self._charpoly_coefficients.set_cache(a)
-            else:
-                self._rational_algebra._charpoly_coefficients.set_cache(a)
+            self.rational_algebra()._charpoly_coefficients.set_cache(a)
 
     @staticmethod
     def _max_random_instance_size(max_dimension):
@@ -2201,13 +2247,15 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         return ZZ(int(ZZ(max_dimension).sqrt()))
 
     @classmethod
-    def random_instance(cls, max_dimension=None, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this type of algebra.
         """
-        if max_dimension is None:
-            max_dimension = cls._max_random_instance_dimension()
-        n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
 
@@ -2280,10 +2328,7 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         from mjo.eja.eja_cache import quaternion_hermitian_eja_coeffs
         a = quaternion_hermitian_eja_coeffs(self)
         if a is not None:
-            if self._rational_algebra is None:
-                self._charpoly_coefficients.set_cache(a)
-            else:
-                self._rational_algebra._charpoly_coefficients.set_cache(a)
+            self.rational_algebra()._charpoly_coefficients.set_cache(a)
 
 
 
@@ -2297,13 +2342,15 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/4 + 1/4))
 
     @classmethod
-    def random_instance(cls, max_dimension=None, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this type of algebra.
         """
-        if max_dimension is None:
-            max_dimension = cls._max_random_instance_dimension()
-        n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
 class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
@@ -2410,13 +2457,15 @@ class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
             return 0
 
     @classmethod
-    def random_instance(cls, max_dimension=None, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this type of algebra.
         """
-        if max_dimension is None:
-            max_dimension = cls._max_random_instance_dimension()
-        n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
     def __init__(self, n, field=AA, **kwargs):
@@ -2435,10 +2484,7 @@ class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         from mjo.eja.eja_cache import octonion_hermitian_eja_coeffs
         a = octonion_hermitian_eja_coeffs(self)
         if a is not None:
-            if self._rational_algebra is None:
-                self._charpoly_coefficients.set_cache(a)
-            else:
-                self._rational_algebra._charpoly_coefficients.set_cache(a)
+            self.rational_algebra()._charpoly_coefficients.set_cache(a)
 
 
 class AlbertEJA(OctonionHermitianEJA):
@@ -2552,13 +2598,15 @@ class HadamardEJA(RationalBasisEJA, ConcreteEJA):
         return max_dimension
 
     @classmethod
-    def random_instance(cls, max_dimension=None, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this type of algebra.
         """
-        if max_dimension is None:
-            max_dimension = cls._max_random_instance_dimension()
-        n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
 
@@ -2713,13 +2761,16 @@ class BilinearFormEJA(RationalBasisEJA, ConcreteEJA):
         return max_dimension
 
     @classmethod
-    def random_instance(cls, max_dimension=None, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this algebra.
         """
-        if max_dimension is None:
-            max_dimension = cls._max_random_instance_dimension()
-        n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
+
         if n.is_zero():
             B = matrix.identity(ZZ, n)
             return cls(B, **kwargs)
@@ -2730,6 +2781,7 @@ class BilinearFormEJA(RationalBasisEJA, ConcreteEJA):
         alpha = ZZ.zero()
         while alpha.is_zero():
             alpha = ZZ.random_element().abs()
+
         B22 = M.transpose()*M + alpha*I
 
         from sage.matrix.special import block_matrix
@@ -2803,15 +2855,17 @@ class JordanSpinEJA(BilinearFormEJA):
         super().__init__(B, *args, **kwargs)
 
     @classmethod
-    def random_instance(cls, max_dimension=None, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this type of algebra.
 
         Needed here to override the implementation for ``BilinearFormEJA``.
         """
-        if max_dimension is None:
-            max_dimension = cls._max_random_instance_dimension()
-        n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
 
@@ -2868,9 +2922,12 @@ class TrivialEJA(RationalBasisEJA, ConcreteEJA):
         self.one.set_cache( self.zero() )
 
     @classmethod
-    def random_instance(cls, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         # We don't take a "size" argument so the superclass method is
-        # inappropriate for us.
+        # inappropriate for us. The ``max_dimension`` argument is
+        # included so that if this method is called generically with a
+        # ``max_dimension=<whatever>`` argument, we don't try to pass
+        # it on to the algebra constructor.
         return cls(**kwargs)
 
 
@@ -2886,6 +2943,7 @@ class CartesianProductEJA(FiniteDimensionalEJA):
 
         sage: from mjo.eja.eja_algebra import (random_eja,
         ....:                                  CartesianProductEJA,
+        ....:                                  ComplexHermitianEJA,
         ....:                                  HadamardEJA,
         ....:                                  JordanSpinEJA,
         ....:                                  RealSymmetricEJA)
@@ -2997,6 +3055,28 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         | b2 || 0  | 0  | b2 |
         +----++----+----+----+
 
+    The "matrix space" of a Cartesian product always consists of
+    ordered pairs (or triples, or...) whose components are the
+    matrix spaces of its factors::
+
+            sage: J1 = HadamardEJA(2)
+            sage: J2 = ComplexHermitianEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: J.matrix_space()
+            The Cartesian product of (Full MatrixSpace of 2 by 1 dense
+            matrices over Algebraic Real Field, Module of 2 by 2 matrices
+            with entries in Algebraic Field over the scalar ring Algebraic
+            Real Field)
+            sage: J.one().to_matrix()[0]
+            [1]
+            [1]
+            sage: J.one().to_matrix()[1]
+            +---+---+
+            | 1 | 0 |
+            +---+---+
+            | 0 | 1 |
+            +---+---+
+
     TESTS:
 
     All factors must share the same base field::
@@ -3019,11 +3099,7 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         sage: expected = J.one()             # long time
         sage: actual == expected             # long time
         True
-
     """
-    Element = FiniteDimensionalEJAElement
-
-
     def __init__(self, factors, **kwargs):
         m = len(factors)
         if m == 0:
@@ -3035,56 +3111,93 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         if not all( J.base_ring() == field for J in factors ):
             raise ValueError("all factors must share the same base field")
 
+        # Figure out the category to use.
         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.
+        category = EuclideanJordanAlgebras(field)
+        if associative: category = category.Associative()
+        category = category.join([category, category.CartesianProducts()])
+
+        # Compute my matrix space.  We don't simply use the
+        # ``cartesian_product()`` functor here because it acts
+        # differently on SageMath MatrixSpaces and our custom
+        # MatrixAlgebras, which are CombinatorialFreeModules. We
+        # always want the result to be represented (and indexed) as an
+        # ordered tuple. This category isn't perfect, but is good
+        # enough for what we need to do.
         MS_cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
         MS_cat = MS_cat.Unital().CartesianProducts()
         MS_factors = tuple( J.matrix_space() for J in factors )
         from sage.sets.cartesian_product import CartesianProduct
-        MS = CartesianProduct(MS_factors, MS_cat)
+        self._matrix_space = CartesianProduct(MS_factors, MS_cat)
 
-        basis = []
-        zero = MS.zero()
+        self._matrix_basis = []
+        zero = self._matrix_space.zero()
         for i in range(m):
             for b in factors[i].matrix_basis():
                 z = list(zero)
                 z[i] = b
-                basis.append(z)
+                self._matrix_basis.append(z)
 
-        basis = tuple( MS(b) for b in basis )
+        self._matrix_basis = tuple( self._matrix_space(b)
+                                    for b in self._matrix_basis )
+        n = len(self._matrix_basis)
 
-        # Define jordan/inner products that operate on that matrix_basis.
-        def jordan_product(x,y):
-            return MS(tuple(
-                (factors[i](x[i])*factors[i](y[i])).to_matrix()
-                for i in range(m)
-            ))
-
-        def inner_product(x, y):
-            return sum(
-                factors[i](x[i]).inner_product(factors[i](y[i]))
-                for i in range(m)
-            )
+        # We already have what we need for the super-superclass constructor.
+        CombinatorialFreeModule.__init__(self,
+                                         field,
+                                         range(n),
+                                         prefix="b",
+                                         category=category,
+                                         bracket=False)
 
-        # There's no need to check the field since it already came
-        # from an EJA. Likewise the axioms are guaranteed to be
-        # satisfied, unless the guy writing this class sucks.
-        #
-        # If you want the basis to be orthonormalized, orthonormalize
-        # the factors.
-        FiniteDimensionalEJA.__init__(self,
-                                      basis,
-                                      jordan_product,
-                                      inner_product,
-                                      field=field,
-                                      matrix_space=MS,
-                                      orthonormalize=False,
-                                      associative=associative,
-                                      cartesian_product=True,
-                                      check_field=False,
-                                      check_axioms=False)
+        # Now create the vector space for the algebra, which will have
+        # its own set of non-ambient coordinates (in terms of the
+        # supplied basis).
+        degree = sum( f._matrix_span.ambient_vector_space().degree()
+                      for f in factors )
+        V = VectorSpace(field, degree)
+        vector_basis = tuple( V(_all2list(b)) for b in self._matrix_basis )
+
+        # Save the span of our matrix basis (when written out as long
+        # vectors) because otherwise we'll have to reconstruct it
+        # every time we want to coerce a matrix into the algebra.
+        self._matrix_span = V.span_of_basis( vector_basis, check=False)
+
+        # Since we don't (re)orthonormalize the basis, the FDEJA
+        # constructor is going to set self._deortho_matrix to the
+        # identity matrix. Here we set it to the correct value using
+        # the deortho matrices from our factors.
+        self._deortho_matrix = matrix.block_diagonal(
+            [J._deortho_matrix for J in factors]
+        )
+
+        self._inner_product_matrix = matrix.block_diagonal(
+            [J._inner_product_matrix for J in factors]
+        )
+        self._inner_product_matrix._cache = {'hermitian': True}
+        self._inner_product_matrix.set_immutable()
+
+        # Building the multiplication table is a bit more tricky
+        # because we have to embed the entries of the factors'
+        # multiplication tables into the product EJA.
+        zed = self.zero()
+        self._multiplication_table = [ [zed for j in range(i+1)]
+                                       for i in range(n) ]
+
+        # Keep track of an offset that tallies the dimensions of all
+        # previous factors. If the second factor is dim=2 and if the
+        # first one is dim=3, then we want to skip the first 3x3 block
+        # when copying the multiplication table for the second factor.
+        offset = 0
+        for f in range(m):
+            phi_f = self.cartesian_embedding(f)
+            factor_dim = factors[f].dimension()
+            for i in range(factor_dim):
+                for j in range(i+1):
+                    f_ij = factors[f]._multiplication_table[i][j]
+                    e = phi_f(f_ij)
+                    self._multiplication_table[offset+i][offset+j] = e
+            offset += factor_dim
 
         self.rank.set_cache(sum(J.rank() for J in factors))
         ones = tuple(J.one().to_matrix() for J in factors)
@@ -3106,65 +3219,6 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         return cartesian_product.symbol.join("%s" % factor
                                              for factor in self._sets)
 
-    def matrix_space(self):
-        r"""
-        Return the space that our matrix basis lives in as a Cartesian
-        product.
-
-        We don't simply use the ``cartesian_product()`` functor here
-        because it acts differently on SageMath MatrixSpaces and our
-        custom MatrixAlgebras, which are CombinatorialFreeModules. We
-        always want the result to be represented (and indexed) as
-        an ordered tuple.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
-            ....:                                  HadamardEJA,
-            ....:                                  OctonionHermitianEJA,
-            ....:                                  RealSymmetricEJA)
-
-        EXAMPLES::
-
-            sage: J1 = HadamardEJA(1)
-            sage: J2 = RealSymmetricEJA(2)
-            sage: J = cartesian_product([J1,J2])
-            sage: J.matrix_space()
-            The Cartesian product of (Full MatrixSpace of 1 by 1 dense
-            matrices over Algebraic Real Field, Full MatrixSpace of 2
-            by 2 dense matrices over Algebraic Real Field)
-
-        ::
-
-            sage: 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 |
-            +----+
-
-        """
-        return super().matrix_space()
-
 
     @cached_method
     def cartesian_projection(self, i):
@@ -3366,9 +3420,9 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
 
     SETUP::
 
-        sage: from mjo.eja.eja_algebra import (HadamardEJA,
+        sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
+        ....:                                  HadamardEJA,
         ....:                                  JordanSpinEJA,
-        ....:                                  OctonionHermitianEJA,
         ....:                                  RealSymmetricEJA)
 
     EXAMPLES:
@@ -3389,33 +3443,58 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
 
     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
+    check that all factors do indeed have a ``rational_algebra()`` method
+    before we construct an algebra that claims to have a rational basis::
+
+        sage: J1 = HadamardEJA(2)
+        sage: jp = lambda X,Y: X*Y
+        sage: ip = lambda X,Y: X[0,0]*Y[0,0]
+        sage: b1 = matrix(QQ, [[1]])
+        sage: J2 = FiniteDimensionalEJA((b1,), jp, ip)
+        sage: cartesian_product([J2,J1]) # factor one not RationalBasisEJA
+        Euclidean Jordan algebra of dimension 1 over Algebraic Real
+        Field (+) Euclidean Jordan algebra of dimension 2 over Algebraic
+        Real Field
+        sage: cartesian_product([J1,J2]) # factor one is RationalBasisEJA
+        Traceback (most recent call last):
+        ...
+        ValueError: factor not a RationalBasisEJA
 
     """
     def __init__(self, algebras, **kwargs):
+        if not all( hasattr(r, "rational_algebra") for r in algebras ):
+            raise ValueError("factor not a RationalBasisEJA")
+
         CartesianProductEJA.__init__(self, algebras, **kwargs)
 
-        self._rational_algebra = None
-        if self.vector_space().base_field() is not QQ:
-            if all( hasattr(r, "_rational_algebra") for r in algebras ):
-                self._rational_algebra = cartesian_product([
-                    r._rational_algebra for r in algebras
-                ])
+    @cached_method
+    def rational_algebra(self):
+        if self.base_ring() is QQ:
+            return self
+
+        return cartesian_product([
+            r.rational_algebra() for r in self.cartesian_factors()
+        ])
 
 
 RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
 
 def random_eja(max_dimension=None, *args, **kwargs):
+    r"""
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import random_eja
+
+    TESTS::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: J = random_eja(max_dimension=n, field=QQ, orthonormalize=False)
+        sage: J.dimension() <= n
+        True
+
+    """
     # Use the ConcreteEJA default as the total upper bound (regardless
     # of any whether or not any individual factors set a lower limit).
     if max_dimension is None: