eja: factor a few things out of the FDEJA constructor.

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

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

At a minimum, the following are required to construct a Euclidean
Jordan algebra:

  * A basis of matrices, column vectors, or MatrixAlgebra elements
  * A Jordan product defined on the basis
  * Its inner product defined on the basis

The real numbers form a Euclidean Jordan algebra when both the Jordan
and inner products are the usual multiplication. We use this as our
example, and demonstrate a few ways to construct an EJA.

First, we can use one-by-one SageMath matrices with algebraic real
entries to represent real numbers. We define the Jordan and inner
products to be essentially real-number multiplication, with the only
difference being that the Jordan product again returns a one-by-one
matrix, whereas the inner product must return a scalar. Our basis for
the one-by-one matrices is of course the set consisting of a single
matrix with its sole entry non-zero::

    sage: from mjo.eja.eja_algebra import FiniteDimensionalEJA
    sage: jp = lambda X,Y: X*Y
    sage: ip = lambda X,Y: X[0,0]*Y[0,0]
    sage: b1 = matrix(AA, [[1]])
    sage: J1 = FiniteDimensionalEJA((b1,), jp, ip)
    sage: J1
    Euclidean Jordan algebra of dimension 1 over Algebraic Real Field

In fact, any positive scalar multiple of that inner-product would work::

    sage: ip2 = lambda X,Y: 16*ip(X,Y)
    sage: J2 = FiniteDimensionalEJA((b1,), jp, ip2)
    sage: J2
    Euclidean Jordan algebra of dimension 1 over Algebraic Real Field

But beware that your basis will be orthonormalized _with respect to the
given inner-product_ unless you pass ``orthonormalize=False`` to the
constructor. For example::

    sage: J3 = FiniteDimensionalEJA((b1,), jp, ip2, orthonormalize=False)
    sage: J3
    Euclidean Jordan algebra of dimension 1 over Algebraic Real Field

To see the difference, you can take the first and only basis element
of the resulting algebra, and ask for it to be converted back into
matrix form::

    sage: J1.basis()[0].to_matrix()
    [1]
    sage: J2.basis()[0].to_matrix()
    [1/4]
    sage: J3.basis()[0].to_matrix()
    [1]

Since square roots are used in that process, the default scalar field
that we use is the field of algebraic real numbers, ``AA``. You can
also Use rational numbers, but only if you either pass
``orthonormalize=False`` or know that orthonormalizing your basis
won't stray beyond the rational numbers. The example above would
have worked only because ``sqrt(16) == 4`` is rational.

Another option for your basis is to use elemebts of a
:class:`MatrixAlgebra`::

    sage: from mjo.matrix_algebra import MatrixAlgebra
    sage: A = MatrixAlgebra(1,AA,AA)
    sage: J4 = FiniteDimensionalEJA(A.gens(), jp, ip)
    sage: J4
    Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
    sage: J4.basis()[0].to_matrix()
    +---+
    | 1 |
    +---+

An easier way to view the entire EJA basis in its original (but
perhaps orthonormalized) matrix form is to use the ``matrix_basis``
method::

    sage: J4.matrix_basis()
    (+---+
    | 1 |
     +---+,)

In particular, a :class:`MatrixAlgebra` is needed to work around the
fact that matrices in SageMath must have entries in the same
(commutative and associative) ring as its scalars. There are many
Euclidean Jordan algebras whose elements are matrices that violate
those assumptions. The complex, quaternion, and octonion Hermitian
matrices all have entries in a ring (the complex numbers, quaternions,
or octonions...) that differs from the algebra's scalar ring (the real
numbers). Quaternions are also non-commutative; the octonions are
neither commutative nor associative.

SETUP::

    sage: from mjo.eja.eja_algebra import random_eja

EXAMPLES::

    sage: random_eja()
    Euclidean Jordan algebra of dimension...
"""

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.misc.cachefunc import cached_method
from sage.misc.table import table
from sage.modules.free_module import FreeModule, VectorSpace
from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
                            PolynomialRing,
                            QuadraticField)
from mjo.eja.eja_element import FiniteDimensionalEJAElement
from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
from mjo.eja.eja_utils import _all2list, _mat2vec

def EuclideanJordanAlgebras(field):
    r"""
    The category of Euclidean Jordan algebras over ``field``, which
    must be a subfield of the real numbers. For now this is just a
    convenient wrapper around all of the other category axioms that
    apply to all EJAs.
    """
    category = MagmaticAlgebras(field).FiniteDimensional()
    category = category.WithBasis().Unital().Commutative()
    return category

class FiniteDimensionalEJA(CombinatorialFreeModule):
    r"""
    A finite-dimensional Euclidean Jordan algebra.

    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.

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

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

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

    SETUP::

        sage: from mjo.eja.eja_algebra import random_eja

    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(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
    """
    Element = FiniteDimensionalEJAElement

    @staticmethod
    def _check_input_field(field):
        if not field.is_subring(RR):
            # Note: this does return true for the real algebraic
            # field, the rationals, and any quadratic field where
            # we've specified a real embedding.
            raise ValueError("scalar field is not real")

    @staticmethod
    def _check_input_axioms(basis, jordan_product, inner_product):
        if not all( jordan_product(bi,bj) == jordan_product(bj,bi)
                    for bi in basis
                    for bj in basis ):
            raise ValueError("Jordan product is not commutative")

        if not all( inner_product(bi,bj) == inner_product(bj,bi)
                    for bi in basis
                    for bj in basis ):
            raise ValueError("inner-product is not commutative")

    def __init__(self,
                 basis,
                 jordan_product,
                 inner_product,
                 field=AA,
                 matrix_space=None,
                 orthonormalize=True,
                 associative=None,
                 cartesian_product=False,
                 check_field=True,
                 check_axioms=True,
                 prefix="b"):

        n = len(basis)

        if check_field:
            self._check_input_field(field)

        if check_axioms:
            # Check commutativity of the Jordan and inner-products.
            # This has to be done before we build the multiplication
            # and inner-product tables/matrices, because we take
            # advantage of symmetry in the process.
            self._check_input_axioms(basis, jordan_product, inner_product)

        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)

        category = EuclideanJordanAlgebras(field)

        if associative:
            # Element subalgebras can take advantage of this.
            category = category.Associative()
        if cartesian_product:
            # Use join() here because otherwise we only get the
            # "Cartesian product of..." and not the things themselves.
            category = category.join([category,
                                      category.CartesianProducts()])

        # Call the superclass constructor so that we can use its from_vector()
        # method to build our multiplication table.
        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.

        degree = 0
        if n > 0:
            degree = len(_all2list(basis[0]))

        # Build an ambient space that fits our matrix basis when
        # written out as "long vectors."
        V = VectorSpace(field, degree)

        # The matrix that will hold the orthonormal -> unorthonormal
        # coordinate transformation. Default to an identity matrix of
        # the appropriate size to avoid special cases for None
        # everywhere.
        self._deortho_matrix = matrix.identity(field,n)

        if orthonormalize:
            # Save a copy of the un-orthonormalized basis for later.
            # 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 )

            from mjo.eja.eja_utils import gram_schmidt
            basis = tuple(gram_schmidt(basis, inner_product))

        # Save the (possibly orthonormalized) matrix basis for
        # later, as well as the space that its elements live in.
        # In most cases we can deduce the matrix space, but when
        # n == 0 (that is, there are no basis elements) we cannot.
        self._matrix_basis = basis
        if matrix_space is None:
            self._matrix_space = self._matrix_basis[0].parent()
        else:
            self._matrix_space = matrix_space

        # Now create the vector space for the algebra, which will have
        # its own set of non-ambient coordinates (in terms of the
        # supplied basis).
        vector_basis = tuple( V(_all2list(b)) for b in basis )

        # Save the span of our matrix basis (when written out as long
        # vectors) because otherwise we'll have to reconstruct it
        # every time we want to coerce a matrix into the algebra.
        self._matrix_span = V.span_of_basis( vector_basis, check=check_axioms)

        if orthonormalize:
            # Now "self._matrix_span" is the vector space of our
            # algebra coordinates. The variables "X1", "X2",...  refer
            # to the entries of vectors in self._matrix_span. Thus to
            # convert back and forth between the orthonormal
            # coordinates and the given ones, we need to stick the
            # original basis in self._matrix_span.
            U = V.span_of_basis( deortho_vector_basis, check=check_axioms)
            self._deortho_matrix = matrix.column( U.coordinate_vector(q)
                                                  for q in vector_basis )


        # Now we actually compute the multiplication and inner-product
        # tables/matrices using the possibly-orthonormalized basis.
        self._inner_product_matrix = matrix.identity(field, n)
        self._multiplication_table = [ [0 for j in range(i+1)]
                                       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 = basis[i]
                q_j = basis[j]

                # 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 = self._matrix_span.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

        self._inner_product_matrix._cache = {'hermitian': True}
        self._inner_product_matrix.set_immutable()

        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):
            ...
            ValueError: not an element of this algebra

        """
        return None


    def product_on_basis(self, i, j):
        r"""
        Returns the Jordan product of the `i` and `j`th basis elements.

        This completely defines the Jordan product on the algebra, and
        is used direclty by our superclass machinery to implement
        :meth:`product`.

        SETUP::

            sage: from mjo.eja.eja_algebra import random_eja

        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

        """
        # 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]

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

            sage: from mjo.eja.eja_algebra import (random_eja,
            ....:                                  HadamardEJA,
            ....:                                  BilinearFormEJA)

        EXAMPLES:

        Our inner product is "associative," which means the following for
        a symmetric bilinear form::

            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

        TESTS:

        Ensure that this is the usual inner product for the algebras
        over `R^n`::

            sage: set_random_seed()
            sage: J = HadamardEJA.random_instance()
            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

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

            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

        """
        B = self._inner_product_matrix
        return (B*x.to_vector()).inner_product(y.to_vector())


    def is_associative(self):
        r"""
        Return whether or not this algebra's Jordan product is associative.

        SETUP::

            sage: from mjo.eja.eja_algebra import ComplexHermitianEJA

        EXAMPLES::

            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

        """
        return "Associative" in self.category().axioms()

    def _is_commutative(self):
        r"""
        Whether or not this algebra's multiplication table is commutative.

        This method should of course always return ``True``, unless
        this algebra was constructed with ``check_axioms=False`` and
        passed an invalid multiplication table.
        """
        return all( x*y == y*x for x in self.gens() for y in self.gens() )

    def _is_jordanian(self):
        r"""
        Whether or not this algebra's multiplication table respects the
        Jordan identity `(x^{2})(xy) = x(x^{2}y)`.

        We only check one arrangement of `x` and `y`, so for a
        ``True`` result to be truly true, you should also check
        :meth:`_is_commutative`. This method should of course always
        return ``True``, unless this algebra was constructed with
        ``check_axioms=False`` and passed an invalid multiplication table.
        """
        return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
                    ==
                    (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j))
                    for i in range(self.dimension())
                    for j in range(self.dimension()) )

    def _jordan_product_is_associative(self):
        r"""
        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 agree with :meth:`is_associative` unless
        you lied about the value of the ``associative`` parameter
        when you constructed the algebra.

        SETUP::

            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)
                    y = self.monomial(j)
                    z = self.monomial(k)
                    diff = (x*y).inner_product(z) - x.inner_product(y*z)

                    if diff.abs() > epsilon:
                        return False

        return True

    def _element_constructor_(self, elt):
        """
        Construct an element of this algebra or a subalgebra from its
        EJA element, vector, or matrix representation.

        This gets called only after the parent element _call_ method
        fails to find a coercion for the argument.

        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

        Subalgebra elements are embedded into the superalgebra::

            sage: J = JordanSpinEJA(3)
            sage: J.one()
            b0
            sage: x = sum(J.gens())
            sage: A = x.subalgebra_generated_by()
            sage: J(A.one())
            b0

        TESTS:

        Ensure that we can convert any element back and forth
        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)

        if hasattr(elt, 'superalgebra_element'):
            # Handle subalgebra elements
            if elt.parent().superalgebra() == self:
                return elt.superalgebra_element()

        if hasattr(elt, 'sparse_vector'):
            # Convert a vector into a column-matrix. We check for
            # "sparse_vector" and not "column" because matrices also
            # have a "column" method.
            elt = elt.column()

        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.
        #
        elt = self._matrix_span.ambient_vector_space()(_all2list(elt))

        try:
            coords = self._matrix_span.coordinate_vector(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):
        """
        Return the algebra's "characteristic polynomial of" function,
        which is itself a multivariate polynomial that, when evaluated
        at the coordinates of some algebra element, returns that
        element's characteristic polynomial.

        The resulting polynomial has `n+1` variables, where `n` is the
        dimension of this algebra. The first `n` variables correspond to
        the coordinates of an algebra element: when evaluated at the
        coordinates of an algebra element with respect to a certain
        basis, the result is a univariate polynomial (in the one
        remaining variable ``t``), namely the characteristic polynomial
        of that element.

        SETUP::

            sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA

        EXAMPLES:

        The characteristic polynomial in the spin algebra is given in
        Alizadeh, Example 11.11::

            sage: J = JordanSpinEJA(3)
            sage: p = J.characteristic_polynomial_of(); p
            X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
            sage: xvec = J.one().to_vector()
            sage: p(*xvec)
            t^2 - 2*t + 1

        By definition, the characteristic polynomial is a monic
        degree-zero polynomial in a rank-zero algebra. Note that
        Cayley-Hamilton is indeed satisfied since the polynomial
        ``1`` evaluates to the identity element of the algebra on
        any argument::

            sage: J = TrivialEJA()
            sage: J.characteristic_polynomial_of()
            1

        """
        r = self.rank()
        n = self.dimension()

        # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
        a = self._charpoly_coefficients()

        # We go to a bit of trouble here to reorder the
        # indeterminates, so that it's easier to evaluate the
        # characteristic polynomial at x's coordinates and get back
        # something in terms of t, which is what we want.
        S = PolynomialRing(self.base_ring(),'t')
        t = S.gen(0)
        if r > 0:
            R = a[0].parent()
            S = PolynomialRing(S, R.variable_names())
            t = S(t)

        return (t**r + sum( a[k]*(t**k) for k in range(r) ))

    def coordinate_polynomial_ring(self):
        r"""
        The multivariate polynomial ring in which this algebra's
        :meth:`characteristic_polynomial_of` lives.

        SETUP::

            sage: from mjo.eja.eja_algebra import (HadamardEJA,
            ....:                                  RealSymmetricEJA)

        EXAMPLES::

            sage: J = HadamardEJA(2)
            sage: J.coordinate_polynomial_ring()
            Multivariate Polynomial Ring in X1, X2...
            sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
            sage: J.coordinate_polynomial_ring()
            Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...

        """
        var_names = tuple( "X%d" % z for z in range(1, self.dimension()+1) )
        return PolynomialRing(self.base_ring(), var_names)

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

            sage: from mjo.eja.eja_algebra import (random_eja,
            ....:                                  HadamardEJA,
            ....:                                  BilinearFormEJA)

        EXAMPLES:

        Our inner product is "associative," which means the following for
        a symmetric bilinear form::

            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

        TESTS:

        Ensure that this is the usual inner product for the algebras
        over `R^n`::

            sage: set_random_seed()
            sage: J = HadamardEJA.random_instance()
            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

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

            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
        """
        B = self._inner_product_matrix
        return (B*x.to_vector()).inner_product(y.to_vector())


    def is_trivial(self):
        """
        Return whether or not this algebra is trivial.

        A trivial algebra contains only the zero element.

        SETUP::

            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
            ....:                                  TrivialEJA)

        EXAMPLES::

            sage: J = ComplexHermitianEJA(3)
            sage: J.is_trivial()
            False

        ::

            sage: J = TrivialEJA()
            sage: J.is_trivial()
            True

        """
        return self.dimension() == 0


    def multiplication_table(self):
        """
        Return a visual representation of this algebra's multiplication
        table (on basis elements).

        SETUP::

            sage: from mjo.eja.eja_algebra import JordanSpinEJA

        EXAMPLES::

            sage: J = JordanSpinEJA(4)
            sage: J.multiplication_table()
            +----++----+----+----+----+
            | *  || b0 | b1 | b2 | b3 |
            +====++====+====+====+====+
            | b0 || b0 | b1 | b2 | b3 |
            +----++----+----+----+----+
            | b1 || b1 | b0 | 0  | 0  |
            +----++----+----+----+----+
            | b2 || b2 | 0  | b0 | 0  |
            +----++----+----+----+----+
            | b3 || b3 | 0  | 0  | b0 |
            +----++----+----+----+----+

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

        # 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) ]

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


    def matrix_basis(self):
        """
        Return an (often more natural) representation of this algebras
        basis as an ordered tuple of matrices.

        Every finite-dimensional Euclidean Jordan Algebra is a, up to
        Jordan isomorphism, a direct sum of five simple
        algebras---four of which comprise Hermitian matrices. And the
        last type of algebra can of course be thought of as `n`-by-`1`
        column matrices (ambiguusly called column vectors) to avoid
        special cases. As a result, matrices (and column vectors) are
        a natural representation format for Euclidean Jordan algebra
        elements.

        But, when we construct an algebra from a basis of matrices,
        those matrix representations are lost in favor of coordinate
        vectors *with respect to* that basis. We could eventually
        convert back if we tried hard enough, but having the original
        representations handy is valuable enough that we simply store
        them and return them from this method.

        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:`CartesianProductEJA` can be displayed
        nicely, without having to have special classes for direct sums
        one of whose components was a matrix algebra.

        SETUP::

            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
            ....:                                  RealSymmetricEJA)

        EXAMPLES::

            sage: J = RealSymmetricEJA(2)
            sage: J.basis()
            Finite family {0: b0, 1: b1, 2: b2}
            sage: J.matrix_basis()
            (
            [1 0]  [                  0 0.7071067811865475?]  [0 0]
            [0 0], [0.7071067811865475?                   0], [0 1]
            )

        ::

            sage: J = JordanSpinEJA(2)
            sage: J.basis()
            Finite family {0: b0, 1: b1}
            sage: J.matrix_basis()
            (
            [1]  [0]
            [0], [1]
            )
        """
        return self._matrix_basis


    def matrix_space(self):
        """
        Return the matrix space in which this algebra's elements live, if
        we think of them as matrices (including column vectors of the
        appropriate size).

        "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
        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()
            Module of 2 by 2 matrices with entries in Algebraic Field over
            the scalar ring Rational Field
            sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
            sage: J.matrix_space()
            Module of 1 by 1 matrices with entries in Quaternion
            Algebra (-1, -1) with base ring Rational Field over
            the scalar ring Rational Field

        """
        return self._matrix_space


    @cached_method
    def one(self):
        """
        Return the unit element of this algebra.

        SETUP::

            sage: from mjo.eja.eja_algebra import (HadamardEJA,
            ....:                                  random_eja)

        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()
            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:

        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: 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,
        regardless of the base field and whether or not we
        orthonormalize::

            sage: set_random_seed()
            sage: J = random_eja()
            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()
            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::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: cached = J.one()
            sage: J.one.clear_cache()
            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.
        # If some linear combination of those basis elements is the
        # algebra identity, then the same linear combination of
        # their matrices has to be the identity matrix.
        #
        # Of course, matrices aren't vectors in sage, so we have to
        # appeal to the "long vectors" isometry.
        oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]

        # Now we use basic linear algebra to find the coefficients,
        # of the matrices-as-vectors-linear-combination, which should
        # work for the original algebra basis too.
        A = matrix(self.base_ring(), oper_vecs)

        # We used the isometry on the left-hand side already, but we
        # still need to do it for the right-hand side. Recall that we
        # wanted something that summed to the identity matrix.
        b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) )

        # Now if there's an identity element in the algebra, this
        # should work. We solve on the left to avoid having to
        # transpose the matrix "A".
        return self.from_vector(A.solve_left(b))


    def peirce_decomposition(self, c):
        """
        The Peirce decomposition of this algebra relative to the
        idempotent ``c``.

        In the future, this can be extended to a complete system of
        orthogonal idempotents.

        INPUT:

          - ``c`` -- an idempotent of this algebra.

        OUTPUT:

        A triple (J0, J5, J1) containing two subalgebras and one subspace
        of this algebra,

          - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
            corresponding to the eigenvalue zero.

          - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
            corresponding to the eigenvalue one-half.

          - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
            corresponding to the eigenvalue one.

        These are the only possible eigenspaces for that operator, and this
        algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
        orthogonal, and are subalgebras of this algebra with the appropriate
        restrictions.

        SETUP::

            sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA

        EXAMPLES:

        The canonical example comes from the symmetric matrices, which
        decompose into diagonal and off-diagonal parts::

            sage: J = RealSymmetricEJA(3)
            sage: C = matrix(QQ, [ [1,0,0],
            ....:                  [0,1,0],
            ....:                  [0,0,0] ])
            sage: c = J(C)
            sage: J0,J5,J1 = J.peirce_decomposition(c)
            sage: J0
            Euclidean Jordan algebra of dimension 1...
            sage: J5
            Vector space of degree 6 and dimension 2...
            sage: J1
            Euclidean Jordan algebra of dimension 3...
            sage: J0.one().to_matrix()
            [0 0 0]
            [0 0 0]
            [0 0 1]
            sage: orig_df = AA.options.display_format
            sage: AA.options.display_format = 'radical'
            sage: J.from_vector(J5.basis()[0]).to_matrix()
            [          0           0 1/2*sqrt(2)]
            [          0           0           0]
            [1/2*sqrt(2)           0           0]
            sage: J.from_vector(J5.basis()[1]).to_matrix()
            [          0           0           0]
            [          0           0 1/2*sqrt(2)]
            [          0 1/2*sqrt(2)           0]
            sage: AA.options.display_format = orig_df
            sage: J1.one().to_matrix()
            [1 0 0]
            [0 1 0]
            [0 0 0]

        TESTS:

        Every algebra decomposes trivially with respect to its identity
        element::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: J0,J5,J1 = J.peirce_decomposition(J.one())
            sage: J0.dimension() == 0 and J5.dimension() == 0
            True
            sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
            True

        The decomposition is into eigenspaces, and its components are
        therefore necessarily orthogonal. Moreover, the identity
        elements in the two subalgebras are the projections onto their
        respective subspaces of the superalgebra's identity element::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: x = J.random_element()
            sage: if not J.is_trivial():
            ....:     while x.is_nilpotent():
            ....:         x = J.random_element()
            sage: c = x.subalgebra_idempotent()
            sage: J0,J5,J1 = J.peirce_decomposition(c)
            sage: ipsum = 0
            sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
            ....:     w = w.superalgebra_element()
            ....:     y = J.from_vector(y)
            ....:     z = z.superalgebra_element()
            ....:     ipsum += w.inner_product(y).abs()
            ....:     ipsum += w.inner_product(z).abs()
            ....:     ipsum += y.inner_product(z).abs()
            sage: ipsum
            0
            sage: J1(c) == J1.one()
            True
            sage: J0(J.one() - c) == J0.one()
            True

        """
        if not c.is_idempotent():
            raise ValueError("element is not idempotent: %s" % c)

        # 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 = self.subalgebra((), check_axioms=False)
        J0 = trivial                          # eigenvalue zero
        J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
        J1 = trivial                          # eigenvalue one

        for (eigval, eigspace) in c.operator().matrix().right_eigenspaces():
            if eigval == ~(self.base_ring()(2)):
                J5 = eigspace
            else:
                gens = tuple( self.from_vector(b) for b in eigspace.basis() )
                subalg = self.subalgebra(gens, check_axioms=False)
                if eigval == 0:
                    J0 = subalg
                elif eigval == 1:
                    J1 = subalg
                else:
                    raise ValueError("unexpected eigenvalue: %s" % eigval)

        return (J0, J5, J1)


    def random_element(self, thorough=False):
        r"""
        Return a random element of this algebra.

        Our algebra superclass method only returns a linear
        combination of at most two basis elements. We instead
        want the vector space "random element" method that
        returns a more diverse selection.

        INPUT:

        - ``thorough`` -- (boolean; default False) whether or not we
          should generate irrational coefficients for the random
          element when our base ring is irrational; this slows the
          algebra operations to a crawl, but any truly random method
          should include them

        """
        # For a general base ring... maybe we can trust this to do the
        # right thing? Unlikely, but.
        V = self.vector_space()
        v = V.random_element()

        if self.base_ring() is AA:
            # The "random element" method of the algebraic reals is
            # stupid at the moment, and only returns integers between
            # -2 and 2, inclusive:
            #
            #   https://trac.sagemath.org/ticket/30875
            #
            # Instead, we implement our own "random vector" method,
            # and then coerce that into the algebra. We use the vector
            # space degree here instead of the dimension because a
            # subalgebra could (for example) be spanned by only two
            # vectors, each with five coordinates.  We need to
            # generate all five coordinates.
            if thorough:
                v *= QQbar.random_element().real()
            else:
                v *= QQ.random_element()

        return self.from_vector(V.coordinate_vector(v))

    def random_elements(self, count, thorough=False):
        """
        Return ``count`` random elements as a tuple.

        INPUT:

        - ``thorough`` -- (boolean; default False) whether or not we
          should generate irrational coefficients for the random
          elements when our base ring is irrational; this slows the
          algebra operations to a crawl, but any truly random method
          should include them

        SETUP::

            sage: from mjo.eja.eja_algebra import JordanSpinEJA

        EXAMPLES::

            sage: J = JordanSpinEJA(3)
            sage: x,y,z = J.random_elements(3)
            sage: all( [ x in J, y in J, z in J ])
            True
            sage: len( J.random_elements(10) ) == 10
            True

        """
        return tuple( self.random_element(thorough)
                      for idx in range(count) )


    @cached_method
    def _charpoly_coefficients(self):
        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()
        vars = R.gens()
        F = R.fraction_field()

        def L_x_i_j(i,j):
            # From a result in my book, these are the entries of the
            # basis representation of L_x.
            return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
                        for k in range(n) )

        L_x = matrix(F, n, n, L_x_i_j)

        r = None
        if self.rank.is_in_cache():
            r = self.rank()
            # There's no need to pad the system with redundant
            # columns if we *know* they'll be redundant.
            n = r

        # Compute an extra power in case the rank is equal to
        # the dimension (otherwise, we would stop at x^(r-1)).
        x_powers = [ (L_x**k)*self.one().to_vector()
                     for k in range(n+1) ]
        A = matrix.column(F, x_powers[:n])
        AE = A.extended_echelon_form()
        E = AE[:,n:]
        A_rref = AE[:,:n]
        if r is None:
            r = A_rref.rank()
        b = x_powers[r]

        # 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. 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):
        r"""
        Return the rank of this EJA.

        This is a cached method because we know the rank a priori for
        all of the algebras we can construct. Thus we can avoid the
        expensive ``_charpoly_coefficients()`` call unless we truly
        need to compute the whole characteristic polynomial.

        SETUP::

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

        EXAMPLES:

        The rank of the Jordan spin algebra is always two::

            sage: JordanSpinEJA(2).rank()
            2
            sage: JordanSpinEJA(3).rank()
            2
            sage: JordanSpinEJA(4).rank()
            2

        The rank of the `n`-by-`n` Hermitian real, complex, or
        quaternion matrices is `n`::

            sage: RealSymmetricEJA(4).rank()
            4
            sage: ComplexHermitianEJA(3).rank()
            3
            sage: QuaternionHermitianEJA(2).rank()
            2

        TESTS:

        Ensure that every EJA that we know how to construct has a
        positive integer rank, unless the algebra is trivial in
        which case its rank will be zero::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: r = J.rank()
            sage: r in ZZ
            True
            sage: r > 0 or (r == 0 and J.is_trivial())
            True

        Ensure that computing the rank actually works, since the ranks
        of all simple algebras are known and will be cached by default::

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

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

        SETUP::

            sage: from mjo.eja.eja_algebra import RealSymmetricEJA

        EXAMPLES::

            sage: J = RealSymmetricEJA(2)
            sage: J.vector_space()
            Vector space of dimension 3 over...

        """
        return self.zero().to_vector().parent().ambient_vector_space()



class RationalBasisEJA(FiniteDimensionalEJA):
    r"""
    Algebras whose supplied basis elements have all rational entries.

    SETUP::

        sage: from mjo.eja.eja_algebra import BilinearFormEJA

    EXAMPLES:

    The supplied basis is orthonormalized by default::

        sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
        sage: J = BilinearFormEJA(B)
        sage: J.matrix_basis()
        (
        [1]  [  0]  [   0]
        [0]  [1/5]  [32/5]
        [0], [  0], [   5]
        )

    """
    def __init__(self,
                 basis,
                 jordan_product,
                 inner_product,
                 field=AA,
                 check_field=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.
            # Use _all2list to get the vector coordinates of octonion
            # entries and not the octonions themselves (which are not
            # rational).
            if not all( all(b_i in QQ for b_i in _all2list(b))
                        for b in basis ):
                raise TypeError("basis not rational")

        super().__init__(basis,
                         jordan_product,
                         inner_product,
                         field=field,
                         check_field=check_field,
                         **kwargs)

        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,
                                       matrix_space=self.matrix_space(),
                                       associative=self.is_associative(),
                                       orthonormalize=False,
                                       check_field=False,
                                       check_axioms=False)

    @cached_method
    def _charpoly_coefficients(self):
        r"""
        SETUP::

            sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
            ....:                                  JordanSpinEJA)

        EXAMPLES:

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

            sage: J = JordanSpinEJA(3)
            sage: J._charpoly_coefficients()
            (X1^2 - X2^2 - X3^2, -2*X1)
            sage: a0 = J._charpoly_coefficients()[0]
            sage: J.base_ring()
            Algebraic Real Field
            sage: a0.base_ring()
            Algebraic Real Field

        """
        if self._rational_algebra is None:
            # There's no need to construct *another* algebra over the
            # rationals if this one is already over the
            # rationals. Likewise, if we never orthonormalized our
            # basis, we might as well just use the given one.
            return super()._charpoly_coefficients()

        # Do the computation over the rationals. The answer will be
        # the same, because all we've done is a change of basis.
        # Then, change back from QQ to our real base ring
        a = ( a_i.change_ring(self.base_ring())
              for a_i in self._rational_algebra._charpoly_coefficients() )

        # Otherwise, convert the coordinate variables back to the
        # deorthonormalized ones.
        R = self.coordinate_polynomial_ring()
        from sage.modules.free_module_element import vector
        X = vector(R, R.gens())
        BX = self._deortho_matrix*X

        subs_dict = { X[i]: BX[i] for i in range(len(X)) }
        return tuple( a_i.subs(subs_dict) for a_i in a )

class ConcreteEJA(FiniteDimensionalEJA):
    r"""
    A class for the Euclidean Jordan algebras that we know by name.

    These are the Jordan algebras whose basis, multiplication table,
    rank, and so on are known a priori. More to the point, they are
    the Euclidean Jordan algebras for which we are able to conjure up
    a "random instance."

    SETUP::

        sage: from mjo.eja.eja_algebra import ConcreteEJA

    TESTS:

    Our basis is normalized with respect to the algebra's inner
    product, unless we specify otherwise::

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

    Since our basis is orthonormal with respect to the algebra's inner
    product, and since we know that this algebra is an EJA, any
    left-multiplication operator's matrix will be symmetric because
    natural->EJA basis representation is an isometry and within the
    EJA the operator is self-adjoint by the Jordan axiom::

        sage: set_random_seed()
        sage: J = ConcreteEJA.random_instance()
        sage: x = J.random_element()
        sage: x.operator().is_self_adjoint()
        True
    """

    @staticmethod
    def _max_random_instance_dimension():
        r"""
        The maximum dimension of any random instance. Ten dimensions seems
        to be about the point where everything takes a turn for the
        worse. And dimension ten (but not nine) allows the 4-by-4 real
        Hermitian matrices, the 2-by-2 quaternion Hermitian matrices,
        and the 2-by-2 octonion Hermitian matrices.
        """
        return 10

    @staticmethod
    def _max_random_instance_size(max_dimension):
        """
        Return an integer "size" that is an upper bound on the size of
        this algebra when it is used in a random test case. This size
        (which can be passed to the algebra's constructor) is itself
        based on the ``max_dimension`` parameter.

        This method must be implemented in each subclass.
        """
        raise NotImplementedError

    @classmethod
    def random_instance(cls, max_dimension=None, *args, **kwargs):
        """
        Return a random instance of this type of algebra whose dimension
        is less than or equal to the lesser of ``max_dimension`` and
        the value returned by ``_max_random_instance_dimension()``. If
        the dimension bound is omitted, then only the
        ``_max_random_instance_dimension()`` is used as a bound.

        This method should be implemented in each subclass.

        SETUP::

            sage: from mjo.eja.eja_algebra import ConcreteEJA

        TESTS:

        Both the class bound and the ``max_dimension`` argument are upper
        bounds on the dimension of the algebra returned::

            sage: from sage.misc.prandom import choice
            sage: eja_class = choice(ConcreteEJA.__subclasses__())
            sage: class_max_d = eja_class._max_random_instance_dimension()
            sage: J = eja_class.random_instance(max_dimension=20,
            ....:                               field=QQ,
            ....:                               orthonormalize=False)
            sage: J.dimension() <= class_max_d
            True
            sage: J = eja_class.random_instance(max_dimension=2,
            ....:                               field=QQ,
            ....:                               orthonormalize=False)
            sage: J.dimension() <= 2
            True

        """
        from sage.misc.prandom import choice
        eja_class = choice(cls.__subclasses__())

        # These all bubble up to the RationalBasisEJA superclass
        # constructor, so any (kw)args valid there are also valid
        # here.
        return eja_class.random_instance(max_dimension, *args, **kwargs)


class MatrixEJA(FiniteDimensionalEJA):
    @staticmethod
    def _denormalized_basis(A):
        """
        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.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

        """
        # 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 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.

        SETUP::

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

        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

        """
        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()


    def __init__(self, matrix_space, **kwargs):
        # We know this is a valid EJA, but will double-check
        # if the user passes check_axioms=True.
        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False


        super().__init__(self._denormalized_basis(matrix_space),
                         self.jordan_product,
                         self.trace_inner_product,
                         field=matrix_space.base_ring(),
                         matrix_space=matrix_space,
                         **kwargs)

        self.rank.set_cache(matrix_space.nrows())
        self.one.set_cache( self(matrix_space.one()) )

class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
    """
    The rank-n simple EJA consisting of real symmetric n-by-n
    matrices, the usual symmetric Jordan product, and the trace inner
    product. It has dimension `(n^2 + n)/2` over the reals.

    SETUP::

        sage: from mjo.eja.eja_algebra import RealSymmetricEJA

    EXAMPLES::

        sage: J = RealSymmetricEJA(2)
        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::

        sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
        Euclidean Jordan algebra of dimension 3 over Real Double Field
        sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
        Euclidean Jordan algebra of dimension 3 over Real Field with
        53 bits of precision

    TESTS:

    The dimension of this algebra is `(n^2 + n) / 2`::

        sage: set_random_seed()
        sage: d = RealSymmetricEJA._max_random_instance_dimension()
        sage: n = RealSymmetricEJA._max_random_instance_size(d)
        sage: J = RealSymmetricEJA(n)
        sage: J.dimension() == (n^2 + n)/2
        True

    The Jordan multiplication is what we think it is::

        sage: set_random_seed()
        sage: J = RealSymmetricEJA.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: RealSymmetricEJA(3, prefix='q').gens()
        (q0, q1, q2, q3, q4, q5)

    We can construct the (trivial) algebra of rank zero::

        sage: RealSymmetricEJA(0)
        Euclidean Jordan algebra of dimension 0 over Algebraic Real Field

    """
    @staticmethod
    def _max_random_instance_size(max_dimension):
        # Obtained by solving d = (n^2 + n)/2.
        # The ZZ-int-ZZ thing is just "floor."
        return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/2 - 1/2))

    @classmethod
    def random_instance(cls, max_dimension=None, *args, **kwargs):
        """
        Return a random instance of this type of algebra.
        """
        class_max_d = cls._max_random_instance_dimension()
        if (max_dimension is None or max_dimension > class_max_d):
            max_dimension = class_max_d
        max_size = cls._max_random_instance_size(max_dimension)
        n = ZZ.random_element(max_size + 1)
        return cls(n, **kwargs)

    def __init__(self, n, field=AA, **kwargs):
        # 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__(A, **kwargs)

        from mjo.eja.eja_cache import real_symmetric_eja_coeffs
        a = real_symmetric_eja_coeffs(self)
        if a is not None:
            if self._rational_algebra is None:
                self._charpoly_coefficients.set_cache(a)
            else:
                self._rational_algebra._charpoly_coefficients.set_cache(a)



class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
    """
    The rank-n simple EJA consisting of complex Hermitian n-by-n
    matrices over the real numbers, the usual symmetric Jordan product,
    and the real-part-of-trace inner product. It has dimension `n^2` over
    the reals.

    SETUP::

        sage: from mjo.eja.eja_algebra import ComplexHermitianEJA

    EXAMPLES:

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

        sage: QQbar(RDF(1))
        Traceback (most recent call last):
        ...
        TypeError: Illegal initializer for algebraic number
        sage: AA(RR(1))
        Traceback (most recent call last):
        ...
        TypeError: Illegal initializer for algebraic number

    This causes the following error when we try to scale a matrix of
    complex numbers by an inexact real number::

        sage: ComplexHermitianEJA(2,field=RR)
        Traceback (most recent call last):
        ...
        TypeError: Unable to coerce entries (=(1.00000000000000,
        -0.000000000000000)) to coefficients in Algebraic Real Field

    TESTS:

    The dimension of this algebra is `n^2`::

        sage: set_random_seed()
        sage: d = ComplexHermitianEJA._max_random_instance_dimension()
        sage: n = ComplexHermitianEJA._max_random_instance_size(d)
        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

    """
    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__(A, **kwargs)

        from mjo.eja.eja_cache import complex_hermitian_eja_coeffs
        a = complex_hermitian_eja_coeffs(self)
        if a is not None:
            if self._rational_algebra is None:
                self._charpoly_coefficients.set_cache(a)
            else:
                self._rational_algebra._charpoly_coefficients.set_cache(a)

    @staticmethod
    def _max_random_instance_size(max_dimension):
        # Obtained by solving d = n^2.
        # The ZZ-int-ZZ thing is just "floor."
        return ZZ(int(ZZ(max_dimension).sqrt()))

    @classmethod
    def random_instance(cls, max_dimension=None, *args, **kwargs):
        """
        Return a random instance of this type of algebra.
        """
        class_max_d = cls._max_random_instance_dimension()
        if (max_dimension is None or max_dimension > class_max_d):
            max_dimension = class_max_d
        max_size = cls._max_random_instance_size(max_dimension)
        n = ZZ.random_element(max_size + 1)
        return cls(n, **kwargs)


class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
    r"""
    The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
    matrices, the usual symmetric Jordan product, and the
    real-part-of-trace inner product. It has dimension `2n^2 - n` over
    the reals.

    SETUP::

        sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA

    EXAMPLES:

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

        sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
        Euclidean Jordan algebra of dimension 6 over Real Double Field
        sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
        Euclidean Jordan algebra of dimension 6 over Real Field with
        53 bits of precision

    TESTS:

    The dimension of this algebra is `2*n^2 - n`::

        sage: set_random_seed()
        sage: d = QuaternionHermitianEJA._max_random_instance_dimension()
        sage: n = QuaternionHermitianEJA._max_random_instance_size(d)
        sage: J = QuaternionHermitianEJA(n)
        sage: J.dimension() == 2*(n^2) - n
        True

    The Jordan multiplication is what we think it is::

        sage: set_random_seed()
        sage: J = QuaternionHermitianEJA.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: QuaternionHermitianEJA(2, prefix='a').gens()
        (a0, a1, a2, a3, a4, a5)

    We can construct the (trivial) algebra of rank zero::

        sage: QuaternionHermitianEJA(0)
        Euclidean Jordan algebra of dimension 0 over Algebraic Real Field

    """
    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__(A, **kwargs)

        from mjo.eja.eja_cache import quaternion_hermitian_eja_coeffs
        a = quaternion_hermitian_eja_coeffs(self)
        if a is not None:
            if self._rational_algebra is None:
                self._charpoly_coefficients.set_cache(a)
            else:
                self._rational_algebra._charpoly_coefficients.set_cache(a)



    @staticmethod
    def _max_random_instance_size(max_dimension):
        r"""
        The maximum rank of a random QuaternionHermitianEJA.
        """
        # Obtained by solving d = 2n^2 - n.
        # The ZZ-int-ZZ thing is just "floor."
        return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/4 + 1/4))

    @classmethod
    def random_instance(cls, max_dimension=None, *args, **kwargs):
        """
        Return a random instance of this type of algebra.
        """
        class_max_d = cls._max_random_instance_dimension()
        if (max_dimension is None or max_dimension > class_max_d):
            max_dimension = class_max_d
        max_size = cls._max_random_instance_size(max_dimension)
        n = ZZ.random_element(max_size + 1)
        return cls(n, **kwargs)

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

        sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
        ....:                                  OctonionHermitianEJA)
        sage: from mjo.hurwitz import Octonions, OctonionMatrixAlgebra

    EXAMPLES:

    The 3-by-3 algebra satisfies the axioms of an EJA::

        sage: OctonionHermitianEJA(3,                    # long time
        ....:                      field=QQ,             # long time
        ....:                      orthonormalize=False, # long time
        ....:                      check_axioms=True)    # long time
        Euclidean Jordan algebra of dimension 27 over Rational Field

    After a change-of-basis, the 2-by-2 algebra has the same
    multiplication table as the ten-dimensional Jordan spin algebra::

        sage: A = OctonionMatrixAlgebra(2,Octonions(QQ),QQ)
        sage: b = OctonionHermitianEJA._denormalized_basis(A)
        sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],)
        sage: jp = OctonionHermitianEJA.jordan_product
        sage: ip = OctonionHermitianEJA.trace_inner_product
        sage: J = FiniteDimensionalEJA(basis,
        ....:                          jp,
        ....:                          ip,
        ....:                          field=QQ,
        ....:                          orthonormalize=False)
        sage: J.multiplication_table()
        +----++----+----+----+----+----+----+----+----+----+----+
        | *  || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
        +====++====+====+====+====+====+====+====+====+====+====+
        | b0 || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
        +----++----+----+----+----+----+----+----+----+----+----+
        | b1 || b1 | b0 | 0  | 0  | 0  | 0  | 0  | 0  | 0  | 0  |
        +----++----+----+----+----+----+----+----+----+----+----+
        | b2 || b2 | 0  | b0 | 0  | 0  | 0  | 0  | 0  | 0  | 0  |
        +----++----+----+----+----+----+----+----+----+----+----+
        | b3 || b3 | 0  | 0  | b0 | 0  | 0  | 0  | 0  | 0  | 0  |
        +----++----+----+----+----+----+----+----+----+----+----+
        | b4 || b4 | 0  | 0  | 0  | b0 | 0  | 0  | 0  | 0  | 0  |
        +----++----+----+----+----+----+----+----+----+----+----+
        | b5 || b5 | 0  | 0  | 0  | 0  | b0 | 0  | 0  | 0  | 0  |
        +----++----+----+----+----+----+----+----+----+----+----+
        | b6 || b6 | 0  | 0  | 0  | 0  | 0  | b0 | 0  | 0  | 0  |
        +----++----+----+----+----+----+----+----+----+----+----+
        | b7 || b7 | 0  | 0  | 0  | 0  | 0  | 0  | b0 | 0  | 0  |
        +----++----+----+----+----+----+----+----+----+----+----+
        | b8 || b8 | 0  | 0  | 0  | 0  | 0  | 0  | 0  | b0 | 0  |
        +----++----+----+----+----+----+----+----+----+----+----+
        | b9 || b9 | 0  | 0  | 0  | 0  | 0  | 0  | 0  | 0  | b0 |
        +----++----+----+----+----+----+----+----+----+----+----+

    TESTS:

    We can actually construct the 27-dimensional Albert algebra,
    and we get the right unit element if we recompute it::

        sage: J = OctonionHermitianEJA(3,                    # long time
        ....:                          field=QQ,             # long time
        ....:                          orthonormalize=False) # long time
        sage: J.one.clear_cache()                            # long time
        sage: J.one()                                        # long time
        b0 + b9 + b26
        sage: J.one().to_matrix()                            # long time
        +----+----+----+
        | e0 | 0  | 0  |
        +----+----+----+
        | 0  | e0 | 0  |
        +----+----+----+
        | 0  | 0  | e0 |
        +----+----+----+

    The 2-by-2 algebra is isomorphic to the ten-dimensional Jordan
    spin algebra, but just to be sure, we recompute its rank::

        sage: J = OctonionHermitianEJA(2,                    # long time
        ....:                          field=QQ,             # long time
        ....:                          orthonormalize=False) # long time
        sage: J.rank.clear_cache()                           # long time
        sage: J.rank()                                       # long time
        2

    """
    @staticmethod
    def _max_random_instance_size(max_dimension):
        r"""
        The maximum rank of a random QuaternionHermitianEJA.
        """
        # There's certainly a formula for this, but with only four
        # cases to worry about, I'm not that motivated to derive it.
        if max_dimension >= 27:
            return 3
        elif max_dimension >= 10:
            return 2
        elif max_dimension >= 1:
            return 1
        else:
            return 0

    @classmethod
    def random_instance(cls, max_dimension=None, *args, **kwargs):
        """
        Return a random instance of this type of algebra.
        """
        class_max_d = cls._max_random_instance_dimension()
        if (max_dimension is None or max_dimension > class_max_d):
            max_dimension = class_max_d
        max_size = cls._max_random_instance_size(max_dimension)
        n = ZZ.random_element(max_size + 1)
        return cls(n, **kwargs)

    def __init__(self, n, field=AA, **kwargs):
        if n > 3:
            # Otherwise we don't get an EJA.
            raise ValueError("n cannot exceed 3")

        # We know this is a valid EJA, but will double-check
        # if the user passes check_axioms=True.
        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False

        from mjo.hurwitz import OctonionMatrixAlgebra
        A = OctonionMatrixAlgebra(n, scalars=field)
        super().__init__(A, **kwargs)

        from mjo.eja.eja_cache import octonion_hermitian_eja_coeffs
        a = octonion_hermitian_eja_coeffs(self)
        if a is not None:
            if self._rational_algebra is None:
                self._charpoly_coefficients.set_cache(a)
            else:
                self._rational_algebra._charpoly_coefficients.set_cache(a)


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

    """
    def __init__(self, *args, **kwargs):
        super().__init__(3, *args, **kwargs)


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

        sage: from mjo.eja.eja_algebra import HadamardEJA

    EXAMPLES:

    This multiplication table can be verified by hand::

        sage: J = HadamardEJA(3)
        sage: b0,b1,b2 = J.gens()
        sage: b0*b0
        b0
        sage: b0*b1
        0
        sage: b0*b2
        0
        sage: b1*b1
        b1
        sage: b1*b2
        0
        sage: b2*b2
        b2

    TESTS:

    We can change the generator prefix::

        sage: HadamardEJA(3, prefix='r').gens()
        (r0, r1, r2)
    """
    def __init__(self, n, field=AA, **kwargs):
        MS = MatrixSpace(field, n, 1)

        if n == 0:
            jordan_product = lambda x,y: x
            inner_product = lambda x,y: x
        else:
            def jordan_product(x,y):
                return MS( 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( MS(b) for b in FreeModule(field, n).basis() )
        super().__init__(column_basis,
                         jordan_product,
                         inner_product,
                         field=field,
                         matrix_space=MS,
                         associative=True,
                         **kwargs)
        self.rank.set_cache(n)

        self.one.set_cache( self.sum(self.gens()) )

    @staticmethod
    def _max_random_instance_dimension():
        r"""
        There's no reason to go higher than five here. That's
        enough to get the point across.
        """
        return 5

    @staticmethod
    def _max_random_instance_size(max_dimension):
        r"""
        The maximum size (=dimension) of a random HadamardEJA.
        """
        return max_dimension

    @classmethod
    def random_instance(cls, max_dimension=None, *args, **kwargs):
        """
        Return a random instance of this type of algebra.
        """
        class_max_d = cls._max_random_instance_dimension()
        if (max_dimension is None or max_dimension > class_max_d):
            max_dimension = class_max_d
        max_size = cls._max_random_instance_size(max_dimension)
        n = ZZ.random_element(max_size + 1)
        return cls(n, **kwargs)


class 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 =
    (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
    a symmetric positive-definite "bilinear form" matrix. Its
    dimension is the size of `B`, and it has rank two in dimensions
    larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
    the identity matrix of order ``n``.

    We insist that the one-by-one upper-left identity block of `B` be
    passed in as well so that we can be passed a matrix of size zero
    to construct a trivial algebra.

    SETUP::

        sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
        ....:                                  JordanSpinEJA)

    EXAMPLES:

    When no bilinear form is specified, the identity matrix is used,
    and the resulting algebra is the Jordan spin algebra::

        sage: B = matrix.identity(AA,3)
        sage: J0 = BilinearFormEJA(B)
        sage: J1 = JordanSpinEJA(3)
        sage: J0.multiplication_table() == J0.multiplication_table()
        True

    An error is raised if the matrix `B` does not correspond to a
    positive-definite bilinear form::

        sage: B = matrix.random(QQ,2,3)
        sage: J = BilinearFormEJA(B)
        Traceback (most recent call last):
        ...
        ValueError: bilinear form is not positive-definite
        sage: B = matrix.zero(QQ,3)
        sage: J = BilinearFormEJA(B)
        Traceback (most recent call last):
        ...
        ValueError: bilinear form is not positive-definite

    TESTS:

    We can create a zero-dimensional algebra::

        sage: B = matrix.identity(AA,0)
        sage: J = BilinearFormEJA(B)
        sage: J.basis()
        Finite family {}

    We can check the multiplication condition given in the Jordan, von
    Neumann, and Wigner paper (and also discussed on my "On the
    symmetry..." paper). Note that this relies heavily on the standard
    choice of basis, as does anything utilizing the bilinear form
    matrix.  We opt not to orthonormalize the basis, because if we
    did, we would have to normalize the `s_{i}` in a similar manner::

        sage: set_random_seed()
        sage: n = ZZ.random_element(5)
        sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
        sage: B11 = matrix.identity(QQ,1)
        sage: B22 = M.transpose()*M
        sage: B = block_matrix(2,2,[ [B11,0  ],
        ....:                        [0, B22 ] ])
        sage: J = BilinearFormEJA(B, orthonormalize=False)
        sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
        sage: V = J.vector_space()
        sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
        ....:         for ei in eis ]
        sage: actual = [ sis[i]*sis[j]
        ....:            for i in range(n-1)
        ....:            for j in range(n-1) ]
        sage: expected = [ J.one() if i == j else J.zero()
        ....:              for i in range(n-1)
        ....:              for j in range(n-1) ]
        sage: actual == expected
        True

    """
    def __init__(self, B, field=AA, **kwargs):
        # The matrix "B" is supplied by the user in most cases,
        # so it makes sense to check whether or not its positive-
        # definite unless we are specifically asked not to...
        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

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

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

        def jordan_product(x,y):
            x0 = x[0,0]
            xbar = x[1:,0]
            y0 = y[0,0]
            ybar = y[1:,0]
            z0 = inner_product(y,x)
            zbar = y0*xbar + x0*ybar
            return MS([z0] + zbar.list())

        column_basis = tuple( MS(b) for b in FreeModule(field, n).basis() )

        # TODO: I haven't actually checked this, but it seems legit.
        associative = False
        if n <= 2:
            associative = True

        super().__init__(column_basis,
                         jordan_product,
                         inner_product,
                         field=field,
                         matrix_space=MS,
                         associative=associative,
                         **kwargs)

        # The rank of this algebra is two, unless we're in a
        # one-dimensional ambient space (because the rank is bounded
        # by the ambient dimension).
        self.rank.set_cache(min(n,2))
        if n == 0:
            self.one.set_cache( self.zero() )
        else:
            self.one.set_cache( self.monomial(0) )

    @staticmethod
    def _max_random_instance_dimension():
        r"""
        There's no reason to go higher than five here. That's
        enough to get the point across.
        """
        return 5

    @staticmethod
    def _max_random_instance_size(max_dimension):
        r"""
        The maximum size (=dimension) of a random BilinearFormEJA.
        """
        return max_dimension

    @classmethod
    def random_instance(cls, max_dimension=None, *args, **kwargs):
        """
        Return a random instance of this algebra.
        """
        class_max_d = cls._max_random_instance_dimension()
        if (max_dimension is None or max_dimension > class_max_d):
            max_dimension = class_max_d
        max_size = cls._max_random_instance_size(max_dimension)
        n = ZZ.random_element(max_size + 1)

        if n.is_zero():
            B = matrix.identity(ZZ, n)
            return cls(B, **kwargs)

        B11 = matrix.identity(ZZ, 1)
        M = matrix.random(ZZ, n-1)
        I = matrix.identity(ZZ, n-1)
        alpha = ZZ.zero()
        while alpha.is_zero():
            alpha = ZZ.random_element().abs()

        B22 = M.transpose()*M + alpha*I

        from sage.matrix.special import block_matrix
        B = block_matrix(2,2, [ [B11,   ZZ(0) ],
                                [ZZ(0), B22 ] ])

        return cls(B, **kwargs)


class JordanSpinEJA(BilinearFormEJA):
    """
    The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
    with the usual inner product and jordan product ``x*y =
    (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
    the reals.

    SETUP::

        sage: from mjo.eja.eja_algebra import JordanSpinEJA

    EXAMPLES:

    This multiplication table can be verified by hand::

        sage: J = JordanSpinEJA(4)
        sage: b0,b1,b2,b3 = J.gens()
        sage: b0*b0
        b0
        sage: b0*b1
        b1
        sage: b0*b2
        b2
        sage: b0*b3
        b3
        sage: b1*b2
        0
        sage: b1*b3
        0
        sage: b2*b3
        0

    We can change the generator prefix::

        sage: JordanSpinEJA(2, prefix='B').gens()
        (B0, B1)

    TESTS:

        Ensure that we have the usual inner product on `R^n`::

            sage: set_random_seed()
            sage: J = JordanSpinEJA.random_instance()
            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

    """
    def __init__(self, n, *args, **kwargs):
        # This is a special case of the BilinearFormEJA with the
        # identity matrix as its bilinear form.
        B = matrix.identity(ZZ, n)

        # Don't orthonormalize because our basis is already
        # orthonormal with respect to our inner-product.
        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)

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

        Needed here to override the implementation for ``BilinearFormEJA``.
        """
        class_max_d = cls._max_random_instance_dimension()
        if (max_dimension is None or max_dimension > class_max_d):
            max_dimension = class_max_d
        max_size = cls._max_random_instance_size(max_dimension)
        n = ZZ.random_element(max_size + 1)
        return cls(n, **kwargs)


class TrivialEJA(RationalBasisEJA, ConcreteEJA):
    """
    The trivial Euclidean Jordan algebra consisting of only a zero element.

    SETUP::

        sage: from mjo.eja.eja_algebra import TrivialEJA

    EXAMPLES::

        sage: J = TrivialEJA()
        sage: J.dimension()
        0
        sage: J.zero()
        0
        sage: J.one()
        0
        sage: 7*J.one()*12*J.one()
        0
        sage: J.one().inner_product(J.one())
        0
        sage: J.one().norm()
        0
        sage: J.one().subalgebra_generated_by()
        Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
        sage: J.rank()
        0

    """
    def __init__(self, field=AA, **kwargs):
        jordan_product = lambda x,y: x
        inner_product = lambda x,y: field.zero()
        basis = ()
        MS = MatrixSpace(field,0)

        # 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,
                         field=field,
                         matrix_space=MS,
                         **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
    def random_instance(cls, max_dimension=None, *args, **kwargs):
        # We don't take a "size" argument so the superclass method is
        # inappropriate for us. The ``max_dimension`` argument is
        # included so that if this method is called generically with a
        # ``max_dimension=<whatever>`` argument, we don't try to pass
        # it on to the algebra constructor.
        return cls(**kwargs)


class CartesianProductEJA(FiniteDimensionalEJA):
    r"""
    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,
        ....:                                  CartesianProductEJA,
        ....:                                  HadamardEJA,
        ....:                                  JordanSpinEJA,
        ....:                                  RealSymmetricEJA)

    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: 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()
        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:

    All factors must share the same base field::

        sage: J1 = HadamardEJA(2, field=QQ)
        sage: J2 = RealSymmetricEJA(2)
        sage: CartesianProductEJA((J1,J2))
        Traceback (most recent call last):
        ...
        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

    """
    Element = FiniteDimensionalEJAElement


    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 )

        # Compute my matrix space. This category isn't perfect, but
        # is good enough for what we need to do.
        MS_cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
        MS_cat = MS_cat.Unital().CartesianProducts()
        MS_factors = tuple( J.matrix_space() for J in factors )
        from sage.sets.cartesian_product import CartesianProduct
        MS = CartesianProduct(MS_factors, MS_cat)

        basis = []
        zero = MS.zero()
        for i in range(m):
            for b in factors[i].matrix_basis():
                z = list(zero)
                z[i] = b
                basis.append(z)

        basis = tuple( MS(b) for b in basis )

        # Define jordan/inner products that operate on that matrix_basis.
        def jordan_product(x,y):
            return MS(tuple(
                (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,
                                      matrix_space=MS,
                                      orthonormalize=False,
                                      associative=associative,
                                      cartesian_product=True,
                                      check_field=False,
                                      check_axioms=False)

        # Since we don't (re)orthonormalize the basis, the FDEJA
        # constructor is going to set self._deortho_matrix to the
        # identity matrix. Here we set it to the correct value using
        # the deortho matrices from our factors.
        self._deortho_matrix = matrix.block_diagonal( [J._deortho_matrix
                                                       for J in factors] )

        self.rank.set_cache(sum(J.rank() for J in factors))
        ones = tuple(J.one().to_matrix() for J in factors)
        self.one.set_cache(self(ones))

    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"""
        Return the space that our matrix basis lives in as a Cartesian
        product.

        We don't simply use the ``cartesian_product()`` functor here
        because it acts differently on SageMath MatrixSpaces and our
        custom MatrixAlgebras, which are CombinatorialFreeModules. We
        always want the result to be represented (and indexed) as
        an ordered tuple.

        SETUP::

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

        EXAMPLES::

            sage: J1 = HadamardEJA(1)
            sage: J2 = RealSymmetricEJA(2)
            sage: J = cartesian_product([J1,J2])
            sage: J.matrix_space()
            The Cartesian product of (Full MatrixSpace of 1 by 1 dense
            matrices over Algebraic Real Field, Full MatrixSpace of 2
            by 2 dense matrices over Algebraic Real Field)

        ::

            sage: J1 = ComplexHermitianEJA(1)
            sage: J2 = ComplexHermitianEJA(1)
            sage: J = cartesian_product([J1,J2])
            sage: J.one().to_matrix()[0]
            +---+
            | 1 |
            +---+
            sage: J.one().to_matrix()[1]
            +---+
            | 1 |
            +---+

        ::

            sage: J1 = OctonionHermitianEJA(1)
            sage: J2 = OctonionHermitianEJA(1)
            sage: J = cartesian_product([J1,J2])
            sage: J.one().to_matrix()[0]
            +----+
            | e0 |
            +----+
            sage: J.one().to_matrix()[1]
            +----+
            | e0 |
            +----+

        """
        return super().matrix_space()


    @cached_method
    def cartesian_projection(self, i):
        r"""
        SETUP::

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

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

        """
        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,
            ....:                                  HadamardEJA,
            ....:                                  RealSymmetricEJA)

        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: 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()
            True
            sage: J1.one().to_vector()
            (1, 0, 0)
            sage: iota_left(J1.one()).to_vector()
            (1, 0, 0, 0, 0, 0)
            sage: J2.one().to_vector()
            (1, 0, 1)
            sage: iota_right(J2.one()).to_vector()
            (0, 0, 0, 1, 0, 1)
            sage: J.one().to_vector()
            (1, 0, 0, 1, 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: 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::

            sage: set_random_seed()
            sage: J1 = random_eja()
            sage: J2 = random_eja()
            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()
            True
            sage: (pi_left*iota_right).is_zero()
            True
            sage: (pi_right*iota_left).is_zero()
            True

        """
        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())



FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA

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

    SETUP::

        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

    TESTS:

    The ``cartesian_product()`` function only uses the first factor to
    decide where the result will live; thus we have to be careful to
    check that all factors do indeed have a `_rational_algebra` member
    before we try to access it::

        sage: J1 = OctonionHermitianEJA(1) # no rational basis
        sage: J2 = HadamardEJA(2)
        sage: cartesian_product([J1,J2])
        Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
        (+) Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
        sage: cartesian_product([J2,J1])
        Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
        (+) Euclidean Jordan algebra of dimension 1 over Algebraic Real Field

    """
    def __init__(self, algebras, **kwargs):
        CartesianProductEJA.__init__(self, algebras, **kwargs)

        self._rational_algebra = None
        if self.vector_space().base_field() is not QQ:
            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(max_dimension=None, *args, **kwargs):
    r"""

    SETUP::

        sage: from mjo.eja.eja_algebra import random_eja

    TESTS::

        sage: set_random_seed()
        sage: n = ZZ.random_element(1,5)
        sage: J = random_eja(max_dimension=n, field=QQ, orthonormalize=False)
        sage: J.dimension() <= n
        True

    """
    # Use the ConcreteEJA default as the total upper bound (regardless
    # of any whether or not any individual factors set a lower limit).
    if max_dimension is None:
        max_dimension = ConcreteEJA._max_random_instance_dimension()
    J1 = ConcreteEJA.random_instance(max_dimension, *args, **kwargs)


    # Roll the dice to see if we attempt a Cartesian product.
    dice_roll = ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1)
    new_max_dimension = max_dimension - J1.dimension()
    if new_max_dimension == 0 or dice_roll != 0:
        # If it's already as big as we're willing to tolerate, just
        # return it and don't worry about Cartesian products.
        return J1
    else:
        # Use random_eja() again so we can get more than two factors
        # if the sub-call also Decides on a cartesian product.
        J2 = random_eja(new_max_dimension, *args, **kwargs)
        return cartesian_product([J1,J2])