]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: add the AlbertEJA class.
[sage.d.git] / mjo / eja / eja_algebra.py
index 558ff6b21f62e664ca1bf6e8730eaeb81a0f25d7..ca3df5fea55949225a51f8255a2c13a33c37d9f2 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,13 +62,10 @@ 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
-from sage.combinat.free_module import (CombinatorialFreeModule,
-                                       CombinatorialFreeModule_CartesianProduct)
+from sage.combinat.free_module import CombinatorialFreeModule
 from sage.matrix.constructor import matrix
 from sage.matrix.matrix_space import MatrixSpace
 from sage.misc.cachefunc import cached_method
@@ -142,19 +142,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                  cartesian_product=False,
                  check_field=True,
                  check_axioms=True,
-                 prefix='e'):
-
-        # Keep track of whether or not the matrix basis consists of
-        # tuples, since we need special cases for them damned near
-        # everywhere.  This is INDEPENDENT of whether or not the
-        # algebra is a cartesian product, since a subalgebra of a
-        # cartesian product will have a basis of tuples, but will not
-        # in general itself be a cartesian product algebra.
-        self._matrix_basis_is_cartesian = False
+                 prefix="b"):
+
         n = len(basis)
-        if n > 0:
-            if hasattr(basis[0], 'cartesian_factors'):
-                self._matrix_basis_is_cartesian = True
 
         if check_field:
             if not field.is_subring(RR):
@@ -163,21 +153,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                 # we've specified a real embedding.
                 raise ValueError("scalar field is not real")
 
-        # If the basis given to us wasn't over the field that it's
-        # supposed to be over, fix that. Or, you know, crash.
-        if not cartesian_product:
-            # The field for a cartesian product algebra comes from one
-            # of its factors and is the same for all factors, so
-            # there's no need to "reapply" it on product algebras.
-            if self._matrix_basis_is_cartesian:
-                # OK since if n == 0, the basis does not consist of tuples.
-                P = basis[0].parent()
-                basis = tuple( P(tuple(b_i.change_ring(field) for b_i in b))
-                               for b in basis )
-            else:
-                basis = tuple( b.change_ring(field) for b in basis )
-
-
         if check_axioms:
             # Check commutativity of the Jordan and inner-products.
             # This has to be done before we build the multiplication
@@ -231,7 +206,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # ambient vector space V that our (vectorized) basis lives in,
         # as well as a subspace W of V spanned by those (vectorized)
         # basis elements. The W-coordinates are the coefficients that
-        # we see in things like x = 1*e1 + 2*e2.
+        # we see in things like x = 1*b1 + 2*b2.
         vector_basis = basis
 
         degree = 0
@@ -358,16 +333,16 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: set_random_seed()
             sage: J = random_eja()
             sage: n = J.dimension()
-            sage: ei = J.zero()
-            sage: ej = J.zero()
-            sage: ei_ej = J.zero()*J.zero()
+            sage: bi = J.zero()
+            sage: bj = J.zero()
+            sage: bi_bj = J.zero()*J.zero()
             sage: if n > 0:
             ....:     i = ZZ.random_element(n)
             ....:     j = ZZ.random_element(n)
-            ....:     ei = J.monomial(i)
-            ....:     ej = J.monomial(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
 
         """
@@ -639,7 +614,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: J2 = RealSymmetricEJA(2)
             sage: J = cartesian_product([J1,J2])
             sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
-            e1 + e5
+            b1 + b5
 
         TESTS:
 
@@ -914,15 +889,15 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: J = JordanSpinEJA(4)
             sage: J.multiplication_table()
             +----++----+----+----+----+
-            | *  || e0 | e1 | e2 | e3 |
+            | *  || b0 | b1 | b2 | b3 |
             +====++====+====+====+====+
-            | e0 || e0 | e1 | e2 | e3 |
+            | b0 || b0 | b1 | b2 | b3 |
             +----++----+----+----+----+
-            | e1 || e1 | e0 | 0  | 0  |
+            | b1 || b1 | b0 | 0  | 0  |
             +----++----+----+----+----+
-            | e2 || e2 | 0  | e0 | 0  |
+            | b2 || b2 | 0  | b0 | 0  |
             +----++----+----+----+----+
-            | e3 || e3 | 0  | 0  | e0 |
+            | b3 || b3 | 0  | 0  | b0 |
             +----++----+----+----+----+
 
         """
@@ -976,7 +951,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: J = RealSymmetricEJA(2)
             sage: J.basis()
-            Finite family {0: e0, 1: e1, 2: e2}
+            Finite family {0: b0, 1: b1, 2: b2}
             sage: J.matrix_basis()
             (
             [1 0]  [                  0 0.7071067811865475?]  [0 0]
@@ -987,7 +962,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: J = JordanSpinEJA(2)
             sage: J.basis()
-            Finite family {0: e0, 1: e1}
+            Finite family {0: b0, 1: b1}
             sage: J.matrix_basis()
             (
             [1]  [0]
@@ -1069,20 +1044,20 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: J = HadamardEJA(5)
             sage: J.one()
-            e0 + e1 + e2 + e3 + e4
+            b0 + b1 + b2 + b3 + b4
 
         The unit element in the Hadamard EJA is inherited in the
         subalgebras generated by its elements::
 
             sage: J = HadamardEJA(5)
             sage: J.one()
-            e0 + e1 + e2 + e3 + e4
+            b0 + b1 + b2 + b3 + b4
             sage: x = sum(J.gens())
             sage: A = x.subalgebra_generated_by(orthonormalize=False)
             sage: A.one()
-            f0
+            c0
             sage: A.one().superalgebra_element()
-            e0 + e1 + e2 + e3 + e4
+            b0 + b1 + b2 + b3 + b4
 
         TESTS:
 
@@ -1569,7 +1544,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::
 
@@ -1682,7 +1657,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 )
 
-class ConcreteEJA(RationalBasisEJA):
+class ConcreteEJA(FiniteDimensionalEJA):
     r"""
     A class for the Euclidean Jordan algebras that we know by name.
 
@@ -1751,6 +1726,19 @@ class ConcreteEJA(RationalBasisEJA):
 
 
 class MatrixEJA:
+    @staticmethod
+    def jordan_product(X,Y):
+        return (X*Y + Y*X)/2
+
+    @staticmethod
+    def trace_inner_product(X,Y):
+        r"""
+        A trace inner-product for matrices that aren't embedded in the
+        reals. It takes MATRICES as arguments, not EJA elements.
+        """
+        return (X*Y).trace().real()
+
+class RealEmbeddedMatrixEJA(MatrixEJA):
     @staticmethod
     def dimension_over_reals():
         r"""
@@ -1796,9 +1784,6 @@ class MatrixEJA:
             raise ValueError("the matrix 'M' must be a real embedding")
         return M
 
-    @staticmethod
-    def jordan_product(X,Y):
-        return (X*Y + Y*X)/2
 
     @classmethod
     def trace_inner_product(cls,X,Y):
@@ -1807,29 +1792,11 @@ class MatrixEJA:
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
-            ....:                                  ComplexHermitianEJA,
+            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
             ....:                                  QuaternionHermitianEJA)
 
         EXAMPLES::
 
-        This gives the same answer as it would if we computed the trace
-        from the unembedded (original) matrices::
-
-            sage: set_random_seed()
-            sage: J = RealSymmetricEJA.random_instance()
-            sage: x,y = J.random_elements(2)
-            sage: Xe = x.to_matrix()
-            sage: Ye = y.to_matrix()
-            sage: X = J.real_unembed(Xe)
-            sage: Y = J.real_unembed(Ye)
-            sage: expected = (X*Y).trace()
-            sage: actual = J.trace_inner_product(Xe,Ye)
-            sage: actual == expected
-            True
-
-        ::
-
             sage: set_random_seed()
             sage: J = ComplexHermitianEJA.random_instance()
             sage: x,y = J.random_elements(2)
@@ -1857,27 +1824,15 @@ class MatrixEJA:
             True
 
         """
-        Xu = cls.real_unembed(X)
-        Yu = cls.real_unembed(Y)
-        tr = (Xu*Yu).trace()
-
-        try:
-            # Works in QQ, AA, RDF, et cetera.
-            return tr.real()
-        except AttributeError:
-            # A quaternion doesn't have a real() method, but does
-            # have coefficient_tuple() method that returns the
-            # coefficients of 1, i, j, and k -- in that order.
-            return tr.coefficient_tuple()[0]
-
-
-class RealMatrixEJA(MatrixEJA):
-    @staticmethod
-    def dimension_over_reals():
-        return 1
-
+        # This does in fact compute the real part of the trace.
+        # If we compute the trace of e.g. a complex matrix M,
+        # then we do so by adding up its diagonal entries --
+        # call them z_1 through z_n. The real embedding of z_1
+        # will be a 2-by-2 REAL matrix [a, b; -b, a] whose trace
+        # as a REAL matrix will be 2*a = 2*Re(z_1). And so forth.
+        return (X*Y).trace()/cls.dimension_over_reals()
 
-class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
+class RealSymmetricEJA(ConcreteEJA, RationalBasisEJA, MatrixEJA):
     """
     The rank-n simple EJA consisting of real symmetric n-by-n
     matrices, the usual symmetric Jordan product, and the trace inner
@@ -1890,13 +1845,13 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
     EXAMPLES::
 
         sage: J = RealSymmetricEJA(2)
-        sage: e0, e1, e2 = J.gens()
-        sage: e0*e0
-        e0
-        sage: e1*e1
-        1/2*e0 + 1/2*e2
-        sage: e2*e2
-        e2
+        sage: b0, b1, b2 = J.gens()
+        sage: b0*b0
+        b0
+        sage: b1*b1
+        1/2*b0 + 1/2*b2
+        sage: b2*b2
+        b2
 
     In theory, our "field" can be any subfield of the reals::
 
@@ -1943,7 +1898,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
 
     """
     @classmethod
-    def _denormalized_basis(cls, n):
+    def _denormalized_basis(cls, n, field):
         """
         Return a basis for the space of real symmetric n-by-n matrices.
 
@@ -1955,7 +1910,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
-            sage: B = RealSymmetricEJA._denormalized_basis(n)
+            sage: B = RealSymmetricEJA._denormalized_basis(n,ZZ)
             sage: all( M.is_symmetric() for M in  B)
             True
 
@@ -1965,7 +1920,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
         S = []
         for i in range(n):
             for j in range(i+1):
-                Eij = matrix(ZZ, n, lambda k,l: k==i and l==j)
+                Eij = matrix(field, n, lambda k,l: k==i and l==j)
                 if i == j:
                     Sij = Eij
                 else:
@@ -1986,7 +1941,7 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
         return cls(n, **kwargs)
 
-    def __init__(self, n, **kwargs):
+    def __init__(self, n, field=AA, **kwargs):
         # We know this is a valid EJA, but will double-check
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
@@ -1995,9 +1950,10 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
         if n <= 1:
             associative = True
 
-        super().__init__(self._denormalized_basis(n),
+        super().__init__(self._denormalized_basis(n,field),
                          self.jordan_product,
                          self.trace_inner_product,
+                         field=field,
                          associative=associative,
                          **kwargs)
 
@@ -2005,12 +1961,12 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         self.rank.set_cache(n)
-        idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
+        idV = self.matrix_space().one()
         self.one.set_cache(self(idV))
 
 
 
-class ComplexMatrixEJA(MatrixEJA):
+class ComplexMatrixEJA(RealEmbeddedMatrixEJA):
     # A manual dictionary-cache for the complex_extension() method,
     # since apparently @classmethods can't also be @cached_methods.
     _complex_extension = {}
@@ -2158,7 +2114,7 @@ class ComplexMatrixEJA(MatrixEJA):
         return matrix(F, n/d, elements)
 
 
-class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
+class ComplexHermitianEJA(ConcreteEJA, RationalBasisEJA, ComplexMatrixEJA):
     """
     The rank-n simple EJA consisting of complex Hermitian n-by-n
     matrices over the real numbers, the usual symmetric Jordan product,
@@ -2217,7 +2173,7 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
     """
 
     @classmethod
-    def _denormalized_basis(cls, n):
+    def _denormalized_basis(cls, n, field):
         """
         Returns a basis for the space of complex Hermitian n-by-n matrices.
 
@@ -2235,15 +2191,14 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
-            sage: B = ComplexHermitianEJA._denormalized_basis(n)
+            sage: B = ComplexHermitianEJA._denormalized_basis(n,ZZ)
             sage: all( M.is_symmetric() for M in  B)
             True
 
         """
-        field = ZZ
-        R = PolynomialRing(field, 'z')
+        R = PolynomialRing(ZZ, 'z')
         z = R.gen()
-        F = field.extension(z**2 + 1, 'I')
+        F = ZZ.extension(z**2 + 1, 'I')
         I = F.gen(1)
 
         # This is like the symmetric case, but we need to be careful:
@@ -2274,12 +2229,12 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
                 # "erase" E_ij
                 Eij[i,j] = 0
 
-        # Since we embedded these, we can drop back to the "field" that we
-        # started with instead of the complex extension "F".
+        # Since we embedded the entries, we can drop back to the
+        # desired real "field" instead of the extension "F".
         return tuple( s.change_ring(field) for s in S )
 
 
-    def __init__(self, n, **kwargs):
+    def __init__(self, n, field=AA, **kwargs):
         # We know this is a valid EJA, but will double-check
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
@@ -2288,9 +2243,10 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
         if n <= 1:
             associative = True
 
-        super().__init__(self._denormalized_basis(n),
+        super().__init__(self._denormalized_basis(n,field),
                          self.jordan_product,
                          self.trace_inner_product,
+                         field=field,
                          associative=associative,
                          **kwargs)
         # TODO: this could be factored out somehow, but is left here
@@ -2312,7 +2268,7 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
         return cls(n, **kwargs)
 
-class QuaternionMatrixEJA(MatrixEJA):
+class QuaternionMatrixEJA(RealEmbeddedMatrixEJA):
 
     # A manual dictionary-cache for the quaternion_extension() method,
     # since apparently @classmethods can't also be @cached_methods.
@@ -2460,7 +2416,9 @@ class QuaternionMatrixEJA(MatrixEJA):
         return matrix(Q, n/d, elements)
 
 
-class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
+class QuaternionHermitianEJA(ConcreteEJA,
+                             RationalBasisEJA,
+                             QuaternionMatrixEJA):
     r"""
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the
@@ -2518,7 +2476,7 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
 
     """
     @classmethod
-    def _denormalized_basis(cls, n):
+    def _denormalized_basis(cls, n, field):
         """
         Returns a basis for the space of quaternion Hermitian n-by-n matrices.
 
@@ -2536,12 +2494,11 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
-            sage: B = QuaternionHermitianEJA._denormalized_basis(n)
+            sage: B = QuaternionHermitianEJA._denormalized_basis(n,ZZ)
             sage: all( M.is_symmetric() for M in B )
             True
 
         """
-        field = ZZ
         Q = QuaternionAlgebra(QQ,-1,-1)
         I,J,K = Q.gens()
 
@@ -2585,12 +2542,12 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
                 # "erase" E_ij
                 Eij[i,j] = 0
 
-        # Since we embedded these, we can drop back to the "field" that we
-        # started with instead of the quaternion algebra "Q".
+        # Since we embedded the entries, we can drop back to the
+        # desired real "field" instead of the quaternion algebra "Q".
         return tuple( s.change_ring(field) for s in S )
 
 
-    def __init__(self, n, **kwargs):
+    def __init__(self, n, field=AA, **kwargs):
         # We know this is a valid EJA, but will double-check
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
@@ -2599,9 +2556,10 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
         if n <= 1:
             associative = True
 
-        super().__init__(self._denormalized_basis(n),
+        super().__init__(self._denormalized_basis(n,field),
                          self.jordan_product,
                          self.trace_inner_product,
+                         field=field,
                          associative=associative,
                          **kwargs)
 
@@ -2628,8 +2586,215 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
         return cls(n, **kwargs)
 
+class OctonionHermitianEJA(ConcreteEJA, MatrixEJA):
+    r"""
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
+        ....:                                  OctonionHermitianEJA)
+
+    EXAMPLES:
+
+    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: b = OctonionHermitianEJA._denormalized_basis(2,QQ)
+        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
 
-class HadamardEJA(ConcreteEJA):
+    """
+    @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.
+            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
+
+        super().__init__(self._denormalized_basis(n,field),
+                         self.jordan_product,
+                         self.trace_inner_product,
+                         field=field,
+                         **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))
+
+
+    @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
+
+        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 )
+
+    @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::
+
+            sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
+
+        TESTS::
+
+            sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2
+
+        """
+        return (X*Y).trace().real().coefficient(0)
+
+
+class AlbertEJA(OctonionHermitianEJA):
+    r"""
+    The Albert algebra is the algebra of three-by-three Hermitian
+    matrices whose entries are octonions.
+
+    SETUP::
+
+        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
+
+    """
+    def __init__(self, *args, **kwargs):
+        super().__init__(3, *args, **kwargs)
+
+
+class HadamardEJA(ConcreteEJA, RationalBasisEJA):
     """
     Return the Euclidean Jordan Algebra corresponding to the set
     `R^n` under the Hadamard product.
@@ -2647,19 +2812,19 @@ class HadamardEJA(ConcreteEJA):
     This multiplication table can be verified by hand::
 
         sage: J = HadamardEJA(3)
-        sage: e0,e1,e2 = J.gens()
-        sage: e0*e0
-        e0
-        sage: e0*e1
+        sage: b0,b1,b2 = J.gens()
+        sage: b0*b0
+        b0
+        sage: b0*b1
         0
-        sage: e0*e2
+        sage: b0*b2
         0
-        sage: e1*e1
-        e1
-        sage: e1*e2
+        sage: b1*b1
+        b1
+        sage: b1*b2
         0
-        sage: e2*e2
-        e2
+        sage: b2*b2
+        b2
 
     TESTS:
 
@@ -2669,7 +2834,7 @@ class HadamardEJA(ConcreteEJA):
         (r0, r1, r2)
 
     """
-    def __init__(self, n, **kwargs):
+    def __init__(self, n, field=AA, **kwargs):
         if n == 0:
             jordan_product = lambda x,y: x
             inner_product = lambda x,y: x
@@ -2690,10 +2855,12 @@ class HadamardEJA(ConcreteEJA):
         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
-        column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )
+        column_basis = tuple( b.column()
+                              for b in FreeModule(field, n).basis() )
         super().__init__(column_basis,
                          jordan_product,
                          inner_product,
+                         field=field,
                          associative=True,
                          **kwargs)
         self.rank.set_cache(n)
@@ -2719,7 +2886,7 @@ class HadamardEJA(ConcreteEJA):
         return cls(n, **kwargs)
 
 
-class BilinearFormEJA(ConcreteEJA):
+class BilinearFormEJA(ConcreteEJA, RationalBasisEJA):
     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 =
@@ -2801,7 +2968,7 @@ class BilinearFormEJA(ConcreteEJA):
         True
 
     """
-    def __init__(self, B, **kwargs):
+    def __init__(self, B, field=AA, **kwargs):
         # The matrix "B" is supplied by the user in most cases,
         # so it makes sense to check whether or not its positive-
         # definite unless we are specifically asked not to...
@@ -2829,7 +2996,8 @@ class BilinearFormEJA(ConcreteEJA):
             return P([z0] + zbar.list())
 
         n = B.nrows()
-        column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )
+        column_basis = tuple( b.column()
+                              for b in FreeModule(field, n).basis() )
 
         # TODO: I haven't actually checked this, but it seems legit.
         associative = False
@@ -2839,6 +3007,7 @@ class BilinearFormEJA(ConcreteEJA):
         super().__init__(column_basis,
                          jordan_product,
                          inner_product,
+                         field=field,
                          associative=associative,
                          **kwargs)
 
@@ -2900,20 +3069,20 @@ class JordanSpinEJA(BilinearFormEJA):
     This multiplication table can be verified by hand::
 
         sage: J = JordanSpinEJA(4)
-        sage: e0,e1,e2,e3 = J.gens()
-        sage: e0*e0
-        e0
-        sage: e0*e1
-        e1
-        sage: e0*e2
-        e2
-        sage: e0*e3
-        e3
-        sage: e1*e2
+        sage: b0,b1,b2,b3 = J.gens()
+        sage: b0*b0
+        b0
+        sage: b0*b1
+        b1
+        sage: b0*b2
+        b2
+        sage: b0*b3
+        b3
+        sage: b1*b2
         0
-        sage: e1*e3
+        sage: b1*b3
         0
-        sage: e2*e3
+        sage: b2*b3
         0
 
     We can change the generator prefix::
@@ -2934,7 +3103,7 @@ class JordanSpinEJA(BilinearFormEJA):
             True
 
     """
-    def __init__(self, n, **kwargs):
+    def __init__(self, n, *args, **kwargs):
         # This is a special case of the BilinearFormEJA with the
         # identity matrix as its bilinear form.
         B = matrix.identity(ZZ, n)
@@ -2945,7 +3114,7 @@ class JordanSpinEJA(BilinearFormEJA):
 
         # 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():
@@ -2965,7 +3134,7 @@ class JordanSpinEJA(BilinearFormEJA):
         return cls(n, **kwargs)
 
 
-class TrivialEJA(ConcreteEJA):
+class TrivialEJA(ConcreteEJA, RationalBasisEJA):
     """
     The trivial Euclidean Jordan algebra consisting of only a zero element.
 
@@ -3125,23 +3294,23 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
         sage: J.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 |
         +----++----+----+----+
         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:
@@ -3250,9 +3419,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::
@@ -3265,10 +3442,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):
@@ -3470,7 +3681,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:
@@ -3487,30 +3700,45 @@ 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
 
-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