eja: refactor all matrix classes upwards (note: everything broken).

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

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index e9f971f42164b2fc960b82724aa1cf29466c29d9..106a0cddec06355a952e71d75699677b25dd7da9 100644 (file)
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -1,9 +1,56 @@
 """
 """

-Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
-specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
-are used in optimization, and have some additional nice methods beyond
-what can be supported in a general Jordan Algebra.
-
+Representations and constructions for Euclidean Jordan algebras.
+
+A Euclidean Jordan algebra is a Jordan algebra that has some
+additional properties:
+
+  1.   It is finite-dimensional.
+  2.   Its scalar field is the real numbers.
+  3a.  An inner product is defined on it, and...
+  3b.  That inner product is compatible with the Jordan product
+       in the sense that `<x*y,z> = <y,x*z>` for all elements
+       `x,y,z` in the algebra.
+
+Every Euclidean Jordan algebra is formally-real: for any two elements
+`x` and `y` in the algebra, `x^{2} + y^{2} = 0` implies that `x = y =
+0`. Conversely, every finite-dimensional formally-real Jordan algebra
+can be made into a Euclidean Jordan algebra with an appropriate choice
+of inner-product.
+
+Formally-real Jordan algebras were originally studied as a framework
+for quantum mechanics. Today, Euclidean Jordan algebras are crucial in
+symmetric cone optimization, since every symmetric cone arises as the
+cone of squares in some Euclidean Jordan algebra.
+
+It is known that every Euclidean Jordan algebra decomposes into an
+orthogonal direct sum (essentially, a Cartesian product) of simple
+algebras, and that moreover, up to Jordan-algebra isomorphism, there
+are only five families of simple algebras. We provide constructions
+for these simple algebras:
+
+  * :class:`BilinearFormEJA`
+  * :class:`RealSymmetricEJA`
+  * :class:`ComplexHermitianEJA`
+  * :class:`QuaternionHermitianEJA`
+  * :class:`OctonionHermitianEJA`
+
+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. The Albert EJA is simply a special
+case of the :class:`OctonionHermitianEJA` where the matrices are
+three-by-three and the resulting space has dimension 27. And
+last/least, the trivial EJA is exactly what you think it is; it could
+also be obtained by constructing a dimension-zero instance of any of
+the other algebras. Cartesian products of these are also supported
+using the usual ``cartesian_product()`` function; as a result, we
+support (up to isomorphism) all Euclidean Jordan algebras.

 
 SETUP::
 
 
 SETUP::
 


@@ -13,13 +60,11 @@ EXAMPLES::
 
     sage: random_eja()
     Euclidean Jordan algebra of dimension...
 
     sage: random_eja()
     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.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
 from sage.matrix.constructor import matrix
 from sage.matrix.matrix_space import MatrixSpace
 from sage.combinat.free_module import CombinatorialFreeModule
 from sage.matrix.constructor import matrix
 from sage.matrix.matrix_space import MatrixSpace


@@ -29,328 +74,366 @@ from sage.modules.free_module import FreeModule, VectorSpace
 from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
                             PolynomialRing,
                             QuadraticField)
 from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
                             PolynomialRing,
                             QuadraticField)

-from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement
-from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
-from mjo.eja.eja_utils import _mat2vec
+from mjo.eja.eja_element import FiniteDimensionalEJAElement
+from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
+from mjo.eja.eja_utils import _all2list, _mat2vec

 
 

-class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
+class FiniteDimensionalEJA(CombinatorialFreeModule):

     r"""
     r"""

-    The lowest-level class for representing a Euclidean Jordan algebra.
-    """
-    def _coerce_map_from_base_ring(self):
-        """
-        Disable the map from the base ring into the algebra.
+    A finite-dimensional Euclidean Jordan algebra.
+
+    INPUT:
+
+      - ``basis`` -- a tuple; a tuple of basis elements in "matrix
+        form," which must be the same form as the arguments to
+        ``jordan_product`` and ``inner_product``. In reality, "matrix
+        form" can be either vectors, matrices, or a Cartesian product
+        (ordered tuple) of vectors or matrices. All of these would
+        ideally be vector spaces in sage with no special-casing
+        needed; but in reality we turn vectors into column-matrices
+        and Cartesian products `(a,b)` into column matrices
+        `(a,b)^{T}` after converting `a` and `b` themselves.
+
+      - ``jordan_product`` -- a function; afunction of two ``basis``
+        elements (in matrix form) that returns their jordan product,
+        also in matrix form; this will be applied to ``basis`` to
+        compute a multiplication table for the algebra.
+
+      - ``inner_product`` -- a function; a function of two ``basis``
+        elements (in matrix form) that returns their inner
+        product. This will be applied to ``basis`` to compute an
+        inner-product table (basically a matrix) for this algebra.
+
+      - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
+        field for the algebra.
+
+      - ``orthonormalize`` -- boolean (default: ``True``); whether or
+        not to orthonormalize the basis. Doing so is expensive and
+        generally rules out using the rationals as your ``field``, but
+        is required for spectral decompositions.

 
 

-        Performing a nonsense conversion like this automatically
-        is counterpedagogical. The fallback is to try the usual
-        element constructor, which should also fail.
+    SETUP::

 
 

-        SETUP::
+        sage: from mjo.eja.eja_algebra import random_eja

 
 

-            sage: from mjo.eja.eja_algebra import random_eja
+    TESTS:

 
 

-        TESTS::
+    We should compute that an element subalgebra is associative even
+    if we circumvent the element method::

 
 

-            sage: set_random_seed()
-            sage: J = random_eja()
-            sage: J(1)
-            Traceback (most recent call last):
-            ...
-            ValueError: not an element of this algebra
+        sage: set_random_seed()
+        sage: J = random_eja(field=QQ,orthonormalize=False)
+        sage: x = J.random_element()
+        sage: A = x.subalgebra_generated_by(orthonormalize=False)
+        sage: basis = tuple(b.superalgebra_element() for b in A.basis())
+        sage: J.subalgebra(basis, orthonormalize=False).is_associative()
+        True

 
 

-        """
-        return None
+    """
+    Element = FiniteDimensionalEJAElement

 
     def __init__(self,
 
     def __init__(self,

-                 field,
-                 multiplication_table,
-                 inner_product_table,
-                 prefix='e',
-                 category=None,
-                 matrix_basis=None,
+                 basis,
+                 jordan_product,
+                 inner_product,
+                 field=AA,
+                 orthonormalize=True,
+                 associative=None,
+                 cartesian_product=False,

                  check_field=True,
                  check_field=True,

-                 check_axioms=True):
-        """
-        INPUT:
+                 check_axioms=True,
+                 prefix="b"):

 
 

-          * field -- the scalar field for this algebra (must be real)
+        n = len(basis)

 
 

-          * multiplication_table -- the multiplication table for this
-            algebra's implicit basis. Only the lower-triangular portion
-            of the table is used, since the multiplication is assumed
-            to be commutative.
+        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")

 
 

-        SETUP::
+        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")

 
 

-            sage: from mjo.eja.eja_algebra import (
-            ....:   FiniteDimensionalEuclideanJordanAlgebra,
-            ....:   JordanSpinEJA,
-            ....:   random_eja)
+            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")

 
 

-        EXAMPLES:

 
 

-        By definition, Jordan multiplication commutes::
+        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()
+            # method because we need to know the category before we
+            # initialize the algebra.
+            associative = all( jordan_product(jordan_product(bi,bj),bk)
+                               ==
+                               jordan_product(bi,jordan_product(bj,bk))
+                               for bi in basis
+                               for bj in basis
+                               for bk in basis)
+
+        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.
+        CombinatorialFreeModule.__init__(self,
+                                         field,
+                                         range(n),
+                                         prefix=prefix,
+                                         category=category,
+                                         bracket=False)
+
+        # Now comes all of the hard work. We'll be constructing an
+        # 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*b1 + 2*b2.
+        vector_basis = basis

 
 

-            sage: set_random_seed()
-            sage: J = random_eja()
-            sage: x,y = J.random_elements(2)
-            sage: x*y == y*x
-            True
+        degree = 0
+        if n > 0:
+            degree = len(_all2list(basis[0]))

 
 

-        An error is raised if the Jordan product is not commutative::
+        # Build an ambient space that fits our matrix basis when
+        # written out as "long vectors."
+        V = VectorSpace(field, degree)

 
 

-            sage: JP = ((1,2),(0,0))
-            sage: IP = ((1,0),(0,1))
-            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
-            Traceback (most recent call last):
-            ...
-            ValueError: Jordan product is not commutative
+        # The matrix that will hole the orthonormal -> unorthonormal
+        # coordinate transformation.
+        self._deortho_matrix = None

 
 

-        An error is raised if the inner-product is not commutative::
+        if orthonormalize:
+            # Save a copy of the un-orthonormalized basis for later.
+            # Convert it to ambient V (vector) coordinates while we're
+            # at it, because we'd have to do it later anyway.
+            deortho_vector_basis = tuple( V(_all2list(b)) for b in basis )

 
 

-            sage: JP = ((1,0),(0,1))
-            sage: IP = ((1,2),(0,0))
-            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
-            Traceback (most recent call last):
-            ...
-            ValueError: inner-product is not commutative
+            from mjo.eja.eja_utils import gram_schmidt
+            basis = tuple(gram_schmidt(basis, inner_product))

 
 

-        TESTS:
+        # Save the (possibly orthonormalized) matrix basis for
+        # later...
+        self._matrix_basis = basis

 
 

-        The ``field`` we're given must be real with ``check_field=True``::
+        # Now create the vector space for the algebra, which will have
+        # 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)

 
 

-            sage: JordanSpinEJA(2,QQbar)
-            Traceback (most recent call last):
-            ...
-            ValueError: scalar field is not real
+        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.
+            U = V.span_of_basis( deortho_vector_basis, check=check_axioms)
+            self._deortho_matrix = matrix( 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)]
+                                       for i in range(n) ]

 
 

-        The multiplication table must be square with ``check_axioms=True``::
+        # Note: the Jordan and inner-products are defined in terms
+        # of the ambient basis. It's important that their arguments
+        # are in ambient coordinates as well.
+        for i in range(n):
+            for j in range(i+1):
+                # ortho basis w.r.t. ambient coords
+                q_i = basis[i]
+                q_j = basis[j]

 
 

-            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()),((1,),))
-            Traceback (most recent call last):
-            ...
-            ValueError: multiplication table is not square
+                # 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)))
+                self._multiplication_table[i][j] = self.from_vector(elt)
+
+                if not orthonormalize:
+                    # If we're orthonormalizing the basis with respect
+                    # to an inner-product, then the inner-product
+                    # matrix with respect to the resulting basis is
+                    # just going to be the identity.
+                    ip = inner_product(q_i, q_j)
+                    self._inner_product_matrix[i,j] = ip
+                    self._inner_product_matrix[j,i] = ip

 
 

-        The multiplication and inner-product tables must be the same
-        size (and in particular, the inner-product table must also be
-        square) with ``check_axioms=True``::
+        self._inner_product_matrix._cache = {'hermitian': True}
+        self._inner_product_matrix.set_immutable()

 
 

-            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),(()))
-            Traceback (most recent call last):
-            ...
-            ValueError: multiplication and inner-product tables are
-            different sizes
-            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),((1,2),))
+        if check_axioms:
+            if not self._is_jordanian():
+                raise ValueError("Jordan identity does not hold")
+            if not self._inner_product_is_associative():
+                raise ValueError("inner product is not associative")
+
+
+    def _coerce_map_from_base_ring(self):
+        """
+        Disable the map from the base ring into the algebra.
+
+        Performing a nonsense conversion like this automatically
+        is counterpedagogical. The fallback is to try the usual
+        element constructor, which should also fail.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import random_eja
+
+        TESTS::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: J(1)

             Traceback (most recent call last):
             ...
             Traceback (most recent call last):
             ...

-            ValueError: multiplication and inner-product tables are
-            different sizes
+            ValueError: not an element of this algebra

 
         """
 
         """

-        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")
+        return None

 
 
 
 

-        # The multiplication and inner-product tables should be square
-        # if the user wants us to verify them. And we verify them as
-        # soon as possible, because we want to exploit their symmetry.
-        n = len(multiplication_table)
-        if check_axioms:
-            if not all( len(l) == n for l in multiplication_table ):
-                raise ValueError("multiplication table is not square")
-
-            # If the multiplication table is square, we can check if
-            # the inner-product table is square by comparing it to the
-            # multiplication table's dimensions.
-            msg = "multiplication and inner-product tables are different sizes"
-            if not len(inner_product_table) == n:
-                raise ValueError(msg)
-
-            if not all( len(l) == n for l in inner_product_table ):
-                raise ValueError(msg)
-
-            # Check commutativity of the Jordan product (symmetry of
-            # the multiplication table) and the commutativity of the
-            # inner-product (symmetry of the inner-product table)
-            # first if we're going to check them at all.. This has to
-            # be done before we define product_on_basis(), because
-            # that method assumes that self._multiplication_table is
-            # symmetric. And it has to be done before we build
-            # self._inner_product_matrix, because the process used to
-            # construct it assumes symmetry as well.
-            if not all(    multiplication_table[j][i]
-                        == multiplication_table[i][j]
-                        for i in range(n)
-                        for j in range(i+1) ):
-                raise ValueError("Jordan product is not commutative")
+    def product_on_basis(self, i, j):
+        r"""
+        Returns the Jordan product of the `i` and `j`th basis elements.

 
 

-            if not all( inner_product_table[j][i]
-                        == inner_product_table[i][j]
-                        for i in range(n)
-                        for j in range(i+1) ):
-                raise ValueError("inner-product is not commutative")
+        This completely defines the Jordan product on the algebra, and
+        is used direclty by our superclass machinery to implement
+        :meth:`product`.

 
 

-        self._matrix_basis = matrix_basis
-
-        if category is None:
-            category = MagmaticAlgebras(field).FiniteDimensional()
-            category = category.WithBasis().Unital()
-
-        fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
-        fda.__init__(field,
-                     range(n),
-                     prefix=prefix,
-                     category=category)
-        self.print_options(bracket='')
-
-        # The multiplication table we're given is necessarily in terms
-        # of vectors, because we don't have an algebra yet for
-        # anything to be an element of. However, it's faster in the
-        # long run to have the multiplication table be in terms of
-        # algebra elements. We do this after calling the superclass
-        # constructor so that from_vector() knows what to do.
-        #
-        # Note: we take advantage of symmetry here, and only store
-        # the lower-triangular portion of the table.
-        self._multiplication_table = [ [ self.vector_space().zero()
-                                         for j in range(i+1) ]
-                                       for i in range(n) ]
+        SETUP::

 
 

-        for i in range(n):
-            for j in range(i+1):
-                elt = self.from_vector(multiplication_table[i][j])
-                self._multiplication_table[i][j] = elt
-
-        self._multiplication_table = tuple(map(tuple, self._multiplication_table))
-
-        # Save our inner product as a matrix, since the efficiency of
-        # matrix multiplication will usually outweigh the fact that we
-        # have to store a redundant upper- or lower-triangular part.
-        # Pre-cache the fact that these are Hermitian (real symmetric,
-        # in fact) in case some e.g. matrix multiplication routine can
-        # take advantage of it.
-        ip_matrix_constructor = lambda i,j: inner_product_table[i][j] if j <= i else inner_product_table[j][i]
-        self._inner_product_matrix = matrix(field, n, ip_matrix_constructor)
-        self._inner_product_matrix._cache = {'hermitian': True}
-        self._inner_product_matrix.set_immutable()
+            sage: from mjo.eja.eja_algebra import random_eja

 
 

-        if check_axioms:
-            if not self._is_jordanian():
-                raise ValueError("Jordan identity does not hold")
-            if not self._inner_product_is_associative():
-                raise ValueError("inner product is not associative")
+        TESTS::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: n = J.dimension()
+            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)
+            ....:     bi = J.monomial(i)
+            ....:     bj = J.monomial(j)
+            ....:     bi_bj = J.product_on_basis(i,j)
+            sage: bi*bj == bi_bj
+            True

 
 

-    def _element_constructor_(self, elt):

         """
         """

-        Construct an element of this algebra from its vector or matrix
-        representation.
+        # We only stored the lower-triangular portion of the
+        # multiplication table.
+        if j <= i:
+            return self._multiplication_table[i][j]
+        else:
+            return self._multiplication_table[j][i]

 
 

-        This gets called only after the parent element _call_ method
-        fails to find a coercion for the argument.
+    def inner_product(self, x, y):
+        """
+        The inner product associated with this Euclidean Jordan algebra.
+
+        Defaults to the trace inner product, but can be overridden by
+        subclasses if they are sure that the necessary properties are
+        satisfied.

 
         SETUP::
 
 
         SETUP::
 

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

             ....:                                  HadamardEJA,
             ....:                                  HadamardEJA,

-            ....:                                  RealSymmetricEJA)
+            ....:                                  BilinearFormEJA)

 
         EXAMPLES:
 
 
         EXAMPLES:
 

-        The identity in `S^n` is converted to the identity in the EJA::
+        Our inner product is "associative," which means the following for
+        a symmetric bilinear form::

 
 

-            sage: J = RealSymmetricEJA(3)
-            sage: I = matrix.identity(QQ,3)
-            sage: J(I) == J.one()
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: x,y,z = J.random_elements(3)
+            sage: (x*y).inner_product(z) == y.inner_product(x*z)

             True
 
             True
 

-        This skew-symmetric matrix can't be represented in the EJA::
-
-            sage: J = RealSymmetricEJA(3)
-            sage: A = matrix(QQ,3, lambda i,j: i-j)
-            sage: J(A)
-            Traceback (most recent call last):
-            ...
-            ValueError: not an element of this algebra
-

         TESTS:
 
         TESTS:
 

-        Ensure that we can convert any element of the two non-matrix
-        simple algebras (whose matrix representations are columns)
-        back and forth faithfully::
+        Ensure that this is the usual inner product for the algebras
+        over `R^n`::

 
             sage: set_random_seed()
             sage: J = HadamardEJA.random_instance()
 
             sage: set_random_seed()
             sage: J = HadamardEJA.random_instance()

-            sage: x = J.random_element()
-            sage: J(x.to_vector().column()) == x
-            True
-            sage: J = JordanSpinEJA.random_instance()
-            sage: x = J.random_element()
-            sage: J(x.to_vector().column()) == x
+            sage: x,y = J.random_elements(2)
+            sage: actual = x.inner_product(y)
+            sage: expected = x.to_vector().inner_product(y.to_vector())
+            sage: actual == expected

             True
             True

-        """
-        msg = "not an element of this algebra"
-        if elt == 0:
-            # The superclass implementation of random_element()
-            # needs to be able to coerce "0" into the algebra.
-            return self.zero()
-        elif elt in self.base_ring():
-            # Ensure that no base ring -> algebra coercion is performed
-            # by this method. There's some stupidity in sage that would
-            # otherwise propagate to this method; for example, sage thinks
-            # that the integer 3 belongs to the space of 2-by-2 matrices.
-            raise ValueError(msg)

 
 

-        if elt not in self.matrix_space():
-            raise ValueError(msg)
+        Ensure that this is one-half of the trace inner-product in a
+        BilinearFormEJA that isn't just the reals (when ``n`` isn't
+        one). This is in Faraut and Koranyi, and also my "On the
+        symmetry..." paper::

 
 

-        # Thanks for nothing! Matrix spaces aren't vector spaces in
-        # Sage, so we have to figure out its matrix-basis coordinates
-        # ourselves. We use the basis space's ring instead of the
-        # element's ring because the basis space might be an algebraic
-        # closure whereas the base ring of the 3-by-3 identity matrix
-        # could be QQ instead of QQbar.
-        V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols())
-        W = V.span_of_basis( _mat2vec(s) for s in self.matrix_basis() )
+            sage: set_random_seed()
+            sage: J = BilinearFormEJA.random_instance()
+            sage: n = J.dimension()
+            sage: x = J.random_element()
+            sage: y = J.random_element()
+            sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
+            True

 
 

-        try:
-            coords =  W.coordinate_vector(_mat2vec(elt))
-        except ArithmeticError:  # vector is not in free module
-            raise ValueError(msg)
+        """
+        B = self._inner_product_matrix
+        return (B*x.to_vector()).inner_product(y.to_vector())

 
 

-        return self.from_vector(coords)

 
 

-    def _repr_(self):
-        """
-        Return a string representation of ``self``.
+    def is_associative(self):
+        r"""
+        Return whether or not this algebra's Jordan product is associative.

 
         SETUP::
 
 
         SETUP::
 

-            sage: from mjo.eja.eja_algebra import JordanSpinEJA
-
-        TESTS:
+            sage: from mjo.eja.eja_algebra import ComplexHermitianEJA

 
 

-        Ensure that it says what we think it says::
+        EXAMPLES::

 
 

-            sage: JordanSpinEJA(2, field=AA)
-            Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
-            sage: JordanSpinEJA(3, field=RDF)
-            Euclidean Jordan algebra of dimension 3 over Real Double Field
+            sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
+            sage: J.is_associative()
+            False
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
+            sage: A.is_associative()
+            True

 
         """
 
         """

-        fmt = "Euclidean Jordan algebra of dimension {} over {}"
-        return fmt.format(self.dimension(), self.base_ring())
-
-    def product_on_basis(self, i, j):
-        # We only stored the lower-triangular portion of the
-        # multiplication table.
-        if j <= i:
-            return self._multiplication_table[i][j]
-        else:
-            return self._multiplication_table[j][i]
+        return "Associative" in self.category().axioms()

 
     def _is_commutative(self):
         r"""
 
     def _is_commutative(self):
         r"""


@@ -360,9 +443,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         this algebra was constructed with ``check_axioms=False`` and
         passed an invalid multiplication table.
         """
         this algebra was constructed with ``check_axioms=False`` and
         passed an invalid multiplication table.
         """

-        return all( self.product_on_basis(i,j) == self.product_on_basis(i,j)
-                    for i in range(self.dimension())
-                    for j in range(self.dimension()) )
+        return all( x*y == y*x for x in self.gens() for y in self.gens() )

 
     def _is_jordanian(self):
         r"""
 
     def _is_jordanian(self):
         r"""


@@ -381,22 +462,111 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
                     for i in range(self.dimension())
                     for j in range(self.dimension()) )
 
                     for i in range(self.dimension())
                     for j in range(self.dimension()) )
 

-    def _inner_product_is_associative(self):
+    def _jordan_product_is_associative(self):

         r"""
         r"""

-        Return whether or not this algebra's inner product `B` is
-        associative; that is, whether or not `B(xy,z) = B(x,yz)`.
+        Return whether or not this algebra's Jordan product is
+        associative; that is, whether or not `x*(y*z) = (x*y)*z`
+        for all `x,y,x`.

 
 

-        This method should of course always return ``True``, unless
-        this algebra was constructed with ``check_axioms=False`` and
-        passed an invalid multiplication table.
-        """
+        This method should agree with :meth:`is_associative` unless
+        you lied about the value of the ``associative`` parameter
+        when you constructed the algebra.

 
 

-        # Used to check whether or not something is zero in an inexact
-        # ring. This number is sufficient to allow the construction of
-        # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
-        epsilon = 1e-16
+        SETUP::

 
 

-        for i in range(self.dimension()):
+            sage: from mjo.eja.eja_algebra import (random_eja,
+            ....:                                  RealSymmetricEJA,
+            ....:                                  ComplexHermitianEJA,
+            ....:                                  QuaternionHermitianEJA)
+
+        EXAMPLES::
+
+            sage: J = RealSymmetricEJA(4, orthonormalize=False)
+            sage: J._jordan_product_is_associative()
+            False
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by()
+            sage: A._jordan_product_is_associative()
+            True
+
+        ::
+
+            sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: J._jordan_product_is_associative()
+            False
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
+            sage: A._jordan_product_is_associative()
+            True
+
+        ::
+
+            sage: J = QuaternionHermitianEJA(2)
+            sage: J._jordan_product_is_associative()
+            False
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by()
+            sage: A._jordan_product_is_associative()
+            True
+
+        TESTS:
+
+        The values we've presupplied to the constructors agree with
+        the computation::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: J.is_associative() == J._jordan_product_is_associative()
+            True
+
+        """
+        R = self.base_ring()
+
+        # Used to check whether or not something is zero.
+        epsilon = R.zero()
+        if not R.is_exact():
+            # I don't know of any examples that make this magnitude
+            # necessary because I don't know how to make an
+            # associative algebra when the element subalgebra
+            # construction is unreliable (as it is over RDF; we can't
+            # find the degree of an element because we can't compute
+            # the rank of a matrix). But even multiplication of floats
+            # is non-associative, so *some* epsilon is needed... let's
+            # just take the one from _inner_product_is_associative?
+            epsilon = 1e-15
+
+        for i in range(self.dimension()):
+            for j in range(self.dimension()):
+                for k in range(self.dimension()):
+                    x = self.monomial(i)
+                    y = self.monomial(j)
+                    z = self.monomial(k)
+                    diff = (x*y)*z - x*(y*z)
+
+                    if diff.norm() > epsilon:
+                        return False
+
+        return True
+
+    def _inner_product_is_associative(self):
+        r"""
+        Return whether or not this algebra's inner product `B` is
+        associative; that is, whether or not `B(xy,z) = B(x,yz)`.
+
+        This method should of course always return ``True``, unless
+        this algebra was constructed with ``check_axioms=False`` and
+        passed an invalid Jordan or inner-product.
+        """
+        R = self.base_ring()
+
+        # Used to check whether or not something is zero.
+        epsilon = R.zero()
+        if not R.is_exact():
+            # This choice is sufficient to allow the construction of
+            # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
+            epsilon = 1e-15
+
+        for i in range(self.dimension()):

             for j in range(self.dimension()):
                 for k in range(self.dimension()):
                     x = self.monomial(i)
             for j in range(self.dimension()):
                 for k in range(self.dimension()):
                     x = self.monomial(i)


@@ -404,15 +574,145 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
                     z = self.monomial(k)
                     diff = (x*y).inner_product(z) - x.inner_product(y*z)
 
                     z = self.monomial(k)
                     diff = (x*y).inner_product(z) - x.inner_product(y*z)
 

-                    if self.base_ring().is_exact():
-                        if diff != 0:
-                            return False
-                    else:
-                        if diff.abs() > epsilon:
-                            return False
+                    if diff.abs() > epsilon:
+                        return False

 
         return True
 
 
         return True
 

+    def _element_constructor_(self, elt):
+        """
+        Construct an element of this algebra from its vector or matrix
+        representation.
+
+        This gets called only after the parent element _call_ method
+        fails to find a coercion for the argument.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (random_eja,
+            ....:                                  JordanSpinEJA,
+            ....:                                  HadamardEJA,
+            ....:                                  RealSymmetricEJA)
+
+        EXAMPLES:
+
+        The identity in `S^n` is converted to the identity in the EJA::
+
+            sage: J = RealSymmetricEJA(3)
+            sage: I = matrix.identity(QQ,3)
+            sage: J(I) == J.one()
+            True
+
+        This skew-symmetric matrix can't be represented in the EJA::
+
+            sage: J = RealSymmetricEJA(3)
+            sage: A = matrix(QQ,3, lambda i,j: i-j)
+            sage: J(A)
+            Traceback (most recent call last):
+            ...
+            ValueError: not an element of this algebra
+
+        Tuples work as well, provided that the matrix basis for the
+        algebra consists of them::
+
+            sage: J1 = HadamardEJA(3)
+            sage: J2 = RealSymmetricEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
+            b1 + b5
+
+        TESTS:
+
+        Ensure that we can convert any element back and forth
+        faithfully between its matrix and algebra representations::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: x = J.random_element()
+            sage: J(x.to_matrix()) == x
+            True
+
+        We cannot coerce elements between algebras just because their
+        matrix representations are compatible::
+
+            sage: J1 = HadamardEJA(3)
+            sage: J2 = JordanSpinEJA(3)
+            sage: J2(J1.one())
+            Traceback (most recent call last):
+            ...
+            ValueError: not an element of this algebra
+            sage: J1(J2.zero())
+            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():
+            # Ensure that no base ring -> algebra coercion is performed
+            # by this method. There's some stupidity in sage that would
+            # otherwise propagate to this method; for example, sage thinks
+            # 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...
+            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)
+
+        # Thanks for nothing! Matrix spaces aren't vector spaces in
+        # Sage, so we have to figure out its matrix-basis coordinates
+        # ourselves. We use the basis space's ring instead of the
+        # element's ring because the basis space might be an algebraic
+        # closure whereas the base ring of the 3-by-3 identity matrix
+        # could be QQ instead of QQbar.
+        #
+        # And, we also have to handle Cartesian product bases (when
+        # the matrix basis consists of tuples) here. The "good news"
+        # 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)
+
+        try:
+            coords = W.coordinate_vector(V(elt))
+        except ArithmeticError:  # vector is not in free module
+            raise ValueError(msg)
+
+        return self.from_vector(coords)
+
+    def _repr_(self):
+        """
+        Return a string representation of ``self``.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import JordanSpinEJA
+
+        TESTS:
+
+        Ensure that it says what we think it says::
+
+            sage: JordanSpinEJA(2, field=AA)
+            Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
+            sage: JordanSpinEJA(3, field=RDF)
+            Euclidean Jordan algebra of dimension 3 over Real Double Field
+
+        """
+        fmt = "Euclidean Jordan algebra of dimension {} over {}"
+        return fmt.format(self.dimension(), self.base_ring())
+
+

     @cached_method
     def characteristic_polynomial_of(self):
         """
     @cached_method
     def characteristic_polynomial_of(self):
         """


@@ -490,7 +790,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: J = HadamardEJA(2)
             sage: J.coordinate_polynomial_ring()
             Multivariate Polynomial Ring in X1, X2...
             sage: J = HadamardEJA(2)
             sage: J.coordinate_polynomial_ring()
             Multivariate Polynomial Ring in X1, X2...

-            sage: J = RealSymmetricEJA(3,QQ,orthonormalize=False)
+            sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)

             sage: J.coordinate_polynomial_ring()
             Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
 
             sage: J.coordinate_polynomial_ring()
             Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
 


@@ -594,33 +894,28 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: J = JordanSpinEJA(4)
             sage: J.multiplication_table()
             +----++----+----+----+----+
             sage: J = JordanSpinEJA(4)
             sage: J.multiplication_table()
             +----++----+----+----+----+

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

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

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

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

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

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

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

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

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

             +----++----+----+----+----+
 
         """
         n = self.dimension()
             +----++----+----+----+----+
 
         """
         n = self.dimension()

-        M = [ [ self.zero() for j in range(n) ]
-              for i in range(n) ]
-        for i in range(n):
-            for j in range(i+1):
-                M[i][j] = self._multiplication_table[i][j]
-                M[j][i] = M[i][j]
+        # Prepend the header row.
+        M = [["*"] + list(self.gens())]

 
 

-        for i in range(n):
-            # Prepend the left "header" column entry Can't do this in
-            # the loop because it messes up the symmetry.
-            M[i] = [self.monomial(i)] + M[i]
+        # And to each subsequent row, prepend an entry that belongs to
+        # the left-side "header column."
+        M += [ [self.monomial(i)] + [ self.monomial(i)*self.monomial(j)
+                                    for j in range(n) ]
+               for i in range(n) ]

 
 

-        # Prepend the header row.
-        M = [["*"] + list(self.gens())] + M

         return table(M, header_row=True, header_column=True, frame=True)
 
 
         return table(M, header_row=True, header_column=True, frame=True)
 
 


@@ -648,7 +943,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         Why implement this for non-matrix algebras? Avoiding special
         cases for the :class:`BilinearFormEJA` pays with simplicity in
         its own right. But mainly, we would like to be able to assume
         Why implement this for non-matrix algebras? Avoiding special
         cases for the :class:`BilinearFormEJA` pays with simplicity in
         its own right. But mainly, we would like to be able to assume

-        that elements of a :class:`DirectSumEJA` can be displayed
+        that elements of a :class:`CartesianProductEJA` can be displayed

         nicely, without having to have special classes for direct sums
         one of whose components was a matrix algebra.
 
         nicely, without having to have special classes for direct sums
         one of whose components was a matrix algebra.
 


@@ -661,7 +956,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
             sage: J = RealSymmetricEJA(2)
             sage: J.basis()
 
             sage: J = RealSymmetricEJA(2)
             sage: J.basis()

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

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


@@ -672,18 +967,14 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
             sage: J = JordanSpinEJA(2)
             sage: J.basis()
 
             sage: J = JordanSpinEJA(2)
             sage: J.basis()

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

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

-        if self._matrix_basis is None:
-            M = self.matrix_space()
-            return tuple( M(b.to_vector()) for b in self.basis() )
-        else:
-            return self._matrix_basis
+        return self._matrix_basis

 
 
     def matrix_space(self):
 
 
     def matrix_space(self):


@@ -692,19 +983,56 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         we think of them as matrices (including column vectors of the
         appropriate size).
 
         we think of them as matrices (including column vectors of the
         appropriate size).
 

-        Generally this will be an `n`-by-`1` column-vector space,
+        "By default" this will be an `n`-by-`1` column-matrix space,

         except when the algebra is trivial. There it's `n`-by-`n`
         (where `n` is zero), to ensure that two elements of the matrix
         except when the algebra is trivial. There it's `n`-by-`n`
         (where `n` is zero), to ensure that two elements of the matrix

-        space (empty matrices) can be multiplied.
+        space (empty matrices) can be multiplied. For algebras of
+        matrices, this returns the space in which their
+        real embeddings live.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
+            ....:                                  JordanSpinEJA,
+            ....:                                  QuaternionHermitianEJA,
+            ....:                                  TrivialEJA)
+
+        EXAMPLES:
+
+        By default, the matrix representation is just a column-matrix
+        equivalent to the vector representation::
+
+            sage: J = JordanSpinEJA(3)
+            sage: J.matrix_space()
+            Full MatrixSpace of 3 by 1 dense matrices over Algebraic
+            Real Field
+
+        The matrix representation in the trivial algebra is
+        zero-by-zero instead of the usual `n`-by-one::
+
+            sage: J = TrivialEJA()
+            sage: J.matrix_space()
+            Full MatrixSpace of 0 by 0 dense matrices over Algebraic
+            Real Field
+
+        The matrix space for complex/quaternion Hermitian matrix EJA
+        is the space in which their real-embeddings live, not the
+        original complex/quaternion matrix space::
+
+            sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: J.matrix_space()
+            Full MatrixSpace of 4 by 4 dense matrices over Rational Field
+            sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
+            sage: J.matrix_space()
+            Module of 1 by 1 matrices with entries in Quaternion
+            Algebra (-1, -1) with base ring Rational Field over
+            the scalar ring Rational Field

 
 

-        Matrix algebras override this with something more useful.

         """
         if self.is_trivial():
             return MatrixSpace(self.base_ring(), 0)
         """
         if self.is_trivial():
             return MatrixSpace(self.base_ring(), 0)

-        elif self._matrix_basis is None or len(self._matrix_basis) == 0:
-            return MatrixSpace(self.base_ring(), self.dimension(), 1)

         else:
         else:

-            return self._matrix_basis[0].matrix_space()
+            return self.matrix_basis()[0].parent()

 
 
     @cached_method
 
 
     @cached_method


@@ -717,23 +1045,57 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: from mjo.eja.eja_algebra import (HadamardEJA,
             ....:                                  random_eja)
 
             sage: from mjo.eja.eja_algebra import (HadamardEJA,
             ....:                                  random_eja)
 

-        EXAMPLES::
+        EXAMPLES:
+
+        We can compute unit element in the Hadamard EJA::
+
+            sage: J = HadamardEJA(5)
+            sage: J.one()
+            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()
 
             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()
+            c0
+            sage: A.one().superalgebra_element()
+            b0 + b1 + b2 + b3 + b4

 
         TESTS:
 
 
         TESTS:
 

-        The identity element acts like the identity::
+        The identity element acts like the identity, regardless of
+        whether or not we orthonormalize::

 
             sage: set_random_seed()
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: J.one()*x == x and x*J.one() == x
             True
 
             sage: set_random_seed()
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: J.one()*x == x and x*J.one() == x
             True

+            sage: A = x.subalgebra_generated_by()
+            sage: y = A.random_element()
+            sage: A.one()*y == y and y*A.one() == y
+            True
+
+        ::
+
+            sage: set_random_seed()
+            sage: J = random_eja(field=QQ, orthonormalize=False)
+            sage: x = J.random_element()
+            sage: J.one()*x == x and x*J.one() == x
+            True
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
+            sage: y = A.random_element()
+            sage: A.one()*y == y and y*A.one() == y
+            True

 
 

-        The matrix of the unit element's operator is the identity::
+        The matrix of the unit element's operator is the identity,
+        regardless of the base field and whether or not we
+        orthonormalize::

 
             sage: set_random_seed()
             sage: J = random_eja()
 
             sage: set_random_seed()
             sage: J = random_eja()


@@ -741,6 +1103,27 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: expected = matrix.identity(J.base_ring(), J.dimension())
             sage: actual == expected
             True
             sage: expected = matrix.identity(J.base_ring(), J.dimension())
             sage: actual == expected
             True

+            sage: x = J.random_element()
+            sage: A = x.subalgebra_generated_by()
+            sage: actual = A.one().operator().matrix()
+            sage: expected = matrix.identity(A.base_ring(), A.dimension())
+            sage: actual == expected
+            True
+
+        ::
+
+            sage: set_random_seed()
+            sage: J = random_eja(field=QQ, orthonormalize=False)
+            sage: actual = J.one().operator().matrix()
+            sage: expected = matrix.identity(J.base_ring(), J.dimension())
+            sage: actual == expected
+            True
+            sage: x = J.random_element()
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
+            sage: actual = A.one().operator().matrix()
+            sage: expected = matrix.identity(A.base_ring(), A.dimension())
+            sage: actual == expected
+            True

 
         Ensure that the cached unit element (often precomputed by
         hand) agrees with the computed one::
 
         Ensure that the cached unit element (often precomputed by
         hand) agrees with the computed one::


@@ -752,6 +1135,15 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: J.one() == cached
             True
 
             sage: J.one() == cached
             True
 

+        ::
+
+            sage: set_random_seed()
+            sage: J = random_eja(field=QQ, orthonormalize=False)
+            sage: cached = J.one()
+            sage: J.one.clear_cache()
+            sage: J.one() == cached
+            True
+

         """
         # We can brute-force compute the matrices of the operators
         # that correspond to the basis elements of this algebra.
         """
         # We can brute-force compute the matrices of the operators
         # that correspond to the basis elements of this algebra.


@@ -896,14 +1288,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         if not c.is_idempotent():
             raise ValueError("element is not idempotent: %s" % c)
 
         if not c.is_idempotent():
             raise ValueError("element is not idempotent: %s" % c)
 

-        from mjo.eja.eja_subalgebra import FiniteDimensionalEuclideanJordanSubalgebra
-

         # Default these to what they should be if they turn out to be
         # trivial, because eigenspaces_left() won't return eigenvalues
         # 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).
         # Default these to what they should be if they turn out to be
         # trivial, because eigenspaces_left() won't return eigenvalues
         # 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 = FiniteDimensionalEuclideanJordanSubalgebra(self, ())
+        trivial = self.subalgebra(())

         J0 = trivial                          # eigenvalue zero
         J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
         J1 = trivial                          # eigenvalue one
         J0 = trivial                          # eigenvalue zero
         J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
         J1 = trivial                          # eigenvalue one


@@ -913,9 +1303,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
                 J5 = eigspace
             else:
                 gens = tuple( self.from_vector(b) for b in eigspace.basis() )
                 J5 = eigspace
             else:
                 gens = tuple( self.from_vector(b) for b in eigspace.basis() )

-                subalg = FiniteDimensionalEuclideanJordanSubalgebra(self,
-                                                                    gens,
-                                                                    check_axioms=False)
+                subalg = self.subalgebra(gens, check_axioms=False)

                 if eigval == 0:
                     J0 = subalg
                 elif eigval == 1:
                 if eigval == 0:
                     J0 = subalg
                 elif eigval == 1:


@@ -1004,6 +1392,21 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         r"""
         The `r` polynomial coefficients of the "characteristic polynomial
         of" function.
         r"""
         The `r` polynomial coefficients of the "characteristic polynomial
         of" function.

+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import random_eja
+
+        TESTS:
+
+        The theory shows that these are all homogeneous polynomials of
+        a known degree::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
+            True
+

         """
         n = self.dimension()
         R = self.coordinate_polynomial_ring()
         """
         n = self.dimension()
         R = self.coordinate_polynomial_ring()


@@ -1039,10 +1442,17 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         # The theory says that only the first "r" coefficients are
         # nonzero, and they actually live in the original polynomial
 
         # The theory says that only the first "r" coefficients are
         # nonzero, and they actually live in the original polynomial

-        # ring and not the fraction field. We negate them because
-        # in the actual characteristic polynomial, they get moved
-        # to the other side where x^r lives.
-        return -A_rref.solve_right(E*b).change_ring(R)[:r]
+        # ring and not the fraction field. We negate them because in
+        # the actual characteristic polynomial, they get moved to the
+        # other side where x^r lives. We don't bother to trim A_rref
+        # down to a square matrix and solve the resulting system,
+        # because the upper-left r-by-r portion of A_rref is
+        # guaranteed to be the identity matrix, so e.g.
+        #
+        #   A_rref.solve_right(Y)
+        #
+        # would just be returning Y.
+        return (-E*b)[:r].change_ring(R)

 
     @cached_method
     def rank(self):
 
     @cached_method
     def rank(self):


@@ -1103,7 +1513,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
             sage: set_random_seed()    # long time
             sage: J = random_eja()     # long time
 
             sage: set_random_seed()    # long time
             sage: J = random_eja()     # long time

-            sage: caches = J.rank()    # long time
+            sage: cached = J.rank()    # long time

             sage: J.rank.clear_cache() # long time
             sage: J.rank() == cached   # long time
             True
             sage: J.rank.clear_cache() # long time
             sage: J.rank() == cached   # long time
             True


@@ -1112,6 +1522,14 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return len(self._charpoly_coefficients())
 
 
         return len(self._charpoly_coefficients())
 
 

+    def subalgebra(self, basis, **kwargs):
+        r"""
+        Create a subalgebra of this algebra from the given basis.
+        """
+        from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
+        return FiniteDimensionalEJASubalgebra(self, basis, **kwargs)
+
+

     def vector_space(self):
         """
         Return the vector space that underlies this algebra.
     def vector_space(self):
         """
         Return the vector space that underlies this algebra.


@@ -1130,11 +1548,10 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return self.zero().to_vector().parent().ambient_vector_space()
 
 
         return self.zero().to_vector().parent().ambient_vector_space()
 
 

-    Element = FiniteDimensionalEuclideanJordanAlgebraElement

 
 

-class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
+class RationalBasisEJA(FiniteDimensionalEJA):

     r"""
     r"""

-    New class for algebras whose supplied basis elements have all rational entries.
+    Algebras whose supplied basis elements have all rational entries.

 
     SETUP::
 
 
     SETUP::
 


@@ -1155,191 +1572,57 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr
 
     """
     def __init__(self,
 
     """
     def __init__(self,

-                 field,

                  basis,
                  jordan_product,
                  inner_product,
                  basis,
                  jordan_product,
                  inner_product,

-                 orthonormalize=True,
-                 prefix='e',
-                 category=None,
+                 field=AA,

                  check_field=True,
                  check_field=True,

-                 check_axioms=True):
+                 **kwargs):

 
         if check_field:
             # Abuse the check_field parameter to check that the entries of
             # out basis (in ambient coordinates) are in the field QQ.
 
         if check_field:
             # Abuse the check_field parameter to check that the entries of
             # out basis (in ambient coordinates) are in the field QQ.

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

                 raise TypeError("basis not rational")
 
                 raise TypeError("basis not rational")
 

-        n = len(basis)
-        vector_basis = basis
-
-        from sage.structure.element import is_Matrix
-        basis_is_matrices = False
-
-        degree = 0
-        if n > 0:
-            if is_Matrix(basis[0]):
-                basis_is_matrices = True
-                from mjo.eja.eja_utils import _vec2mat
-                vector_basis = tuple( map(_mat2vec,basis) )
-                degree = basis[0].nrows()**2
-            else:
-                degree = basis[0].degree()
-
-        V = VectorSpace(field, degree)
+        super().__init__(basis,
+                         jordan_product,
+                         inner_product,
+                         field=field,
+                         check_field=check_field,
+                         **kwargs)

 
 

-        # If we were asked to orthonormalize, and if the orthonormal
-        # basis is different from the given one, then we also want to
-        # compute multiplication and inner-product tables for the
-        # deorthonormalized basis. These can be used later to
-        # construct a deorthonormalized copy of this algebra over QQ
-        # in which several operations are much faster.

         self._rational_algebra = None
         self._rational_algebra = None

+        if field is not QQ:
+            # There's no point in constructing the extra algebra if this
+            # one is already rational.
+            #
+            # Note: the same Jordan and inner-products work here,
+            # because they are necessarily defined with respect to
+            # ambient coordinates and not any particular basis.
+            self._rational_algebra = FiniteDimensionalEJA(
+                                       basis,
+                                       jordan_product,
+                                       inner_product,
+                                       field=QQ,
+                                       associative=self.is_associative(),
+                                       orthonormalize=False,
+                                       check_field=False,
+                                       check_axioms=False)

 
 

-        if orthonormalize:
-            if self.base_ring() is not QQ:
-                # There's no point in constructing the extra algebra if this
-                # one is already rational. If the original basis is rational
-                # but normalization would make it irrational, then this whole
-                # constructor will just fail anyway as it tries to stick an
-                # irrational number into a rational algebra.
-                #
-                # Note: the same Jordan and inner-products work here,
-                # because they are necessarily defined with respect to
-                # ambient coordinates and not any particular basis.
-                self._rational_algebra = RationalBasisEuclideanJordanAlgebra(
-                                          QQ,
-                                          basis,
-                                          jordan_product,
-                                          inner_product,
-                                          orthonormalize=False,
-                                          prefix=prefix,
-                                          category=category,
-                                          check_field=False,
-                                          check_axioms=False)
-
-            # Compute the deorthonormalized tables before we orthonormalize
-            # the given basis.
-            W = V.span_of_basis( vector_basis )
-
-            # Note: the Jordan and inner-products are defined in terms
-            # of the ambient basis. It's important that their arguments
-            # are in ambient coordinates as well.
-            for i in range(n):
-                for j in range(i+1):
-                    # given basis w.r.t. ambient coords
-                    q_i = vector_basis[i]
-                    q_j = vector_basis[j]
-
-                    if basis_is_matrices:
-                        q_i = _vec2mat(q_i)
-                        q_j = _vec2mat(q_j)
-
-                    elt = jordan_product(q_i, q_j)
-                    ip = inner_product(q_i, q_j)
-
-                    if basis_is_matrices:
-                        # do another mat2vec because the multiplication
-                        # table is in terms of vectors
-                        elt = _mat2vec(elt)
+    @cached_method
+    def _charpoly_coefficients(self):
+        r"""
+        SETUP::

 
 

-        # We overwrite the name "vector_basis" in a second, but never modify it
-        # in place, to this effectively makes a copy of it.
-        deortho_vector_basis = vector_basis
-        self._deortho_matrix = None
+            sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
+            ....:                                  JordanSpinEJA)

 
 

-        if orthonormalize:
-            from mjo.eja.eja_utils import gram_schmidt
-            if basis_is_matrices:
-                vector_ip = lambda x,y: inner_product(_vec2mat(x), _vec2mat(y))
-                vector_basis = gram_schmidt(vector_basis, vector_ip)
-            else:
-                vector_basis = gram_schmidt(vector_basis, inner_product)
-
-            W = V.span_of_basis( vector_basis )
-
-            # Normalize the "matrix" basis, too!
-            basis = vector_basis
-
-            if basis_is_matrices:
-                basis = tuple( map(_vec2mat,basis) )
-
-        W = V.span_of_basis( vector_basis )
-
-        # 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.
-        U = V.span_of_basis( deortho_vector_basis )
-        self._deortho_matrix = matrix( U.coordinate_vector(q)
-                                       for q in vector_basis )
-
-        # If the superclass constructor is going to verify the
-        # symmetry of this table, it has better at least be
-        # square...
-        if check_axioms:
-            mult_table = [ [0 for j in range(n)] for i in range(n) ]
-            ip_table = [ [0 for j in range(n)] for i in range(n) ]
-        else:
-            mult_table = [ [0 for j in range(i+1)] for i in range(n) ]
-            ip_table = [ [0 for j in range(i+1)] for i in range(n) ]
-
-        # Note: the Jordan and inner-products are defined in terms
-        # of the ambient basis. It's important that their arguments
-        # are in ambient coordinates as well.
-        for i in range(n):
-            for j in range(i+1):
-                # ortho basis w.r.t. ambient coords
-                q_i = vector_basis[i]
-                q_j = vector_basis[j]
-
-                if basis_is_matrices:
-                    q_i = _vec2mat(q_i)
-                    q_j = _vec2mat(q_j)
-
-                elt = jordan_product(q_i, q_j)
-                ip = inner_product(q_i, q_j)
-
-                if basis_is_matrices:
-                    # do another mat2vec because the multiplication
-                    # table is in terms of vectors
-                    elt = _mat2vec(elt)
-
-                elt = W.coordinate_vector(elt)
-                mult_table[i][j] = elt
-                ip_table[i][j] = ip
-                if check_axioms:
-                    # The tables are square if we're verifying that they
-                    # are commutative.
-                    mult_table[j][i] = elt
-                    ip_table[j][i] = ip
-
-        if basis_is_matrices:
-            for m in basis:
-                m.set_immutable()
-        else:
-            basis = tuple( x.column() for x in basis )
-
-        super().__init__(field,
-                         mult_table,
-                         ip_table,
-                         prefix,
-                         category,
-                         basis, # matrix basis
-                         check_field,
-                         check_axioms)
-
-    @cached_method
-    def _charpoly_coefficients(self):
-        r"""
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
-            ....:                                  JordanSpinEJA)
-
-        EXAMPLES:
+        EXAMPLES:

 
         The base ring of the resulting polynomial coefficients is what
         it should be, and not the rationals (unless the algebra was
 
         The base ring of the resulting polynomial coefficients is what
         it should be, and not the rationals (unless the algebra was


@@ -1355,13 +1638,12 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr
             Algebraic Real Field
 
         """
             Algebraic Real Field
 
         """

-        if self.base_ring() is QQ or self._rational_algebra is None:
+        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.
             # 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.

-            superclass = super(RationalBasisEuclideanJordanAlgebra, self)
-            return superclass._charpoly_coefficients()
+            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.
 
         # Do the computation over the rationals. The answer will be
         # the same, because all we've done is a change of basis.


@@ -1369,7 +1651,14 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr
         a = ( a_i.change_ring(self.base_ring())
               for a_i in self._rational_algebra._charpoly_coefficients() )
 
         a = ( a_i.change_ring(self.base_ring())
               for a_i in self._rational_algebra._charpoly_coefficients() )
 

-        # Now convert the coordinate variables back to the
+        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)
+
+        # Otherwise, convert the coordinate variables back to the

         # deorthonormalized ones.
         R = self.coordinate_polynomial_ring()
         from sage.modules.free_module_element import vector
         # deorthonormalized ones.
         R = self.coordinate_polynomial_ring()
         from sage.modules.free_module_element import vector


@@ -1379,7 +1668,7 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr
         subs_dict = { X[i]: BX[i] for i in range(len(X)) }
         return tuple( a_i.subs(subs_dict) for a_i in a )
 
         subs_dict = { X[i]: BX[i] for i in range(len(X)) }
         return tuple( a_i.subs(subs_dict) for a_i in a )
 

-class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra):
+class ConcreteEJA(FiniteDimensionalEJA):

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


@@ -1390,7 +1679,7 @@ class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra):
 
     SETUP::
 
 
     SETUP::
 

-        sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
+        sage: from mjo.eja.eja_algebra import ConcreteEJA

 
     TESTS:
 
 
     TESTS:
 


@@ -1398,7 +1687,7 @@ class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra):
     product, unless we specify otherwise::
 
         sage: set_random_seed()
     product, unless we specify otherwise::
 
         sage: set_random_seed()

-        sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
+        sage: J = ConcreteEJA.random_instance()

         sage: all( b.norm() == 1 for b in J.gens() )
         True
 
         sage: all( b.norm() == 1 for b in J.gens() )
         True
 


@@ -1409,7 +1698,7 @@ class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra):
     EJA the operator is self-adjoint by the Jordan axiom::
 
         sage: set_random_seed()
     EJA the operator is self-adjoint by the Jordan axiom::
 
         sage: set_random_seed()

-        sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
+        sage: J = ConcreteEJA.random_instance()

         sage: x = J.random_element()
         sage: x.operator().is_self_adjoint()
         True
         sage: x = J.random_element()
         sage: x.operator().is_self_adjoint()
         True


@@ -1441,77 +1730,155 @@ class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra):
         from sage.misc.prandom import choice
         eja_class = choice(cls.__subclasses__())
 
         from sage.misc.prandom import choice
         eja_class = choice(cls.__subclasses__())
 

-        # These all bubble up to the RationalBasisEuclideanJordanAlgebra
-        # superclass constructor, so any (kw)args valid there are also
-        # valid here.
+        # 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(*args, **kwargs)
 
 

-class MatrixEuclideanJordanAlgebra:
+class MatrixEJA:

     @staticmethod
     @staticmethod

-    def real_embed(M):
+    def _denormalized_basis(A):

         """
         """

-        Embed the matrix ``M`` into a space of real matrices.
+        Returns a basis for the space of complex Hermitian n-by-n matrices.

 
 

-        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.
+        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.

 
 

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

 
 

-        """
-        raise NotImplementedError
+            sage: from mjo.hurwitz import (ComplexMatrixAlgebra,
+            ....:                          QuaternionMatrixAlgebra,
+            ....:                          OctonionMatrixAlgebra)
+            sage: from mjo.eja.eja_algebra import MatrixEJA

 
 

+        TESTS::
+
+            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
+
+        ::
+
+            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

 
 

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

         """
         """

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

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

-    @classmethod
-    def trace_inner_product(cls,X,Y):
-        Xu = cls.real_unembed(X)
-        Yu = cls.real_unembed(Y)
-        tr = (Xu*Yu).trace()
+    @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.

 
 

-        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]
+        SETUP::

 
 

+            sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
+            ....:                                  ComplexHermitianEJA,
+            ....:                                  QuaternionHermitianEJA,
+            ....:                                  OctonionHermitianEJA)

 
 

-class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
-    @staticmethod
-    def real_embed(M):
-        """
-        The identity function, for embedding real matrices into real
-        matrices.
-        """
-        return M
+        EXAMPLES::
+
+            sage: J = RealSymmetricEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2
+
+        ::
+
+            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

 
 

-    @staticmethod
-    def real_unembed(M):
-        """
-        The identity function, for unembedding real matrices from real
-        matrices.

         """
         """

-        return M
+        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(ConcreteEuclideanJordanAlgebra,
-                       RealMatrixEuclideanJordanAlgebra):
+
+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
     """
     The rank-n simple EJA consisting of real symmetric n-by-n
     matrices, the usual symmetric Jordan product, and the trace inner


@@ -1524,19 +1891,19 @@ class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra,
     EXAMPLES::
 
         sage: J = RealSymmetricEJA(2)
     EXAMPLES::
 
         sage: J = RealSymmetricEJA(2)

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

 
     In theory, our "field" can be any subfield of the reals::
 
 
     In theory, our "field" can be any subfield of the reals::
 

-        sage: RealSymmetricEJA(2, RDF)
+        sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)

         Euclidean Jordan algebra of dimension 3 over Real Double Field
         Euclidean Jordan algebra of dimension 3 over Real Double Field

-        sage: RealSymmetricEJA(2, RR)
+        sage: RealSymmetricEJA(2, field=RR, check_axioms=True)

         Euclidean Jordan algebra of dimension 3 over Real Field with
         53 bits of precision
 
         Euclidean Jordan algebra of dimension 3 over Real Field with
         53 bits of precision
 


@@ -1576,221 +1943,39 @@ class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra,
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """

-    @classmethod
-    def _denormalized_basis(cls, n, 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,QQ)
-            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
 
     @classmethod
     @staticmethod
     def _max_random_instance_size():
         return 4 # Dimension 10
 
     @classmethod

-    def random_instance(cls, field=AA, **kwargs):
+    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 a random instance of this type of algebra.
         """
         n = ZZ.random_element(cls._max_random_instance_size() + 1)

-        return cls(n, field, **kwargs)
+        return cls(n, **kwargs)

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

-        basis = self._denormalized_basis(n, field)
-        super(RealSymmetricEJA, self).__init__(field,
-                                               basis,
-                                               self.jordan_product,
-                                               self.trace_inner_product,
-                                               **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
+
+        A = MatrixSpace(field, n)
+        super().__init__(self._denormalized_basis(A),
+                         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)
         self.rank.set_cache(n)

-        self.one.set_cache(self(matrix.identity(field,n)))
-
-
-class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
-    @staticmethod
-    def real_embed(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 \
-            ....:   ComplexMatrixEuclideanJordanAlgebra
-
-        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: ComplexMatrixEuclideanJordanAlgebra.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 = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
-            sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
-            sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
-            sage: Xe*Ye == XYe
-            True
-
-        """
-        n = M.nrows()
-        if M.ncols() != n:
-            raise ValueError("the matrix 'M' must be square")
-
-        # We don't need any adjoined elements...
-        field = M.base_ring().base_ring()
-
-        blocks = []
-        for z in M.list():
-            a = z.list()[0] # real part, I guess
-            b = z.list()[1] # imag part, I guess
-            blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
-
-        return matrix.block(field, n, blocks)
-
-
-    @staticmethod
-    def real_unembed(M):
-        """
-        The inverse of _embed_complex_matrix().
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import \
-            ....:   ComplexMatrixEuclideanJordanAlgebra
-
-        EXAMPLES::
-
-            sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
-            ....:                 [-2,  1,  -4,  3],
-            ....:                 [ 9,  10, 11, 12],
-            ....:                 [-10, 9, -12, 11] ])
-            sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
-            [  2*I + 1   4*I + 3]
-            [ 10*I + 9 12*I + 11]
+        self.one.set_cache(self(A.one()))

 
 

-        TESTS:

 
 

-        Unembedding is the inverse of embedding::

 
 

-            sage: set_random_seed()
-            sage: F = QuadraticField(-1, 'I')
-            sage: M = random_matrix(F, 3)
-            sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
-            sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
-            True
-
-        """
-        n = ZZ(M.nrows())
-        if M.ncols() != n:
-            raise ValueError("the matrix 'M' must be square")
-        if not n.mod(2).is_zero():
-            raise ValueError("the matrix 'M' must be a complex embedding")
-
-        # If "M" was normalized, its base ring might have roots
-        # adjoined and they can stick around after unembedding.
-        field = M.base_ring()
-        R = PolynomialRing(field, 'z')
-        z = R.gen()
-        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
-        else:
-            F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
-        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/2):
-            for j in range(n/2):
-                submat = M[2*k:2*k+2,2*j:2*j+2]
-                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/2, elements)
-
-
-    @classmethod
-    def trace_inner_product(cls,X,Y):
-        """
-        Compute a matrix inner product in this algebra directly from
-        its real embedding.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
-
-        TESTS:
-
-        This gives the same answer as the slow, default method implemented
-        in :class:`MatrixEuclideanJordanAlgebra`::
-
-            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 = ComplexHermitianEJA.real_unembed(Xe)
-            sage: Y = ComplexHermitianEJA.real_unembed(Ye)
-            sage: expected = (X*Y).trace().real()
-            sage: actual = ComplexHermitianEJA.trace_inner_product(Xe,Ye)
-            sage: actual == expected
-            True
-
-        """
-        return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/2
-
-
-class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra,
-                          ComplexMatrixEuclideanJordanAlgebra):
+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,
     """
     The rank-n simple EJA consisting of complex Hermitian n-by-n
     matrices over the real numbers, the usual symmetric Jordan product,


@@ -1805,289 +1990,81 @@ class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra,
 
     In theory, our "field" can be any subfield of the reals::
 
 
     In theory, our "field" can be any subfield of the reals::
 

-        sage: ComplexHermitianEJA(2, RDF)
+        sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True)

         Euclidean Jordan algebra of dimension 4 over Real Double Field
         Euclidean Jordan algebra of dimension 4 over Real Double Field

-        sage: ComplexHermitianEJA(2, RR)
-        Euclidean Jordan algebra of dimension 4 over Real Field with
-        53 bits of precision
-
-    TESTS:
-
-    The dimension of this algebra is `n^2`::
-
-        sage: set_random_seed()
-        sage: n_max = ComplexHermitianEJA._max_random_instance_size()
-        sage: n = ZZ.random_element(1, n_max)
-        sage: J = ComplexHermitianEJA(n)
-        sage: J.dimension() == n^2
-        True
-
-    The Jordan multiplication is what we think it is::
-
-        sage: set_random_seed()
-        sage: J = ComplexHermitianEJA.random_instance()
-        sage: x,y = J.random_elements(2)
-        sage: actual = (x*y).to_matrix()
-        sage: X = x.to_matrix()
-        sage: Y = y.to_matrix()
-        sage: expected = (X*Y + Y*X)/2
-        sage: actual == expected
-        True
-        sage: J(expected) == x*y
-        True
-
-    We can change the generator prefix::
-
-        sage: ComplexHermitianEJA(2, prefix='z').gens()
-        (z0, z1, z2, z3)
-
-    We can construct the (trivial) algebra of rank zero::
-
-        sage: ComplexHermitianEJA(0)
-        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: field = QuadraticField(2, 'sqrt2')
-            sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
-            sage: all( M.is_symmetric() for M in  B)
-            True
-
-        """
-        R = PolynomialRing(field, 'z')
-        z = R.gen()
-        F = field.extension(z**2 + 1, 'I')
-        I = F.gen()
-
-        # 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 = []
-        for i in range(n):
-            for j in range(i+1):
-                Eij = matrix(F, n, lambda k,l: k==i and l==j)
-                if i == j:
-                    Sij = cls.real_embed(Eij)
-                    S.append(Sij)
-                else:
-                    # The second one has a minus because it's conjugated.
-                    Sij_real = cls.real_embed(Eij + Eij.transpose())
-                    S.append(Sij_real)
-                    Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose())
-                    S.append(Sij_imag)
-
-        # Since we embedded these, we can drop back to the "field" that we
-        # started with instead of the complex extension "F".
-        return tuple( s.change_ring(field) for s in S )
-
-
-    def __init__(self, n, field=AA, **kwargs):
-        basis = self._denormalized_basis(n,field)
-        super(ComplexHermitianEJA, self).__init__(field,
-                                               basis,
-                                               self.jordan_product,
-                                               self.trace_inner_product,
-                                               **kwargs)
-        self.rank.set_cache(n)
-        # TODO: pre-cache the identity!
-
-    @staticmethod
-    def _max_random_instance_size():
-        return 3 # Dimension 9
-
-    @classmethod
-    def random_instance(cls, field=AA, **kwargs):
-        """
-        Return a random instance of this type of algebra.
-        """
-        n = ZZ.random_element(cls._max_random_instance_size() + 1)
-        return cls(n, field, **kwargs)
-
-class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
-    @staticmethod
-    def real_embed(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 \
-            ....:   QuaternionMatrixEuclideanJordanAlgebra
-
-        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: QuaternionMatrixEuclideanJordanAlgebra.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 = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
-            sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
-            sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
-            sage: Xe*Ye == XYe
-            True
-
-        """
-        quaternions = M.base_ring()
-        n = M.nrows()
-        if M.ncols() != n:
-            raise ValueError("the matrix 'M' must be square")
-
-        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 = ComplexMatrixEuclideanJordanAlgebra.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)
-
-
-
-    @staticmethod
-    def real_unembed(M):
-        """
-        The inverse of _embed_quaternion_matrix().
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import \
-            ....:   QuaternionMatrixEuclideanJordanAlgebra
-
-        EXAMPLES::
-
-            sage: M = matrix(QQ, [[ 1,  2,  3,  4],
-            ....:                 [-2,  1, -4,  3],
-            ....:                 [-3,  4,  1, -2],
-            ....:                 [-4, -3,  2,  1]])
-            sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
-            [1 + 2*i + 3*j + 4*k]
+        sage: ComplexHermitianEJA(2, field=RR, check_axioms=True)
+        Euclidean Jordan algebra of dimension 4 over Real Field with
+        53 bits of precision

 
 

-        TESTS:
+    TESTS:

 
 

-        Unembedding is the inverse of embedding::
+    The dimension of this algebra is `n^2`::

 
 

-            sage: set_random_seed()
-            sage: Q = QuaternionAlgebra(QQ, -1, -1)
-            sage: M = random_matrix(Q, 3)
-            sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
-            sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
-            True
+        sage: set_random_seed()
+        sage: n_max = ComplexHermitianEJA._max_random_instance_size()
+        sage: n = ZZ.random_element(1, n_max)
+        sage: J = ComplexHermitianEJA(n)
+        sage: J.dimension() == n^2
+        True

 
 

-        """
-        n = ZZ(M.nrows())
-        if M.ncols() != n:
-            raise ValueError("the matrix 'M' must be square")
-        if not n.mod(4).is_zero():
-            raise ValueError("the matrix 'M' must be a quaternion embedding")
-
-        # Use the base ring of the matrix to ensure that its entries can be
-        # multiplied by elements of the quaternion algebra.
-        field = M.base_ring()
-        Q = QuaternionAlgebra(field,-1,-1)
-        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/4):
-            for m in range(n/4):
-                submat = ComplexMatrixEuclideanJordanAlgebra.real_unembed(
-                    M[4*l:4*l+4,4*m:4*m+4] )
-                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/4, elements)
+    The Jordan multiplication is what we think it is::

 
 

+        sage: set_random_seed()
+        sage: J = ComplexHermitianEJA.random_instance()
+        sage: x,y = J.random_elements(2)
+        sage: actual = (x*y).to_matrix()
+        sage: X = x.to_matrix()
+        sage: Y = y.to_matrix()
+        sage: expected = (X*Y + Y*X)/2
+        sage: actual == expected
+        True
+        sage: J(expected) == x*y
+        True

 
 

-    @classmethod
-    def trace_inner_product(cls,X,Y):
-        """
-        Compute a matrix inner product in this algebra directly from
-        its real embedding.
+    We can change the generator prefix::

 
 

-        SETUP::
+        sage: ComplexHermitianEJA(2, prefix='z').gens()
+        (z0, z1, z2, z3)

 
 

-            sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
+    We can construct the (trivial) algebra of rank zero::

 
 

-        TESTS:
+        sage: ComplexHermitianEJA(0)
+        Euclidean Jordan algebra of dimension 0 over Algebraic Real Field

 
 

-        This gives the same answer as the slow, default method implemented
-        in :class:`MatrixEuclideanJordanAlgebra`::
+    """
+    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
+
+        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,
+                         **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)
+        self.one.set_cache(self(A.one()))

 
 

-            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 = QuaternionHermitianEJA.real_unembed(Xe)
-            sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
-            sage: expected = (X*Y).trace().coefficient_tuple()[0]
-            sage: actual = QuaternionHermitianEJA.trace_inner_product(Xe,Ye)
-            sage: actual == expected
-            True
+    @staticmethod
+    def _max_random_instance_size():
+        return 3 # Dimension 9

 
 

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

         """
         """

-        return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/4
+        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        return cls(n, **kwargs)

 
 
 
 

-class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra,
-                             QuaternionMatrixEuclideanJordanAlgebra):
+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
     r"""
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the


@@ -2102,9 +2079,9 @@ class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra,
 
     In theory, our "field" can be any subfield of the reals::
 
 
     In theory, our "field" can be any subfield of the reals::
 

-        sage: QuaternionHermitianEJA(2, RDF)
+        sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)

         Euclidean Jordan algebra of dimension 6 over Real Double Field
         Euclidean Jordan algebra of dimension 6 over Real Double Field

-        sage: QuaternionHermitianEJA(2, RR)
+        sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)

         Euclidean Jordan algebra of dimension 6 over Real Field with
         53 bits of precision
 
         Euclidean Jordan algebra of dimension 6 over Real Field with
         53 bits of precision
 


@@ -2144,96 +2121,195 @@ class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra,
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """

-    @classmethod
-    def _denormalized_basis(cls, n, field):
-        """
-        Returns a basis for the space of quaternion Hermitian n-by-n matrices.
+    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
+
+        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,
+                         **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)
+        self.one.set_cache(self(A.one()))

 
 

-        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::
+    @staticmethod
+    def _max_random_instance_size():
+        r"""
+        The maximum rank of a random QuaternionHermitianEJA.
+        """
+        return 2 # Dimension 6

 
 

-            sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
+    @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)

 
 

-        TESTS::
+class OctonionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
+    r"""
+    SETUP::

 
 

-            sage: set_random_seed()
-            sage: n = ZZ.random_element(1,5)
-            sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
-            sage: all( M.is_symmetric() for M in B )
-            True
+        sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
+        ....:                                  OctonionHermitianEJA)

 
 

-        """
-        Q = QuaternionAlgebra(QQ,-1,-1)
-        I,J,K = Q.gens()
+    EXAMPLES:

 
 

-        # 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 = []
-        for i in range(n):
-            for j in range(i+1):
-                Eij = matrix(Q, n, lambda k,l: k==i and l==j)
-                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.
-                    Sij_real = cls.real_embed(Eij + Eij.transpose())
-                    S.append(Sij_real)
-                    Sij_I = cls.real_embed(I*Eij - I*Eij.transpose())
-                    S.append(Sij_I)
-                    Sij_J = cls.real_embed(J*Eij - J*Eij.transpose())
-                    S.append(Sij_J)
-                    Sij_K = cls.real_embed(K*Eij - K*Eij.transpose())
-                    S.append(Sij_K)
-
-        # Since we embedded these, we can drop back to the "field" that we
-        # started with instead of the quaternion algebra "Q".
-        return tuple( s.change_ring(field) for s in S )
+    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:

 
 

-    def __init__(self, n, field=AA, **kwargs):
-        basis = self._denormalized_basis(n,field)
-        super(QuaternionHermitianEJA, self).__init__(field,
-                                                     basis,
-                                                     self.jordan_product,
-                                                     self.trace_inner_product,
-                                                     **kwargs)
-        self.rank.set_cache(n)
-        # TODO: cache one()!
+    We can actually construct the 27-dimensional Albert algebra,
+    and we get the right unit element if we recompute it::
+
+        sage: J = OctonionHermitianEJA(3,                    # long time
+        ....:                          field=QQ,             # long time
+        ....:                          orthonormalize=False) # long time
+        sage: J.one.clear_cache()                            # long time
+        sage: J.one()                                        # long time
+        b0 + b9 + b26
+        sage: J.one().to_matrix()                            # long time
+        +----+----+----+
+        | e0 | 0  | 0  |
+        +----+----+----+
+        | 0  | e0 | 0  |
+        +----+----+----+
+        | 0  | 0  | e0 |
+        +----+----+----+
+
+    The 2-by-2 algebra is isomorphic to the ten-dimensional Jordan
+    spin algebra, but just to be sure, we recompute its rank::
+
+        sage: J = OctonionHermitianEJA(2,                    # long time
+        ....:                          field=QQ,             # long time
+        ....:                          orthonormalize=False) # long time
+        sage: J.rank.clear_cache()                           # long time
+        sage: J.rank()                                       # long time
+        2

 
 

+    """

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

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

 
     @classmethod
 
     @classmethod

-    def random_instance(cls, field=AA, **kwargs):
+    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 a random instance of this type of algebra.
         """
         n = ZZ.random_element(cls._max_random_instance_size() + 1)

-        return cls(n, field, **kwargs)
+        return cls(n, **kwargs)
+
+    def __init__(self, n, field=AA, **kwargs):
+        if n > 3:
+            # Otherwise we don't get an EJA.
+            raise ValueError("n cannot exceed 3")
+
+        # We know this is a valid EJA, but will double-check
+        # if the user passes check_axioms=True.
+        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
+
+        from mjo.hurwitz import OctonionMatrixAlgebra
+        A = OctonionMatrixAlgebra(n, scalars=field)
+        super().__init__(self._denormalized_basis(A),
+                         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)
+        self.one.set_cache(self(A.one()))
+
+
+class AlbertEJA(OctonionHermitianEJA):
+    r"""
+    The Albert algebra is the algebra of three-by-three Hermitian
+    matrices whose entries are octonions.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import AlbertEJA
+
+    EXAMPLES::

 
 

+        sage: AlbertEJA(field=QQ, orthonormalize=False)
+        Euclidean Jordan algebra of dimension 27 over Rational Field
+        sage: AlbertEJA() # long time
+        Euclidean Jordan algebra of dimension 27 over Algebraic Real Field

 
 

-class HadamardEJA(ConcreteEuclideanJordanAlgebra):

     """
     """

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

 
 

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

 
     SETUP::
 
 
     SETUP::
 


@@ -2244,19 +2320,19 @@ class HadamardEJA(ConcreteEuclideanJordanAlgebra):
     This multiplication table can be verified by hand::
 
         sage: J = HadamardEJA(3)
     This multiplication table can be verified by hand::
 
         sage: J = HadamardEJA(3)

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

         0
         0

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

         0
         0

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

         0
         0

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

 
     TESTS:
 
 
     TESTS:
 


@@ -2264,27 +2340,36 @@ class HadamardEJA(ConcreteEuclideanJordanAlgebra):
 
         sage: HadamardEJA(3, prefix='r').gens()
         (r0, r1, r2)
 
         sage: HadamardEJA(3, prefix='r').gens()
         (r0, r1, r2)

-

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

-        V = VectorSpace(field, n)
-        basis = V.basis()
-
-        def jordan_product(x,y):
-            return V([ xi*yi for (xi,yi) in zip(x,y) ])
-        def inner_product(x,y):
-            return x.inner_product(y)
-
-        # Don't orthonormalize because our basis is already orthonormal
-        # with respect to our inner-product.
-        super(HadamardEJA, self).__init__(field,
-                                          basis,
-                                          jordan_product,
-                                          inner_product,
-                                          orthonormalize=False,
-                                          check_field=False,
-                                          check_axioms=False,
-                                          **kwargs)
+        if n == 0:
+            jordan_product = lambda x,y: x
+            inner_product = lambda x,y: x
+        else:
+            def jordan_product(x,y):
+                P = x.parent()
+                return P( xi*yi for (xi,yi) in zip(x,y) )
+
+            def inner_product(x,y):
+                return (x.T*y)[0,0]
+
+        # New defaults for keyword arguments. Don't orthonormalize
+        # because our basis is already orthonormal with respect to our
+        # inner-product. Don't check the axioms, because we know this
+        # is a valid EJA... but do double-check if the user passes
+        # check_axioms=True. Note: we DON'T override the "check_field"
+        # default here, because the user can pass in a field!
+        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(field, n).basis() )
+        super().__init__(column_basis,
+                         jordan_product,
+                         inner_product,
+                         field=field,
+                         associative=True,
+                         **kwargs)

         self.rank.set_cache(n)
 
         if n == 0:
         self.rank.set_cache(n)
 
         if n == 0:


@@ -2300,15 +2385,15 @@ class HadamardEJA(ConcreteEuclideanJordanAlgebra):
         return 5
 
     @classmethod
         return 5
 
     @classmethod

-    def random_instance(cls, field=AA, **kwargs):
+    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 a random instance of this type of algebra.
         """
         n = ZZ.random_element(cls._max_random_instance_size() + 1)

-        return cls(n, field, **kwargs)
+        return cls(n, **kwargs)

 
 
 
 

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

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


@@ -2388,31 +2473,50 @@ class BilinearFormEJA(ConcreteEuclideanJordanAlgebra):
         ....:              for j in range(n-1) ]
         sage: actual == expected
         True
         ....:              for j in range(n-1) ]
         sage: actual == expected
         True

+

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

-        if not B.is_positive_definite():
-            raise ValueError("bilinear form is not positive-definite")
-
-        n = B.nrows()
-        V = VectorSpace(field, n)
+        # 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...
+        if ("check_axioms" not in kwargs) or kwargs["check_axioms"]:
+            if not B.is_positive_definite():
+                raise ValueError("bilinear form is not positive-definite")
+
+        # However, all of the other data for this EJA is computed
+        # by us in manner that guarantees the axioms are
+        # satisfied. So, again, unless we are specifically asked to
+        # verify things, we'll skip the rest of the checks.
+        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False

 
         def inner_product(x,y):
 
         def inner_product(x,y):

-            return (B*x).inner_product(y)
+            return (y.T*B*x)[0,0]

 
         def jordan_product(x,y):
 
         def jordan_product(x,y):

-            x0 = x[0]
-            xbar = x[1:]
-            y0 = y[0]
-            ybar = y[1:]
-            z0 = inner_product(x,y)
+            P = x.parent()
+            x0 = x[0,0]
+            xbar = x[1:,0]
+            y0 = y[0,0]
+            ybar = y[1:,0]
+            z0 = inner_product(y,x)

             zbar = y0*xbar + x0*ybar
             zbar = y0*xbar + x0*ybar

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

 
 

-        super(BilinearFormEJA, self).__init__(field,
-                                              V.basis(),
-                                              jordan_product,
-                                              inner_product,
-                                              **kwargs)
+        n = B.nrows()
+        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
+        if n <= 2:
+            associative = True
+
+        super().__init__(column_basis,
+                         jordan_product,
+                         inner_product,
+                         field=field,
+                         associative=associative,
+                         **kwargs)

 
         # The rank of this algebra is two, unless we're in a
         # one-dimensional ambient space (because the rank is bounded
 
         # The rank of this algebra is two, unless we're in a
         # one-dimensional ambient space (because the rank is bounded


@@ -2432,28 +2536,28 @@ class BilinearFormEJA(ConcreteEuclideanJordanAlgebra):
         return 5
 
     @classmethod
         return 5
 
     @classmethod

-    def random_instance(cls, field=AA, **kwargs):
+    def random_instance(cls, **kwargs):

         """
         Return a random instance of this algebra.
         """
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
         if n.is_zero():
         """
         Return a random instance of this algebra.
         """
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
         if n.is_zero():

-            B = matrix.identity(field, n)
-            return cls(B, field, **kwargs)
+            B = matrix.identity(ZZ, n)
+            return cls(B, **kwargs)

 
 

-        B11 = matrix.identity(field,1)
-        M = matrix.random(field, n-1)
-        I = matrix.identity(field, n-1)
-        alpha = field.zero()
+        B11 = matrix.identity(ZZ, 1)
+        M = matrix.random(ZZ, n-1)
+        I = matrix.identity(ZZ, n-1)
+        alpha = ZZ.zero()

         while alpha.is_zero():
         while alpha.is_zero():

-            alpha = field.random_element().abs()
+            alpha = ZZ.random_element().abs()

         B22 = M.transpose()*M + alpha*I
 
         from sage.matrix.special import block_matrix
         B = block_matrix(2,2, [ [B11,   ZZ(0) ],
                                 [ZZ(0), B22 ] ])
 
         B22 = M.transpose()*M + alpha*I
 
         from sage.matrix.special import block_matrix
         B = block_matrix(2,2, [ [B11,   ZZ(0) ],
                                 [ZZ(0), B22 ] ])
 

-        return cls(B, field, **kwargs)
+        return cls(B, **kwargs)

 
 
 class JordanSpinEJA(BilinearFormEJA):
 
 
 class JordanSpinEJA(BilinearFormEJA):


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

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

         0
         0

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

         0
         0

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

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


@@ -2506,19 +2610,18 @@ class JordanSpinEJA(BilinearFormEJA):
             True
 
     """
             True
 
     """

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

         # This is a special case of the BilinearFormEJA with the
         # identity matrix as its bilinear form.
         # This is a special case of the BilinearFormEJA with the
         # identity matrix as its bilinear form.

-        B = matrix.identity(field, n)
+        B = matrix.identity(ZZ, n)

 
         # Don't orthonormalize because our basis is already
         # orthonormal with respect to our inner-product.
 
         # Don't orthonormalize because our basis is already
         # orthonormal with respect to our inner-product.

-        super(JordanSpinEJA, self).__init__(B,
-                                            field,
-                                            orthonormalize=False,
-                                            check_field=False,
-                                            check_axioms=False,
-                                            **kwargs)
+        if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
+
+        # But also don't pass check_field=False here, because the user
+        # can pass in a field!
+        super().__init__(B, *args, **kwargs)

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


@@ -2528,17 +2631,17 @@ class JordanSpinEJA(BilinearFormEJA):
         return 5
 
     @classmethod
         return 5
 
     @classmethod

-    def random_instance(cls, field=AA, **kwargs):
+    def random_instance(cls, **kwargs):

         """
         Return a random instance of this type of algebra.
 
         Needed here to override the implementation for ``BilinearFormEJA``.
         """
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
         """
         Return a random instance of this type of algebra.
 
         Needed here to override the implementation for ``BilinearFormEJA``.
         """
         n = ZZ.random_element(cls._max_random_instance_size() + 1)

-        return cls(n, field, **kwargs)
+        return cls(n, **kwargs)

 
 
 
 

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

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


@@ -2567,175 +2670,455 @@ class TrivialEJA(ConcreteEuclideanJordanAlgebra):
         0
 
     """
         0
 
     """

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

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

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

         basis = ()
         basis = ()

-        super(TrivialEJA, self).__init__(field,
-                                         basis,
-                                         jordan_product,
-                                         inner_product,
-                                         **kwargs)
+
+        # New defaults for keyword arguments
+        if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
+        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
+
+        super().__init__(basis,
+                         jordan_product,
+                         inner_product,
+                         associative=True,
+                         **kwargs)
+

         # The rank is zero using my definition, namely the dimension of the
         # largest subalgebra generated by any element.
         self.rank.set_cache(0)
         self.one.set_cache( self.zero() )
 
     @classmethod
         # The rank is zero using my definition, namely the dimension of the
         # largest subalgebra generated by any element.
         self.rank.set_cache(0)
         self.one.set_cache( self.zero() )
 
     @classmethod

-    def random_instance(cls, field=AA, **kwargs):
+    def random_instance(cls, **kwargs):

         # We don't take a "size" argument so the superclass method is
         # inappropriate for us.
         # We don't take a "size" argument so the superclass method is
         # inappropriate for us.

-        return cls(field, **kwargs)
+        return cls(**kwargs)

 
 

-class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
+
+class CartesianProductEJA(FiniteDimensionalEJA):

     r"""
     r"""

-    The external (orthogonal) direct sum of two other Euclidean Jordan
-    algebras. Essentially the Cartesian product of its two factors.
-    Every Euclidean Jordan algebra decomposes into an orthogonal
-    direct sum of simple Euclidean Jordan algebras, so no generality
-    is lost by providing only this construction.
+    The external (orthogonal) direct sum of two or more Euclidean
+    Jordan algebras. Every Euclidean Jordan algebra decomposes into an
+    orthogonal direct sum of simple Euclidean Jordan algebras which is
+    then isometric to a Cartesian product, so no generality is lost by
+    providing only this construction.

 
     SETUP::
 
         sage: from mjo.eja.eja_algebra import (random_eja,
 
     SETUP::
 
         sage: from mjo.eja.eja_algebra import (random_eja,

+        ....:                                  CartesianProductEJA,

         ....:                                  HadamardEJA,
         ....:                                  HadamardEJA,

-        ....:                                  RealSymmetricEJA,
-        ....:                                  DirectSumEJA)
+        ....:                                  JordanSpinEJA,
+        ....:                                  RealSymmetricEJA)

 
 

-    EXAMPLES::
+    EXAMPLES:
+
+    The Jordan product is inherited from our factors and implemented by
+    our CombinatorialFreeModule Cartesian product superclass::

 
 

+        sage: set_random_seed()

         sage: J1 = HadamardEJA(2)
         sage: J1 = HadamardEJA(2)

-        sage: J2 = RealSymmetricEJA(3)
-        sage: J = DirectSumEJA(J1,J2)
-        sage: J.dimension()
-        8
+        sage: J2 = RealSymmetricEJA(2)
+        sage: J = cartesian_product([J1,J2])
+        sage: x,y = J.random_elements(2)
+        sage: x*y in J
+        True
+
+    The ability to retrieve the original factors is implemented by our
+    CombinatorialFreeModule Cartesian product superclass::
+
+        sage: J1 = HadamardEJA(2, field=QQ)
+        sage: J2 = JordanSpinEJA(3, field=QQ)
+        sage: J = cartesian_product([J1,J2])
+        sage: J.cartesian_factors()
+        (Euclidean Jordan algebra of dimension 2 over Rational Field,
+         Euclidean Jordan algebra of dimension 3 over Rational Field)
+
+    You can provide more than two factors::
+
+        sage: J1 = HadamardEJA(2)
+        sage: J2 = JordanSpinEJA(3)
+        sage: J3 = RealSymmetricEJA(3)
+        sage: cartesian_product([J1,J2,J3])
+        Euclidean Jordan algebra of dimension 2 over Algebraic Real
+        Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
+        Real Field (+) Euclidean Jordan algebra of dimension 6 over
+        Algebraic Real Field
+
+    Rank is additive on a Cartesian product::
+
+        sage: J1 = HadamardEJA(1)
+        sage: J2 = RealSymmetricEJA(2)
+        sage: J = cartesian_product([J1,J2])
+        sage: J1.rank.clear_cache()
+        sage: J2.rank.clear_cache()
+        sage: J.rank.clear_cache()

         sage: J.rank()
         sage: J.rank()

-        5
+        3
+        sage: J.rank() == J1.rank() + J2.rank()
+        True
+
+    The same rank computation works over the rationals, with whatever
+    basis you like::
+
+        sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
+        sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
+        sage: J = cartesian_product([J1,J2])
+        sage: J1.rank.clear_cache()
+        sage: J2.rank.clear_cache()
+        sage: J.rank.clear_cache()
+        sage: J.rank()
+        3
+        sage: J.rank() == J1.rank() + J2.rank()
+        True
+
+    The product algebra will be associative if and only if all of its
+    components are associative::
+
+        sage: J1 = HadamardEJA(2)
+        sage: J1.is_associative()
+        True
+        sage: J2 = HadamardEJA(3)
+        sage: J2.is_associative()
+        True
+        sage: J3 = RealSymmetricEJA(3)
+        sage: J3.is_associative()
+        False
+        sage: CP1 = cartesian_product([J1,J2])
+        sage: CP1.is_associative()
+        True
+        sage: CP2 = cartesian_product([J1,J3])
+        sage: CP2.is_associative()
+        False
+
+    Cartesian products of Cartesian products work::
+
+        sage: J1 = JordanSpinEJA(1)
+        sage: J2 = JordanSpinEJA(1)
+        sage: J3 = JordanSpinEJA(1)
+        sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
+        sage: J.multiplication_table()
+        +----++----+----+----+
+        | *  || b0 | b1 | b2 |
+        +====++====+====+====+
+        | b0 || b0 | 0  | 0  |
+        +----++----+----+----+
+        | b1 || 0  | b1 | 0  |
+        +----++----+----+----+
+        | b2 || 0  | 0  | b2 |
+        +----++----+----+----+
+        sage: HadamardEJA(3).multiplication_table()
+        +----++----+----+----+
+        | *  || b0 | b1 | b2 |
+        +====++====+====+====+
+        | b0 || b0 | 0  | 0  |
+        +----++----+----+----+
+        | b1 || 0  | b1 | 0  |
+        +----++----+----+----+
+        | b2 || 0  | 0  | b2 |
+        +----++----+----+----+

 
     TESTS:
 
 
     TESTS:
 

-    The external direct sum construction is only valid when the two factors
-    have the same base ring; an error is raised otherwise::
+    All factors must share the same base field::

 
 

-        sage: set_random_seed()
-        sage: J1 = random_eja(AA)
-        sage: J2 = random_eja(QQ,orthonormalize=False)
-        sage: J = DirectSumEJA(J1,J2)
+        sage: J1 = HadamardEJA(2, field=QQ)
+        sage: J2 = RealSymmetricEJA(2)
+        sage: CartesianProductEJA((J1,J2))

         Traceback (most recent call last):
         ...
         Traceback (most recent call last):
         ...

-        ValueError: algebras must share the same base field
+        ValueError: all factors must share the same base field
+
+    The cached unit element is the same one that would be computed::
+
+        sage: set_random_seed()              # long time
+        sage: J1 = random_eja()              # long time
+        sage: J2 = random_eja()              # long time
+        sage: J = cartesian_product([J1,J2]) # long time
+        sage: actual = J.one()               # long time
+        sage: J.one.clear_cache()            # long time
+        sage: expected = J.one()             # long time
+        sage: actual == expected             # long time
+        True

 
     """
 
     """

-    def __init__(self, J1, J2, **kwargs):
-        if J1.base_ring() != J2.base_ring():
-            raise ValueError("algebras must share the same base field")
-        field = J1.base_ring()
-
-        self._factors = (J1, J2)
-        n1 = J1.dimension()
-        n2 = J2.dimension()
-        n = n1+n2
-        V = VectorSpace(field, n)
-        mult_table = [ [ V.zero() for j in range(i+1) ]
-                       for i in range(n) ]
-        for i in range(n1):
-            for j in range(i+1):
-                p = (J1.monomial(i)*J1.monomial(j)).to_vector()
-                mult_table[i][j] = V(p.list() + [field.zero()]*n2)
+    Element = FiniteDimensionalEJAElement

 
 

-        for i in range(n2):
-            for j in range(i+1):
-                p = (J2.monomial(i)*J2.monomial(j)).to_vector()
-                mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
-
-        # TODO: build the IP table here from the two constituent IP
-        # matrices (it'll be block diagonal, I think).
-        ip_table = [ [ field.zero() for j in range(i+1) ]
-                       for i in range(n) ]
-        super(DirectSumEJA, self).__init__(field,
-                                           mult_table,
-                                           ip_table,
-                                           check_axioms=False,
-                                           **kwargs)
-        self.rank.set_cache(J1.rank() + J2.rank())
-
-
-    def factors(self):
+
+    def __init__(self, factors, **kwargs):
+        m = len(factors)
+        if m == 0:
+            return TrivialEJA()
+
+        self._sets = factors
+
+        field = factors[0].base_ring()
+        if not all( J.base_ring() == field for J in factors ):
+            raise ValueError("all factors must share the same base field")
+
+        associative = all( f.is_associative() for f in factors )
+
+        MS = self.matrix_space()
+        basis = []
+        zero = MS.zero()
+        for i in range(m):
+            for b in factors[i].matrix_basis():
+                z = list(zero)
+                z[i] = b
+                basis.append(z)
+
+        basis = tuple( MS(b) for b in basis )
+
+        # Define jordan/inner products that operate on that matrix_basis.
+        def jordan_product(x,y):
+            return MS(tuple(
+                (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)
+            )
+
+        # 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,
+                                      orthonormalize=False,
+                                      associative=associative,
+                                      cartesian_product=True,
+                                      check_field=False,
+                                      check_axioms=False)
+
+        ones = tuple(J.one().to_matrix() for J in factors)
+        self.one.set_cache(self(ones))
+        self.rank.set_cache(sum(J.rank() for J in factors))
+
+    def cartesian_factors(self):
+        # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
+        return self._sets
+
+    def cartesian_factor(self, i):
+        r"""
+        Return the ``i``th factor of this algebra.
+        """
+        return self._sets[i]
+
+    def _repr_(self):
+        # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
+        from sage.categories.cartesian_product import cartesian_product
+        return cartesian_product.symbol.join("%s" % factor
+                                             for factor in self._sets)
+
+    def matrix_space(self):

         r"""
         r"""

-        Return the pair of this algebra's factors.
+        Return the space that our matrix basis lives in as a Cartesian
+        product.
+
+        We don't simply use the ``cartesian_product()`` functor here
+        because it acts differently on SageMath MatrixSpaces and our
+        custom MatrixAlgebras, which are CombinatorialFreeModules. We
+        always want the result to be represented (and indexed) as
+        an ordered tuple.

 
         SETUP::
 
 
         SETUP::
 

-            sage: from mjo.eja.eja_algebra import (HadamardEJA,
-            ....:                                  JordanSpinEJA,
-            ....:                                  DirectSumEJA)
+            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
+            ....:                                  HadamardEJA,
+            ....:                                  OctonionHermitianEJA,
+            ....:                                  RealSymmetricEJA)

 
         EXAMPLES::
 
 
         EXAMPLES::
 

-            sage: J1 = HadamardEJA(2,QQ)
-            sage: J2 = JordanSpinEJA(3,QQ)
-            sage: J = DirectSumEJA(J1,J2)
-            sage: J.factors()
-            (Euclidean Jordan algebra of dimension 2 over Rational Field,
-             Euclidean Jordan algebra of dimension 3 over Rational Field)
+            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 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 |
+            +----+

 
         """
 
         """

-        return self._factors
+        scalars = self.cartesian_factor(0).base_ring()

 
 

-    def projections(self):
-        r"""
-        Return a pair of projections onto this algebra's factors.
+        # 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):
+        r"""

         SETUP::
 
         SETUP::
 

-            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
-            ....:                                  ComplexHermitianEJA,
-            ....:                                  DirectSumEJA)
+            sage: from mjo.eja.eja_algebra import (random_eja,
+            ....:                                  JordanSpinEJA,
+            ....:                                  HadamardEJA,
+            ....:                                  RealSymmetricEJA,
+            ....:                                  ComplexHermitianEJA)

 
 

-        EXAMPLES::
+        EXAMPLES:
+
+        The projection morphisms are Euclidean Jordan algebra
+        operators::
+
+            sage: J1 = HadamardEJA(2)
+            sage: J2 = RealSymmetricEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: J.cartesian_projection(0)
+            Linear operator between finite-dimensional Euclidean Jordan
+            algebras represented by the matrix:
+            [1 0 0 0 0]
+            [0 1 0 0 0]
+            Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
+            Real Field (+) Euclidean Jordan algebra of dimension 3 over
+            Algebraic Real Field
+            Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
+            Real Field
+            sage: J.cartesian_projection(1)
+            Linear operator between finite-dimensional Euclidean Jordan
+            algebras represented by the matrix:
+            [0 0 1 0 0]
+            [0 0 0 1 0]
+            [0 0 0 0 1]
+            Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
+            Real Field (+) Euclidean Jordan algebra of dimension 3 over
+            Algebraic Real Field
+            Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
+            Real Field
+
+        The projections work the way you'd expect on the vector
+        representation of an element::

 
             sage: J1 = JordanSpinEJA(2)
             sage: J2 = ComplexHermitianEJA(2)
 
             sage: J1 = JordanSpinEJA(2)
             sage: J2 = ComplexHermitianEJA(2)

-            sage: J = DirectSumEJA(J1,J2)
-            sage: (pi_left, pi_right) = J.projections()
-            sage: J.one().to_vector()
-            (1, 0, 1, 0, 0, 1)
+            sage: J = cartesian_product([J1,J2])
+            sage: pi_left = J.cartesian_projection(0)
+            sage: pi_right = J.cartesian_projection(1)

             sage: pi_left(J.one()).to_vector()
             (1, 0)
             sage: pi_right(J.one()).to_vector()
             (1, 0, 0, 1)
             sage: pi_left(J.one()).to_vector()
             (1, 0)
             sage: pi_right(J.one()).to_vector()
             (1, 0, 0, 1)

+            sage: J.one().to_vector()
+            (1, 0, 1, 0, 0, 1)
+
+        TESTS:
+
+        The answer never changes::
+
+            sage: set_random_seed()
+            sage: J1 = random_eja()
+            sage: J2 = random_eja()
+            sage: J = cartesian_product([J1,J2])
+            sage: P0 = J.cartesian_projection(0)
+            sage: P1 = J.cartesian_projection(0)
+            sage: P0 == P1
+            True

 
         """
 
         """

-        (J1,J2) = self.factors()
-        m = J1.dimension()
-        n = J2.dimension()
-        V_basis = self.vector_space().basis()
-        # Need to specify the dimensions explicitly so that we don't
-        # wind up with a zero-by-zero matrix when we want e.g. a
-        # zero-by-two matrix (important for composing things).
-        P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
-        P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
-        pi_left = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J1,P1)
-        pi_right = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J2,P2)
-        return (pi_left, pi_right)
-
-    def inclusions(self):
-        r"""
-        Return the pair of inclusion maps from our factors into us.
+        offset = sum( self.cartesian_factor(k).dimension()
+                      for k in range(i) )
+        Ji = self.cartesian_factor(i)
+        Pi = self._module_morphism(lambda j: Ji.monomial(j - offset),
+                                   codomain=Ji)

 
 

+        return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix())
+
+    @cached_method
+    def cartesian_embedding(self, i):
+        r"""

         SETUP::
 
             sage: from mjo.eja.eja_algebra import (random_eja,
             ....:                                  JordanSpinEJA,
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (random_eja,
             ....:                                  JordanSpinEJA,

-            ....:                                  RealSymmetricEJA,
-            ....:                                  DirectSumEJA)
+            ....:                                  HadamardEJA,
+            ....:                                  RealSymmetricEJA)

 
 

-        EXAMPLES::
+        EXAMPLES:
+
+        The embedding morphisms are Euclidean Jordan algebra
+        operators::
+
+            sage: J1 = HadamardEJA(2)
+            sage: J2 = RealSymmetricEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: J.cartesian_embedding(0)
+            Linear operator between finite-dimensional Euclidean Jordan
+            algebras represented by the matrix:
+            [1 0]
+            [0 1]
+            [0 0]
+            [0 0]
+            [0 0]
+            Domain: Euclidean Jordan algebra of dimension 2 over
+            Algebraic Real Field
+            Codomain: Euclidean Jordan algebra of dimension 2 over
+            Algebraic Real Field (+) Euclidean Jordan algebra of
+            dimension 3 over Algebraic Real Field
+            sage: J.cartesian_embedding(1)
+            Linear operator between finite-dimensional Euclidean Jordan
+            algebras represented by the matrix:
+            [0 0 0]
+            [0 0 0]
+            [1 0 0]
+            [0 1 0]
+            [0 0 1]
+            Domain: Euclidean Jordan algebra of dimension 3 over
+            Algebraic Real Field
+            Codomain: Euclidean Jordan algebra of dimension 2 over
+            Algebraic Real Field (+) Euclidean Jordan algebra of
+            dimension 3 over Algebraic Real Field
+
+        The embeddings work the way you'd expect on the vector
+        representation of an element::

 
             sage: J1 = JordanSpinEJA(3)
             sage: J2 = RealSymmetricEJA(2)
 
             sage: J1 = JordanSpinEJA(3)
             sage: J2 = RealSymmetricEJA(2)

-            sage: J = DirectSumEJA(J1,J2)
-            sage: (iota_left, iota_right) = J.inclusions()
+            sage: J = cartesian_product([J1,J2])
+            sage: iota_left = J.cartesian_embedding(0)
+            sage: iota_right = J.cartesian_embedding(1)

             sage: iota_left(J1.zero()) == J.zero()
             True
             sage: iota_right(J2.zero()) == J.zero()
             sage: iota_left(J1.zero()) == J.zero()
             True
             sage: iota_right(J2.zero()) == J.zero()


@@ -2753,6 +3136,17 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
 
         TESTS:
 
 
         TESTS:
 

+        The answer never changes::
+
+            sage: set_random_seed()
+            sage: J1 = random_eja()
+            sage: J2 = random_eja()
+            sage: J = cartesian_product([J1,J2])
+            sage: E0 = J.cartesian_embedding(0)
+            sage: E1 = J.cartesian_embedding(0)
+            sage: E0 == E1
+            True
+

         Composing a projection with the corresponding inclusion should
         produce the identity map, and mismatching them should produce
         the zero map::
         Composing a projection with the corresponding inclusion should
         produce the identity map, and mismatching them should produce
         the zero map::


@@ -2760,9 +3154,11 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
             sage: set_random_seed()
             sage: J1 = random_eja()
             sage: J2 = random_eja()
             sage: set_random_seed()
             sage: J1 = random_eja()
             sage: J2 = random_eja()

-            sage: J = DirectSumEJA(J1,J2)
-            sage: (iota_left, iota_right) = J.inclusions()
-            sage: (pi_left, pi_right) = J.projections()
+            sage: J = cartesian_product([J1,J2])
+            sage: iota_left = J.cartesian_embedding(0)
+            sage: iota_right = J.cartesian_embedding(1)
+            sage: pi_left = J.cartesian_projection(0)
+            sage: pi_right = J.cartesian_projection(1)

             sage: pi_left*iota_left == J1.one().operator()
             True
             sage: pi_right*iota_right == J2.one().operator()
             sage: pi_left*iota_left == J1.one().operator()
             True
             sage: pi_right*iota_right == J2.one().operator()


@@ -2773,58 +3169,83 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
             True
 
         """
             True
 
         """

-        (J1,J2) = self.factors()
-        m = J1.dimension()
-        n = J2.dimension()
-        V_basis = self.vector_space().basis()
-        # Need to specify the dimensions explicitly so that we don't
-        # wind up with a zero-by-zero matrix when we want e.g. a
-        # two-by-zero matrix (important for composing things).
-        I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m])
-        I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:])
-        iota_left = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,self,I1)
-        iota_right = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,self,I2)
-        return (iota_left, iota_right)
+        offset = sum( self.cartesian_factor(k).dimension()
+                      for k in range(i) )
+        Ji = self.cartesian_factor(i)
+        Ei = Ji._module_morphism(lambda j: self.monomial(j + offset),
+                                 codomain=self)
+        return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix())

 
 

-    def inner_product(self, x, y):
-        r"""
-        The standard Cartesian inner-product.

 
 

-        We project ``x`` and ``y`` onto our factors, and add up the
-        inner-products from the subalgebras.

 
 

-        SETUP::
+FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA

 
 

+class RationalBasisCartesianProductEJA(CartesianProductEJA,
+                                       RationalBasisEJA):
+    r"""
+    A separate class for products of algebras for which we know a
+    rational basis.

 
 

-            sage: from mjo.eja.eja_algebra import (HadamardEJA,
-            ....:                                  QuaternionHermitianEJA,
-            ....:                                  DirectSumEJA)
-
-        EXAMPLE::
-
-            sage: J1 = HadamardEJA(3,QQ)
-            sage: J2 = QuaternionHermitianEJA(2,QQ,orthonormalize=False)
-            sage: J = DirectSumEJA(J1,J2)
-            sage: x1 = J1.one()
-            sage: x2 = x1
-            sage: y1 = J2.one()
-            sage: y2 = y1
-            sage: x1.inner_product(x2)
-            3
-            sage: y1.inner_product(y2)
-            2
-            sage: J.one().inner_product(J.one())
-            5
+    SETUP::

 
 

-        """
-        (pi_left, pi_right) = self.projections()
-        x1 = pi_left(x)
-        x2 = pi_right(x)
-        y1 = pi_left(y)
-        y2 = pi_right(y)
+        sage: from mjo.eja.eja_algebra import (HadamardEJA,
+        ....:                                  JordanSpinEJA,
+        ....:                                  OctonionHermitianEJA,
+        ....:                                  RealSymmetricEJA)
+
+    EXAMPLES:
+
+    This gives us fast characteristic polynomial computations in
+    product algebras, too::
+
+
+        sage: J1 = JordanSpinEJA(2)
+        sage: J2 = RealSymmetricEJA(3)
+        sage: J = cartesian_product([J1,J2])
+        sage: J.characteristic_polynomial_of().degree()
+        5
+        sage: J.rank()
+        5

 
 

-        return (x1.inner_product(y1) + x2.inner_product(y2))
+    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)

 
 

-random_eja = ConcreteEuclideanJordanAlgebra.random_instance
+        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
+                ])
+
+
+RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
+
+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