]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: refactor all matrix classes upwards (note: everything broken).
[sage.d.git] / mjo / eja / eja_algebra.py
index aef11acd1dbb6dc60b4765be1790f2fb71c5a1be..106a0cddec06355a952e71d75699677b25dd7da9 100644 (file)
@@ -32,22 +32,25 @@ for these simple algebras:
   * :class:`RealSymmetricEJA`
   * :class:`ComplexHermitianEJA`
   * :class:`QuaternionHermitianEJA`
+  * :class:`OctonionHermitianEJA`
 
-Missing from this list is the algebra of three-by-three octononion
-Hermitian matrices, as there is (as of yet) no implementation of the
-octonions in SageMath. In addition to these, we provide two other
-example constructions,
+In addition to these, we provide two other example constructions,
 
+  * :class:`JordanSpinEJA`
   * :class:`HadamardEJA`
+  * :class:`AlbertEJA`
   * :class:`TrivialEJA`
 
 The Jordan spin algebra is a bilinear form algebra where the bilinear
 form is the identity. The Hadamard EJA is simply a Cartesian product
-of one-dimensional spin algebras. And last but not least, the trivial
-EJA is exactly what you think. Cartesian products of these are also
-supported using the usual ``cartesian_product()`` function; as a
-result, we support (up to isomorphism) all Euclidean Jordan algebras
-that don't involve octonions.
+of one-dimensional spin algebras. The Albert EJA is simply a special
+case of the :class:`OctonionHermitianEJA` where the matrices are
+three-by-three and the resulting space has dimension 27. And
+last/least, the trivial EJA is exactly what you think it is; it could
+also be obtained by constructing a dimension-zero instance of any of
+the other algebras. Cartesian products of these are also supported
+using the usual ``cartesian_product()`` function; as a result, we
+support (up to isomorphism) all Euclidean Jordan algebras.
 
 SETUP::
 
@@ -59,8 +62,6 @@ EXAMPLES::
     Euclidean Jordan algebra of dimension...
 """
 
-from itertools import repeat
-
 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
 from sage.categories.magmatic_algebras import MagmaticAlgebras
 from sage.categories.sets_cat import cartesian_product
@@ -171,6 +172,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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()
@@ -1018,7 +1024,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             Full MatrixSpace of 4 by 4 dense matrices over Rational Field
             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():
@@ -1543,7 +1551,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
 class RationalBasisEJA(FiniteDimensionalEJA):
     r"""
-    New class for algebras whose supplied basis elements have all rational entries.
+    Algebras whose supplied basis elements have all rational entries.
 
     SETUP::
 
@@ -1574,7 +1582,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 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,
@@ -1726,112 +1738,147 @@ class ConcreteEJA(FiniteDimensionalEJA):
 
 class MatrixEJA:
     @staticmethod
-    def jordan_product(X,Y):
-        return (X*Y + Y*X)/2
-
-    @staticmethod
-    def trace_inner_product(X,Y):
-        r"""
-        A trace inner-product for matrices that aren't embedded in the
-        reals. It takes MATRICES as arguments, not EJA elements.
+    def _denormalized_basis(A):
         """
-        return (X*Y).trace().real()
+        Returns a basis for the space of complex Hermitian n-by-n matrices.
 
-class RealEmbeddedMatrixEJA(MatrixEJA):
-    @staticmethod
-    def dimension_over_reals():
-        r"""
-        The dimension of this matrix's base ring 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.
 
-        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.
+        SETUP::
 
-        This is used to determine the size of the matrix returned from
-        :meth:`real_embed`, among other things.
-        """
-        raise NotImplementedError
+            sage: from mjo.hurwitz import (ComplexMatrixAlgebra,
+            ....:                          QuaternionMatrixAlgebra,
+            ....:                          OctonionMatrixAlgebra)
+            sage: from mjo.eja.eja_algebra import MatrixEJA
 
-    @classmethod
-    def real_embed(cls,M):
-        """
-        Embed the matrix ``M`` into a space of real matrices.
+        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)
 
-    @classmethod
-    def trace_inner_product(cls,X,Y):
+        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 trace_inner_product(X,Y):
         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 (ComplexHermitianEJA,
-            ....:                                  QuaternionHermitianEJA)
+            sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
+            ....:                                  ComplexHermitianEJA,
+            ....:                                  QuaternionHermitianEJA,
+            ....:                                  OctonionHermitianEJA)
 
         EXAMPLES::
 
-            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 = 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 = 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 = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2
+
+        ::
+
+            sage: J = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2
+
+        ::
+
+            sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2
 
         """
-        # This does in fact compute the real part of the trace.
-        # If we compute the trace of e.g. a complex matrix M,
-        # then we do so by adding up its diagonal entries --
-        # call them z_1 through z_n. The real embedding of z_1
-        # will be a 2-by-2 REAL matrix [a, b; -b, a] whose trace
-        # as a REAL matrix will be 2*a = 2*Re(z_1). And so forth.
-        return (X*Y).trace()/cls.dimension_over_reals()
+        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]
+
+        # Works for real and complex numbers.
+        return tr.real()
 
-class RealSymmetricEJA(ConcreteEJA, RationalBasisEJA, MatrixEJA):
+
+
+class RealSymmetricEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
     """
     The rank-n simple EJA consisting of real symmetric n-by-n
     matrices, the usual symmetric Jordan product, and the trace inner
@@ -1896,38 +1943,6 @@ class RealSymmetricEJA(ConcreteEJA, RationalBasisEJA, MatrixEJA):
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """
-    @classmethod
-    def _denormalized_basis(cls, n, field):
-        """
-        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,ZZ)
-            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(field, 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
@@ -1945,175 +1960,22 @@ class RealSymmetricEJA(ConcreteEJA, RationalBasisEJA, MatrixEJA):
         # 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,field),
+        A = MatrixSpace(field, n)
+        super().__init__(self._denormalized_basis(A),
                          self.jordan_product,
                          self.trace_inner_product,
                          field=field,
-                         associative=associative,
                          **kwargs)
 
         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         self.rank.set_cache(n)
-        idV = self.matrix_space().one()
-        self.one.set_cache(self(idV))
-
-
+        self.one.set_cache(self(A.one()))
 
-class ComplexMatrixEJA(RealEmbeddedMatrixEJA):
-    # 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:
-
-        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, RationalBasisEJA, ComplexMatrixEJA):
+class ComplexHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
     """
     The rank-n simple EJA consisting of complex Hermitian n-by-n
     matrices over the real numbers, the usual symmetric Jordan product,
@@ -2170,90 +2032,24 @@ class ComplexHermitianEJA(ConcreteEJA, RationalBasisEJA, ComplexMatrixEJA):
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """
-
-    @classmethod
-    def _denormalized_basis(cls, n, field):
-        """
-        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,ZZ)
-            sage: all( M.is_symmetric() for M in  B)
-            True
-
-        """
-        R = PolynomialRing(ZZ, 'z')
-        z = R.gen()
-        F = ZZ.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 the entries, we can drop back to the
-        # desired real "field" instead of the extension "F".
-        return tuple( s.change_ring(field) for s in S )
-
-
     def __init__(self, n, 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
 
-        associative = False
-        if n <= 1:
-            associative = True
-
-        super().__init__(self._denormalized_basis(n,field),
+        from mjo.hurwitz import ComplexMatrixAlgebra
+        A = ComplexMatrixAlgebra(n, scalars=field)
+        super().__init__(self._denormalized_basis(A),
                          self.jordan_product,
                          self.trace_inner_product,
                          field=field,
-                         associative=associative,
                          **kwargs)
+
         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         self.rank.set_cache(n)
-        idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
-        self.one.set_cache(self(idV))
+        self.one.set_cache(self(A.one()))
 
     @staticmethod
     def _max_random_instance_size():
@@ -2267,157 +2063,8 @@ class ComplexHermitianEJA(ConcreteEJA, RationalBasisEJA, ComplexMatrixEJA):
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
         return cls(n, **kwargs)
 
-class QuaternionMatrixEJA(RealEmbeddedMatrixEJA):
-
-    # 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,
-                             RationalBasisEJA,
-                             QuaternionMatrixEJA):
+class QuaternionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
     r"""
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the
@@ -2474,100 +2121,24 @@ class QuaternionHermitianEJA(ConcreteEJA,
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """
-    @classmethod
-    def _denormalized_basis(cls, n, field):
-        """
-        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::
-
-            sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
-
-        TESTS::
-
-            sage: set_random_seed()
-            sage: n = ZZ.random_element(1,5)
-            sage: B = QuaternionHermitianEJA._denormalized_basis(n,ZZ)
-            sage: all( M.is_symmetric() for M in B )
-            True
-
-        """
-        Q = QuaternionAlgebra(QQ,-1,-1)
-        I,J,K = Q.gens()
-
-        # 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 the entries, we can drop back to the
-        # desired real "field" instead of the quaternion algebra "Q".
-        return tuple( s.change_ring(field) for s in S )
-
-
     def __init__(self, n, 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
 
-        associative = False
-        if n <= 1:
-            associative = True
-
-        super().__init__(self._denormalized_basis(n,field),
+        from mjo.hurwitz import QuaternionMatrixAlgebra
+        A = QuaternionMatrixAlgebra(n, scalars=field)
+        super().__init__(self._denormalized_basis(A),
                          self.jordan_product,
                          self.trace_inner_product,
                          field=field,
-                         associative=associative,
                          **kwargs)
 
         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         self.rank.set_cache(n)
-        idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
-        self.one.set_cache(self(idV))
+        self.one.set_cache(self(A.one()))
 
 
     @staticmethod
@@ -2585,7 +2156,7 @@ class QuaternionHermitianEJA(ConcreteEJA,
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
         return cls(n, **kwargs)
 
-class OctonionHermitianEJA(FiniteDimensionalEJA, MatrixEJA):
+class OctonionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
     r"""
     SETUP::
 
@@ -2668,7 +2239,23 @@ class OctonionHermitianEJA(FiniteDimensionalEJA, MatrixEJA):
         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.
+        """
+        return 1 # Dimension 1
+
+    @classmethod
+    def random_instance(cls, **kwargs):
+        """
+        Return a random instance of this type of algebra.
+        """
+        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        return cls(n, **kwargs)
+
     def __init__(self, n, field=AA, **kwargs):
         if n > 3:
             # Otherwise we don't get an EJA.
@@ -2678,7 +2265,9 @@ class OctonionHermitianEJA(FiniteDimensionalEJA, MatrixEJA):
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
-        super().__init__(self._denormalized_basis(n,field),
+        from mjo.hurwitz import OctonionMatrixAlgebra
+        A = OctonionMatrixAlgebra(n, scalars=field)
+        super().__init__(self._denormalized_basis(A),
                          self.jordan_product,
                          self.trace_inner_product,
                          field=field,
@@ -2688,81 +2277,39 @@ class OctonionHermitianEJA(FiniteDimensionalEJA, MatrixEJA):
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         self.rank.set_cache(n)
-        idV = self.matrix_space().one()
-        self.one.set_cache(self(idV))
-
-
-    @classmethod
-    def _denormalized_basis(cls, n, field):
-        """
-        Returns a basis for the space of octonion Hermitian n-by-n
-        matrices.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
+        self.one.set_cache(self(A.one()))
 
-        EXAMPLES::
 
-            sage: B = OctonionHermitianEJA._denormalized_basis(3,QQ)
-            sage: all( M.is_hermitian() for M in B )
-            True
-            sage: len(B)
-            27
-
-        """
-        from mjo.octonions import OctonionMatrixAlgebra
-        MS = OctonionMatrixAlgebra(n, scalars=field)
-        es = MS.entry_algebra().gens()
-
-        basis = []
-        for i in range(n):
-            for j in range(i+1):
-                if i == j:
-                    E_ii = MS.monomial( (i,j,es[0]) )
-                    basis.append(E_ii)
-                else:
-                    for e in es:
-                        E_ij  = MS.monomial( (i,j,e)             )
-                        ec = e.conjugate()
-                        # If the conjugate has a negative sign in front
-                        # of it, (j,i,ec) won't be a monomial!
-                        if (j,i,ec) in MS.indices():
-                            E_ij += MS.monomial( (j,i,ec) )
-                        else:
-                            E_ij -= MS.monomial( (j,i,-ec) )
-                        basis.append(E_ij)
-
-        return tuple( basis )
+class AlbertEJA(OctonionHermitianEJA):
+    r"""
+    The Albert algebra is the algebra of three-by-three Hermitian
+    matrices whose entries are octonions.
 
-    @staticmethod
-    def trace_inner_product(X,Y):
-        r"""
-        The octonions don't know that the reals are embedded in them,
-        so we have to take the e0 component ourselves.
+    SETUP::
 
-        SETUP::
+        sage: from mjo.eja.eja_algebra import AlbertEJA
 
-            sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
+    EXAMPLES::
 
-        TESTS::
+        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
 
-            sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
-            sage: I = J.one().to_matrix()
-            sage: J.trace_inner_product(I, -I)
-            -2
+    """
+    def __init__(self, *args, **kwargs):
+        super().__init__(3, *args, **kwargs)
 
-        """
-        return (X*Y).trace().real().coefficient(0)
 
-class HadamardEJA(ConcreteEJA, RationalBasisEJA):
+class HadamardEJA(RationalBasisEJA, ConcreteEJA):
     """
-    Return the Euclidean Jordan Algebra corresponding to the set
-    `R^n` under the Hadamard product.
+    Return the Euclidean Jordan algebra on `R^n` with the Hadamard
+    (pointwise real-number multiplication) Jordan product and the
+    usual inner-product.
 
-    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.
+    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::
 
@@ -2793,7 +2340,6 @@ class HadamardEJA(ConcreteEJA, RationalBasisEJA):
 
         sage: HadamardEJA(3, prefix='r').gens()
         (r0, r1, r2)
-
     """
     def __init__(self, n, field=AA, **kwargs):
         if n == 0:
@@ -2847,7 +2393,7 @@ class HadamardEJA(ConcreteEJA, RationalBasisEJA):
         return cls(n, **kwargs)
 
 
-class BilinearFormEJA(ConcreteEJA, RationalBasisEJA):
+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 =
@@ -3095,7 +2641,7 @@ class JordanSpinEJA(BilinearFormEJA):
         return cls(n, **kwargs)
 
 
-class TrivialEJA(ConcreteEJA, RationalBasisEJA):
+class TrivialEJA(RationalBasisEJA, ConcreteEJA):
     """
     The trivial Euclidean Jordan algebra consisting of only a zero element.
 
@@ -3380,9 +2926,17 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         Return the space that our matrix basis lives in as a Cartesian
         product.
 
+        We don't simply use the ``cartesian_product()`` functor here
+        because it acts differently on SageMath MatrixSpaces and our
+        custom MatrixAlgebras, which are CombinatorialFreeModules. We
+        always want the result to be represented (and indexed) as
+        an ordered tuple.
+
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
+            ....:                                  HadamardEJA,
+            ....:                                  OctonionHermitianEJA,
             ....:                                  RealSymmetricEJA)
 
         EXAMPLES::
@@ -3395,10 +2949,44 @@ class CartesianProductEJA(FiniteDimensionalEJA):
             matrices over Algebraic Real Field, Full MatrixSpace of 2
             by 2 dense matrices over Algebraic Real Field)
 
+        ::
+
+            sage: J1 = ComplexHermitianEJA(1)
+            sage: J2 = ComplexHermitianEJA(1)
+            sage: J = cartesian_product([J1,J2])
+            sage: J.one().to_matrix()[0]
+            [1 0]
+            [0 1]
+            sage: J.one().to_matrix()[1]
+            [1 0]
+            [0 1]
+
+        ::
+
+            sage: J1 = OctonionHermitianEJA(1)
+            sage: J2 = OctonionHermitianEJA(1)
+            sage: J = cartesian_product([J1,J2])
+            sage: J.one().to_matrix()[0]
+            +----+
+            | e0 |
+            +----+
+            sage: J.one().to_matrix()[1]
+            +----+
+            | e0 |
+            +----+
+
         """
-        from sage.categories.cartesian_product import cartesian_product
-        return cartesian_product( [J.matrix_space()
-                                   for J in self.cartesian_factors()] )
+        scalars = self.cartesian_factor(0).base_ring()
+
+        # This category isn't perfect, but is good enough for what we
+        # need to do.
+        cat = MagmaticAlgebras(scalars).FiniteDimensional().WithBasis()
+        cat = cat.Unital().CartesianProducts()
+        factors = tuple( J.matrix_space() for J in self.cartesian_factors() )
+
+        from sage.sets.cartesian_product import CartesianProduct
+        return CartesianProduct(factors, cat)
+
 
     @cached_method
     def cartesian_projection(self, i):
@@ -3600,7 +3188,9 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
 
     SETUP::
 
-        sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+        sage: from mjo.eja.eja_algebra import (HadamardEJA,
+        ....:                                  JordanSpinEJA,
+        ....:                                  OctonionHermitianEJA,
         ....:                                  RealSymmetricEJA)
 
     EXAMPLES:
@@ -3617,15 +3207,32 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
         sage: J.rank()
         5
 
+    TESTS:
+
+    The ``cartesian_product()`` function only uses the first factor to
+    decide where the result will live; thus we have to be careful to
+    check that all factors do indeed have a `_rational_algebra` member
+    before we try to access it::
+
+        sage: J1 = OctonionHermitianEJA(1) # no rational basis
+        sage: J2 = HadamardEJA(2)
+        sage: cartesian_product([J1,J2])
+        Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
+        (+) Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
+        sage: cartesian_product([J2,J1])
+        Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
+        (+) Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
+
     """
     def __init__(self, algebras, **kwargs):
         CartesianProductEJA.__init__(self, algebras, **kwargs)
 
         self._rational_algebra = None
         if self.vector_space().base_field() is not QQ:
-            self._rational_algebra = cartesian_product([
-                r._rational_algebra for r in algebras
-            ])
+            if all( hasattr(r, "_rational_algebra") for r in algebras ):
+                self._rational_algebra = cartesian_product([
+                    r._rational_algebra for r in algebras
+                ])
 
 
 RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA