eja: move the "field" argument to (usually passed through) kwargs.

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

"""
Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
are used in optimization, and have some additional nice methods beyond
what can be supported in a general Jordan Algebra.


SETUP::

    sage: from mjo.eja.eja_algebra import random_eja

EXAMPLES::

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

"""

from itertools import repeat

from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
from sage.categories.magmatic_algebras import MagmaticAlgebras
from sage.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 FiniteDimensionalEuclideanJordanAlgebraElement
from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
from mjo.eja.eja_utils import _mat2vec

class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
    r"""
    The lowest-level class for representing a Euclidean Jordan algebra.
    """
    def _coerce_map_from_base_ring(self):
        """
        Disable the map from the base ring into the algebra.

        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 __init__(self,
                 field,
                 multiplication_table,
                 inner_product_table,
                 prefix='e',
                 category=None,
                 matrix_basis=None,
                 check_field=True,
                 check_axioms=True):
        """
        INPUT:

          * field -- the scalar field for this algebra (must be real)

          * multiplication_table -- the multiplication table for this
            algebra's implicit basis. Only the lower-triangular portion
            of the table is used, since the multiplication is assumed
            to be commutative.

        SETUP::

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

        EXAMPLES:

        By definition, Jordan multiplication commutes::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: x,y = J.random_elements(2)
            sage: x*y == y*x
            True

        An error is raised if the Jordan product is not commutative::

            sage: JP = ((1,2),(0,0))
            sage: IP = ((1,0),(0,1))
            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
            Traceback (most recent call last):
            ...
            ValueError: Jordan product is not commutative

        An error is raised if the inner-product is not commutative::

            sage: JP = ((1,0),(0,1))
            sage: IP = ((1,2),(0,0))
            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
            Traceback (most recent call last):
            ...
            ValueError: inner-product is not commutative

        TESTS:

        The ``field`` we're given must be real with ``check_field=True``::

            sage: JordanSpinEJA(2, field=QQbar)
            Traceback (most recent call last):
            ...
            ValueError: scalar field is not real
            sage: JordanSpinEJA(2, field=QQbar, check_field=False)
            Euclidean Jordan algebra of dimension 2 over Algebraic Field

        The multiplication table must be square with ``check_axioms=True``::

            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()),((1,),))
            Traceback (most recent call last):
            ...
            ValueError: multiplication table is not square

        The multiplication and inner-product tables must be the same
        size (and in particular, the inner-product table must also be
        square) with ``check_axioms=True``::

            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),(()))
            Traceback (most recent call last):
            ...
            ValueError: multiplication and inner-product tables are
            different sizes
            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),((1,2),))
            Traceback (most recent call last):
            ...
            ValueError: multiplication and inner-product tables are
            different sizes

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


        # The multiplication and inner-product tables should be square
        # if the user wants us to verify them. And we verify them as
        # soon as possible, because we want to exploit their symmetry.
        n = len(multiplication_table)
        if check_axioms:
            if not all( len(l) == n for l in multiplication_table ):
                raise ValueError("multiplication table is not square")

            # If the multiplication table is square, we can check if
            # the inner-product table is square by comparing it to the
            # multiplication table's dimensions.
            msg = "multiplication and inner-product tables are different sizes"
            if not len(inner_product_table) == n:
                raise ValueError(msg)

            if not all( len(l) == n for l in inner_product_table ):
                raise ValueError(msg)

            # Check commutativity of the Jordan product (symmetry of
            # the multiplication table) and the commutativity of the
            # inner-product (symmetry of the inner-product table)
            # first if we're going to check them at all.. This has to
            # be done before we define product_on_basis(), because
            # that method assumes that self._multiplication_table is
            # symmetric. And it has to be done before we build
            # self._inner_product_matrix, because the process used to
            # construct it assumes symmetry as well.
            if not all(    multiplication_table[j][i]
                        == multiplication_table[i][j]
                        for i in range(n)
                        for j in range(i+1) ):
                raise ValueError("Jordan product is not commutative")

            if not all( inner_product_table[j][i]
                        == inner_product_table[i][j]
                        for i in range(n)
                        for j in range(i+1) ):
                raise ValueError("inner-product is not commutative")

        self._matrix_basis = matrix_basis

        if category is None:
            category = MagmaticAlgebras(field).FiniteDimensional()
            category = category.WithBasis().Unital()

        fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
        fda.__init__(field,
                     range(n),
                     prefix=prefix,
                     category=category)
        self.print_options(bracket='')

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

        for i in range(n):
            for j in range(i+1):
                elt = self.from_vector(multiplication_table[i][j])
                self._multiplication_table[i][j] = elt

        self._multiplication_table = tuple(map(tuple, self._multiplication_table))

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

        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 _element_constructor_(self, elt):
        """
        Construct an element of this algebra from its vector or matrix
        representation.

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

        SETUP::

            sage: from mjo.eja.eja_algebra import (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

        TESTS:

        Ensure that we can convert any element of the two non-matrix
        simple algebras (whose matrix representations are columns)
        back and forth faithfully::

            sage: set_random_seed()
            sage: J = HadamardEJA.random_instance()
            sage: x = J.random_element()
            sage: J(x.to_vector().column()) == x
            True
            sage: J = JordanSpinEJA.random_instance()
            sage: x = J.random_element()
            sage: J(x.to_vector().column()) == x
            True

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

        if elt not in self.matrix_space():
            raise ValueError(msg)

        # 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.
        #
        # We pass check=False because the matrix basis is "guaranteed"
        # to be linearly independent... right? Ha ha.
        V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols())
        W = V.span_of_basis( (_mat2vec(s) for s in self.matrix_basis()),
                             check=False)

        try:
            coords =  W.coordinate_vector(_mat2vec(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())

    def product_on_basis(self, i, j):
        # We only stored the lower-triangular portion of the
        # multiplication table.
        if j <= i:
            return self._multiplication_table[i][j]
        else:
            return self._multiplication_table[j][i]

    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( self.product_on_basis(i,j) == self.product_on_basis(i,j)
                    for i in range(self.dimension())
                    for j in range(self.dimension()) )

    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 _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 multiplication table.
        """

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

        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 self.base_ring().is_exact():
                        if diff != 0:
                            return False
                    else:
                        if diff.abs() > epsilon:
                            return False

        return True

    @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()
            +----++----+----+----+----+
            | *  || e0 | e1 | e2 | e3 |
            +====++====+====+====+====+
            | e0 || e0 | e1 | e2 | e3 |
            +----++----+----+----+----+
            | e1 || e1 | e0 | 0  | 0  |
            +----++----+----+----+----+
            | e2 || e2 | 0  | e0 | 0  |
            +----++----+----+----+----+
            | e3 || e3 | 0  | 0  | e0 |
            +----++----+----+----+----+

        """
        n = self.dimension()
        M = [ [ self.zero() for j in range(n) ]
              for i in range(n) ]
        for i in range(n):
            for j in range(i+1):
                M[i][j] = self._multiplication_table[i][j]
                M[j][i] = M[i][j]

        for i in range(n):
            # Prepend the left "header" column entry Can't do this in
            # the loop because it messes up the symmetry.
            M[i] = [self.monomial(i)] + M[i]

        # Prepend the header row.
        M = [["*"] + list(self.gens())] + M
        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:`DirectSumEJA` 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: e0, 1: e1, 2: e2}
            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: e0, 1: e1}
            sage: J.matrix_basis()
            (
            [1]  [0]
            [0], [1]
            )
        """
        if self._matrix_basis is None:
            M = self.matrix_space()
            return tuple( M(b.to_vector()) for b in self.basis() )
        else:
            return self._matrix_basis


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

        Generally this will be an `n`-by-`1` column-vector 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.

        Matrix algebras override this with something more useful.
        """
        if self.is_trivial():
            return MatrixSpace(self.base_ring(), 0)
        elif self._matrix_basis is None or len(self._matrix_basis) == 0:
            return MatrixSpace(self.base_ring(), self.dimension(), 1)
        else:
            return self._matrix_basis[0].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::

            sage: J = HadamardEJA(5)
            sage: J.one()
            e0 + e1 + e2 + e3 + e4

        TESTS:

        The identity element acts like the identity::

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

        The matrix of the unit element's operator is the identity::

            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

        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

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

        from mjo.eja.eja_subalgebra import FiniteDimensionalEuclideanJordanSubalgebra

        # Default these to what they should be if they turn out to be
        # trivial, because eigenspaces_left() won't return eigenvalues
        # corresponding to trivial spaces (e.g. it returns only the
        # eigenspace corresponding to lambda=1 if you take the
        # decomposition relative to the identity element).
        trivial = FiniteDimensionalEuclideanJordanSubalgebra(self, ())
        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 = FiniteDimensionalEuclideanJordanSubalgebra(self,
                                                                    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.
        """
        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.
        return -A_rref.solve_right(E*b).change_ring(R)[: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: caches = J.rank()    # long time
            sage: J.rank.clear_cache() # long time
            sage: J.rank() == cached   # long time
            True

        """
        return len(self._charpoly_coefficients())


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


    Element = FiniteDimensionalEuclideanJordanAlgebraElement

class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
    r"""
    New class for 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,
                 orthonormalize=True,
                 prefix='e',
                 category=None,
                 check_field=True,
                 check_axioms=True):

        if check_field:
            # Abuse the check_field parameter to check that the entries of
            # out basis (in ambient coordinates) are in the field QQ.
            if not all( all(b_i in QQ for b_i in b.list()) for b in basis ):
                raise TypeError("basis not rational")

        # Temporary(?) hack to ensure that the matrix and vector bases
        # are over the same ring.
        basis = tuple( b.change_ring(field) for b in basis )

        n = len(basis)
        vector_basis = basis

        from sage.structure.element import is_Matrix
        basis_is_matrices = False

        degree = 0
        if n > 0:
            if is_Matrix(basis[0]):
                basis_is_matrices = True
                from mjo.eja.eja_utils import _vec2mat
                vector_basis = tuple( map(_mat2vec,basis) )
                degree = basis[0].nrows()**2
            else:
                degree = basis[0].degree()

        V = VectorSpace(field, degree)

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

        if orthonormalize:
            if self.base_ring() is not QQ:
                # There's no point in constructing the extra algebra if this
                # one is already rational. If the original basis is rational
                # but normalization would make it irrational, then this whole
                # constructor will just fail anyway as it tries to stick an
                # irrational number into a rational algebra.
                #
                # Note: the same Jordan and inner-products work here,
                # because they are necessarily defined with respect to
                # ambient coordinates and not any particular basis.
                self._rational_algebra = RationalBasisEuclideanJordanAlgebra(
                                          basis,
                                          jordan_product,
                                          inner_product,
                                          field=QQ,
                                          orthonormalize=False,
                                          prefix=prefix,
                                          category=category,
                                          check_field=False,
                                          check_axioms=False)

            # Compute the deorthonormalized tables before we orthonormalize
            # the given basis. The "check" parameter here guarantees that
            # the basis is linearly-independent.
            W = V.span_of_basis( vector_basis, check=check_axioms)

            # Note: the Jordan and inner-products are defined in terms
            # of the ambient basis. It's important that their arguments
            # are in ambient coordinates as well.
            for i in range(n):
                for j in range(i+1):
                    # given basis w.r.t. ambient coords
                    q_i = vector_basis[i]
                    q_j = vector_basis[j]

                    if basis_is_matrices:
                        q_i = _vec2mat(q_i)
                        q_j = _vec2mat(q_j)

                    elt = jordan_product(q_i, q_j)
                    ip = inner_product(q_i, q_j)

                    if basis_is_matrices:
                        # do another mat2vec because the multiplication
                        # table is in terms of vectors
                        elt = _mat2vec(elt)

        # We overwrite the name "vector_basis" in a second, but never modify it
        # in place, to this effectively makes a copy of it.
        deortho_vector_basis = vector_basis
        self._deortho_matrix = None

        if orthonormalize:
            from mjo.eja.eja_utils import gram_schmidt
            if basis_is_matrices:
                vector_ip = lambda x,y: inner_product(_vec2mat(x), _vec2mat(y))
                vector_basis = gram_schmidt(vector_basis, vector_ip)
            else:
                vector_basis = gram_schmidt(vector_basis, inner_product)

            # Normalize the "matrix" basis, too!
            basis = vector_basis

            if basis_is_matrices:
                basis = tuple( map(_vec2mat,basis) )

        W = V.span_of_basis( vector_basis, check=check_axioms)

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

        # If the superclass constructor is going to verify the
        # symmetry of this table, it has better at least be
        # square...
        if check_axioms:
            mult_table = [ [0 for j in range(n)] for i in range(n) ]
            ip_table = [ [0 for j in range(n)] for i in range(n) ]
        else:
            mult_table = [ [0 for j in range(i+1)] for i in range(n) ]
            ip_table = [ [0 for j in range(i+1)] for i in range(n) ]

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

                if basis_is_matrices:
                    q_i = _vec2mat(q_i)
                    q_j = _vec2mat(q_j)

                elt = jordan_product(q_i, q_j)
                ip = inner_product(q_i, q_j)

                if basis_is_matrices:
                    # do another mat2vec because the multiplication
                    # table is in terms of vectors
                    elt = _mat2vec(elt)

                elt = W.coordinate_vector(elt)
                mult_table[i][j] = elt
                ip_table[i][j] = ip
                if check_axioms:
                    # The tables are square if we're verifying that they
                    # are commutative.
                    mult_table[j][i] = elt
                    ip_table[j][i] = ip

        if basis_is_matrices:
            for m in basis:
                m.set_immutable()
        else:
            basis = tuple( x.column() for x in basis )

        super().__init__(field,
                         mult_table,
                         ip_table,
                         prefix,
                         category,
                         basis, # matrix basis
                         check_field,
                         check_axioms)

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

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

        EXAMPLES:

        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.base_ring() is QQ or 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.
            superclass = super(RationalBasisEuclideanJordanAlgebra, self)
            return superclass._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() )

        # Now 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 ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra):
    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 ConcreteEuclideanJordanAlgebra

    TESTS:

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

        sage: set_random_seed()
        sage: J = ConcreteEuclideanJordanAlgebra.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 = ConcreteEuclideanJordanAlgebra.random_instance()
        sage: x = J.random_element()
        sage: x.operator().is_self_adjoint()
        True
    """

    @staticmethod
    def _max_random_instance_size():
        """
        Return an integer "size" that is an upper bound on the size of
        this algebra when it is used in a random test
        case. Unfortunately, the term "size" is ambiguous -- when
        dealing with `R^n` under either the Hadamard or Jordan spin
        product, the "size" refers to the dimension `n`. When dealing
        with a matrix algebra (real symmetric or complex/quaternion
        Hermitian), it refers to the size of the matrix, which is far
        less than the dimension of the underlying vector space.

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

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

        This method should be implemented in each subclass.
        """
        from sage.misc.prandom import choice
        eja_class = choice(cls.__subclasses__())

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


class MatrixEuclideanJordanAlgebra:
    @staticmethod
    def real_embed(M):
        """
        Embed the matrix ``M`` into a space of real matrices.

        The matrix ``M`` can have entries in any field at the moment:
        the real numbers, complex numbers, or quaternions. And although
        they are not a field, we can probably support octonions at some
        point, too. This function returns a real matrix that "acts like"
        the original with respect to matrix multiplication; i.e.

          real_embed(M*N) = real_embed(M)*real_embed(N)

        """
        raise NotImplementedError


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

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

    @classmethod
    def trace_inner_product(cls,X,Y):
        Xu = cls.real_unembed(X)
        Yu = cls.real_unembed(Y)
        tr = (Xu*Yu).trace()

        try:
            # Works in QQ, AA, RDF, et cetera.
            return tr.real()
        except AttributeError:
            # A quaternion doesn't have a real() method, but does
            # have coefficient_tuple() method that returns the
            # coefficients of 1, i, j, and k -- in that order.
            return tr.coefficient_tuple()[0]


class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
    @staticmethod
    def real_embed(M):
        """
        The identity function, for embedding real matrices into real
        matrices.
        """
        return M

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


class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra,
                       RealMatrixEuclideanJordanAlgebra):
    """
    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: e0, e1, e2 = J.gens()
        sage: e0*e0
        e0
        sage: e1*e1
        1/2*e0 + 1/2*e2
        sage: e2*e2
        e2

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

        sage: RealSymmetricEJA(2, field=RDF)
        Euclidean Jordan algebra of dimension 3 over Real Double Field
        sage: RealSymmetricEJA(2, field=RR)
        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: n_max = RealSymmetricEJA._max_random_instance_size()
        sage: n = ZZ.random_element(1, n_max)
        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

    """
    @classmethod
    def _denormalized_basis(cls, n):
        """
        Return a basis for the space of real symmetric n-by-n matrices.

        SETUP::

            sage: from mjo.eja.eja_algebra import RealSymmetricEJA

        TESTS::

            sage: set_random_seed()
            sage: n = ZZ.random_element(1,5)
            sage: B = RealSymmetricEJA._denormalized_basis(n)
            sage: all( M.is_symmetric() for M in  B)
            True

        """
        # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
        # coordinates.
        S = []
        for i in range(n):
            for j in range(i+1):
                Eij = matrix(ZZ, n, lambda k,l: k==i and l==j)
                if i == j:
                    Sij = Eij
                else:
                    Sij = Eij + Eij.transpose()
                S.append(Sij)
        return tuple(S)


    @staticmethod
    def _max_random_instance_size():
        return 4 # Dimension 10

    @classmethod
    def random_instance(cls, **kwargs):
        """
        Return a random instance of this type of algebra.
        """
        n = ZZ.random_element(cls._max_random_instance_size() + 1)
        return cls(n, **kwargs)

    def __init__(self, n, **kwargs):
        super(RealSymmetricEJA, self).__init__(self._denormalized_basis(n),
                                               self.jordan_product,
                                               self.trace_inner_product,
                                               **kwargs)
        self.rank.set_cache(n)
        self.one.set_cache(self(matrix.identity(ZZ,n)))


class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
    @staticmethod
    def real_embed(M):
        """
        Embed the n-by-n complex matrix ``M`` into the space of real
        matrices of size 2n-by-2n via the map the sends each entry `z = a +
        bi` to the block matrix ``[[a,b],[-b,a]]``.

        SETUP::

            sage: from mjo.eja.eja_algebra import \
            ....:   ComplexMatrixEuclideanJordanAlgebra

        EXAMPLES::

            sage: F = QuadraticField(-1, 'I')
            sage: x1 = F(4 - 2*i)
            sage: x2 = F(1 + 2*i)
            sage: x3 = F(-i)
            sage: x4 = F(6)
            sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
            sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
            [ 4 -2| 1  2]
            [ 2  4|-2  1]
            [-----+-----]
            [ 0 -1| 6  0]
            [ 1  0| 0  6]

        TESTS:

        Embedding is a homomorphism (isomorphism, in fact)::

            sage: set_random_seed()
            sage: n = ZZ.random_element(3)
            sage: F = QuadraticField(-1, 'I')
            sage: X = random_matrix(F, n)
            sage: Y = random_matrix(F, n)
            sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
            sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
            sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
            sage: Xe*Ye == XYe
            True

        """
        n = M.nrows()
        if M.ncols() != n:
            raise ValueError("the matrix 'M' must be square")

        # We don't need any adjoined elements...
        field = M.base_ring().base_ring()

        blocks = []
        for z in M.list():
            a = z.list()[0] # real part, I guess
            b = z.list()[1] # imag part, I guess
            blocks.append(matrix(field, 2, [[a,b],[-b,a]]))

        return matrix.block(field, n, blocks)


    @staticmethod
    def real_unembed(M):
        """
        The inverse of _embed_complex_matrix().

        SETUP::

            sage: from mjo.eja.eja_algebra import \
            ....:   ComplexMatrixEuclideanJordanAlgebra

        EXAMPLES::

            sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
            ....:                 [-2,  1,  -4,  3],
            ....:                 [ 9,  10, 11, 12],
            ....:                 [-10, 9, -12, 11] ])
            sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
            [  2*I + 1   4*I + 3]
            [ 10*I + 9 12*I + 11]

        TESTS:

        Unembedding is the inverse of embedding::

            sage: set_random_seed()
            sage: F = QuadraticField(-1, 'I')
            sage: M = random_matrix(F, 3)
            sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
            sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
            True

        """
        n = ZZ(M.nrows())
        if M.ncols() != n:
            raise ValueError("the matrix 'M' must be square")
        if not n.mod(2).is_zero():
            raise ValueError("the matrix 'M' must be a complex embedding")

        # If "M" was normalized, its base ring might have roots
        # adjoined and they can stick around after unembedding.
        field = M.base_ring()
        R = PolynomialRing(field, 'z')
        z = R.gen()
        if field is AA:
            # Sage doesn't know how to embed AA into QQbar, i.e. how
            # to adjoin sqrt(-1) to AA.
            F = QQbar
        else:
            F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
        i = F.gen()

        # Go top-left to bottom-right (reading order), converting every
        # 2-by-2 block we see to a single complex element.
        elements = []
        for k in range(n/2):
            for j in range(n/2):
                submat = M[2*k:2*k+2,2*j:2*j+2]
                if submat[0,0] != submat[1,1]:
                    raise ValueError('bad on-diagonal submatrix')
                if submat[0,1] != -submat[1,0]:
                    raise ValueError('bad off-diagonal submatrix')
                z = submat[0,0] + submat[0,1]*i
                elements.append(z)

        return matrix(F, n/2, elements)


    @classmethod
    def trace_inner_product(cls,X,Y):
        """
        Compute a matrix inner product in this algebra directly from
        its real embedding.

        SETUP::

            sage: from mjo.eja.eja_algebra import ComplexHermitianEJA

        TESTS:

        This gives the same answer as the slow, default method implemented
        in :class:`MatrixEuclideanJordanAlgebra`::

            sage: set_random_seed()
            sage: J = ComplexHermitianEJA.random_instance()
            sage: x,y = J.random_elements(2)
            sage: Xe = x.to_matrix()
            sage: Ye = y.to_matrix()
            sage: X = ComplexHermitianEJA.real_unembed(Xe)
            sage: Y = ComplexHermitianEJA.real_unembed(Ye)
            sage: expected = (X*Y).trace().real()
            sage: actual = ComplexHermitianEJA.trace_inner_product(Xe,Ye)
            sage: actual == expected
            True

        """
        return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/2


class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra,
                          ComplexMatrixEuclideanJordanAlgebra):
    """
    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::

        sage: ComplexHermitianEJA(2, field=RDF)
        Euclidean Jordan algebra of dimension 4 over Real Double Field
        sage: ComplexHermitianEJA(2, field=RR)
        Euclidean Jordan algebra of dimension 4 over Real Field with
        53 bits of precision

    TESTS:

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

        sage: set_random_seed()
        sage: n_max = ComplexHermitianEJA._max_random_instance_size()
        sage: n = ZZ.random_element(1, n_max)
        sage: J = ComplexHermitianEJA(n)
        sage: J.dimension() == n^2
        True

    The Jordan multiplication is what we think it is::

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

    We can change the generator prefix::

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

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

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

    """

    @classmethod
    def _denormalized_basis(cls, n):
        """
        Returns a basis for the space of complex Hermitian n-by-n matrices.

        Why do we embed these? Basically, because all of numerical linear
        algebra assumes that you're working with vectors consisting of `n`
        entries from a field and scalars from the same field. There's no way
        to tell SageMath that (for example) the vectors contain complex
        numbers, while the scalar field is real.

        SETUP::

            sage: from mjo.eja.eja_algebra import ComplexHermitianEJA

        TESTS::

            sage: set_random_seed()
            sage: n = ZZ.random_element(1,5)
            sage: field = QuadraticField(2, 'sqrt2')
            sage: B = ComplexHermitianEJA._denormalized_basis(n)
            sage: all( M.is_symmetric() for M in  B)
            True

        """
        field = ZZ
        R = PolynomialRing(field, 'z')
        z = R.gen()
        F = field.extension(z**2 + 1, 'I')
        I = F.gen(1)

        # This is like the symmetric case, but we need to be careful:
        #
        #   * We want conjugate-symmetry, not just symmetry.
        #   * The diagonal will (as a result) be real.
        #
        S = []
        for i in range(n):
            for j in range(i+1):
                Eij = matrix(F, n, lambda k,l: k==i and l==j)
                if i == j:
                    Sij = cls.real_embed(Eij)
                    S.append(Sij)
                else:
                    # The second one has a minus because it's conjugated.
                    Sij_real = cls.real_embed(Eij + Eij.transpose())
                    S.append(Sij_real)
                    Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose())
                    S.append(Sij_imag)

        # Since we embedded these, we can drop back to the "field" that we
        # started with instead of the complex extension "F".
        return tuple( s.change_ring(field) for s in S )


    def __init__(self, n, **kwargs):
        super(ComplexHermitianEJA, self).__init__(self._denormalized_basis(n),
                                                  self.jordan_product,
                                                  self.trace_inner_product,
                                                  **kwargs)
        self.rank.set_cache(n)
        # TODO: pre-cache the identity!

    @staticmethod
    def _max_random_instance_size():
        return 3 # Dimension 9

    @classmethod
    def random_instance(cls, **kwargs):
        """
        Return a random instance of this type of algebra.
        """
        n = ZZ.random_element(cls._max_random_instance_size() + 1)
        return cls(n, **kwargs)

class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
    @staticmethod
    def real_embed(M):
        """
        Embed the n-by-n quaternion matrix ``M`` into the space of real
        matrices of size 4n-by-4n by first sending each quaternion entry `z
        = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
        c+di],[-c + di, a-bi]]`, and then embedding those into a real
        matrix.

        SETUP::

            sage: from mjo.eja.eja_algebra import \
            ....:   QuaternionMatrixEuclideanJordanAlgebra

        EXAMPLES::

            sage: Q = QuaternionAlgebra(QQ,-1,-1)
            sage: i,j,k = Q.gens()
            sage: x = 1 + 2*i + 3*j + 4*k
            sage: M = matrix(Q, 1, [[x]])
            sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
            [ 1  2  3  4]
            [-2  1 -4  3]
            [-3  4  1 -2]
            [-4 -3  2  1]

        Embedding is a homomorphism (isomorphism, in fact)::

            sage: set_random_seed()
            sage: n = ZZ.random_element(2)
            sage: Q = QuaternionAlgebra(QQ,-1,-1)
            sage: X = random_matrix(Q, n)
            sage: Y = random_matrix(Q, n)
            sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
            sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
            sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
            sage: Xe*Ye == XYe
            True

        """
        quaternions = M.base_ring()
        n = M.nrows()
        if M.ncols() != n:
            raise ValueError("the matrix 'M' must be square")

        F = QuadraticField(-1, 'I')
        i = F.gen()

        blocks = []
        for z in M.list():
            t = z.coefficient_tuple()
            a = t[0]
            b = t[1]
            c = t[2]
            d = t[3]
            cplxM = matrix(F, 2, [[ a + b*i, c + d*i],
                                 [-c + d*i, a - b*i]])
            realM = ComplexMatrixEuclideanJordanAlgebra.real_embed(cplxM)
            blocks.append(realM)

        # We should have real entries by now, so use the realest field
        # we've got for the return value.
        return matrix.block(quaternions.base_ring(), n, blocks)



    @staticmethod
    def real_unembed(M):
        """
        The inverse of _embed_quaternion_matrix().

        SETUP::

            sage: from mjo.eja.eja_algebra import \
            ....:   QuaternionMatrixEuclideanJordanAlgebra

        EXAMPLES::

            sage: M = matrix(QQ, [[ 1,  2,  3,  4],
            ....:                 [-2,  1, -4,  3],
            ....:                 [-3,  4,  1, -2],
            ....:                 [-4, -3,  2,  1]])
            sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
            [1 + 2*i + 3*j + 4*k]

        TESTS:

        Unembedding is the inverse of embedding::

            sage: set_random_seed()
            sage: Q = QuaternionAlgebra(QQ, -1, -1)
            sage: M = random_matrix(Q, 3)
            sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
            sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
            True

        """
        n = ZZ(M.nrows())
        if M.ncols() != n:
            raise ValueError("the matrix 'M' must be square")
        if not n.mod(4).is_zero():
            raise ValueError("the matrix 'M' must be a quaternion embedding")

        # Use the base ring of the matrix to ensure that its entries can be
        # multiplied by elements of the quaternion algebra.
        field = M.base_ring()
        Q = QuaternionAlgebra(field,-1,-1)
        i,j,k = Q.gens()

        # Go top-left to bottom-right (reading order), converting every
        # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
        # quaternion block.
        elements = []
        for l in range(n/4):
            for m in range(n/4):
                submat = ComplexMatrixEuclideanJordanAlgebra.real_unembed(
                    M[4*l:4*l+4,4*m:4*m+4] )
                if submat[0,0] != submat[1,1].conjugate():
                    raise ValueError('bad on-diagonal submatrix')
                if submat[0,1] != -submat[1,0].conjugate():
                    raise ValueError('bad off-diagonal submatrix')
                z  = submat[0,0].real()
                z += submat[0,0].imag()*i
                z += submat[0,1].real()*j
                z += submat[0,1].imag()*k
                elements.append(z)

        return matrix(Q, n/4, elements)


    @classmethod
    def trace_inner_product(cls,X,Y):
        """
        Compute a matrix inner product in this algebra directly from
        its real embedding.

        SETUP::

            sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA

        TESTS:

        This gives the same answer as the slow, default method implemented
        in :class:`MatrixEuclideanJordanAlgebra`::

            sage: set_random_seed()
            sage: J = QuaternionHermitianEJA.random_instance()
            sage: x,y = J.random_elements(2)
            sage: Xe = x.to_matrix()
            sage: Ye = y.to_matrix()
            sage: X = QuaternionHermitianEJA.real_unembed(Xe)
            sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
            sage: expected = (X*Y).trace().coefficient_tuple()[0]
            sage: actual = QuaternionHermitianEJA.trace_inner_product(Xe,Ye)
            sage: actual == expected
            True

        """
        return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/4


class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra,
                             QuaternionMatrixEuclideanJordanAlgebra):
    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)
        Euclidean Jordan algebra of dimension 6 over Real Double Field
        sage: QuaternionHermitianEJA(2, field=RR)
        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: n_max = QuaternionHermitianEJA._max_random_instance_size()
        sage: n = ZZ.random_element(1, n_max)
        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

    """
    @classmethod
    def _denormalized_basis(cls, n):
        """
        Returns a basis for the space of quaternion Hermitian n-by-n matrices.

        Why do we embed these? Basically, because all of numerical
        linear algebra assumes that you're working with vectors consisting
        of `n` entries from a field and scalars from the same field. There's
        no way to tell SageMath that (for example) the vectors contain
        complex numbers, while the scalar field is real.

        SETUP::

            sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA

        TESTS::

            sage: set_random_seed()
            sage: n = ZZ.random_element(1,5)
            sage: B = QuaternionHermitianEJA._denormalized_basis(n)
            sage: all( M.is_symmetric() for M in B )
            True

        """
        field = ZZ
        Q = QuaternionAlgebra(QQ,-1,-1)
        I,J,K = Q.gens()

        # This is like the symmetric case, but we need to be careful:
        #
        #   * We want conjugate-symmetry, not just symmetry.
        #   * The diagonal will (as a result) be real.
        #
        S = []
        for i in range(n):
            for j in range(i+1):
                Eij = matrix(Q, n, lambda k,l: k==i and l==j)
                if i == j:
                    Sij = cls.real_embed(Eij)
                    S.append(Sij)
                else:
                    # The second, third, and fourth ones have a minus
                    # because they're conjugated.
                    Sij_real = cls.real_embed(Eij + Eij.transpose())
                    S.append(Sij_real)
                    Sij_I = cls.real_embed(I*Eij - I*Eij.transpose())
                    S.append(Sij_I)
                    Sij_J = cls.real_embed(J*Eij - J*Eij.transpose())
                    S.append(Sij_J)
                    Sij_K = cls.real_embed(K*Eij - K*Eij.transpose())
                    S.append(Sij_K)

        # Since we embedded these, we can drop back to the "field" that we
        # started with instead of the quaternion algebra "Q".
        return tuple( s.change_ring(field) for s in S )


    def __init__(self, n, **kwargs):
        super(QuaternionHermitianEJA, self).__init__(self._denormalized_basis(n),
                                                     self.jordan_product,
                                                     self.trace_inner_product,
                                                     **kwargs)
        self.rank.set_cache(n)
        # TODO: cache one()!

    @staticmethod
    def _max_random_instance_size():
        r"""
        The maximum rank of a random QuaternionHermitianEJA.
        """
        return 2 # Dimension 6

    @classmethod
    def random_instance(cls, **kwargs):
        """
        Return a random instance of this type of algebra.
        """
        n = ZZ.random_element(cls._max_random_instance_size() + 1)
        return cls(n, **kwargs)


class HadamardEJA(ConcreteEuclideanJordanAlgebra):
    """
    Return the Euclidean Jordan Algebra corresponding to the set
    `R^n` under the Hadamard product.

    Note: this is nothing more than the Cartesian product of ``n``
    copies of the spin algebra. Once Cartesian product algebras
    are implemented, this can go.

    SETUP::

        sage: from mjo.eja.eja_algebra import HadamardEJA

    EXAMPLES:

    This multiplication table can be verified by hand::

        sage: J = HadamardEJA(3)
        sage: e0,e1,e2 = J.gens()
        sage: e0*e0
        e0
        sage: e0*e1
        0
        sage: e0*e2
        0
        sage: e1*e1
        e1
        sage: e1*e2
        0
        sage: e2*e2
        e2

    TESTS:

    We can change the generator prefix::

        sage: HadamardEJA(3, prefix='r').gens()
        (r0, r1, r2)

    """
    def __init__(self, n, **kwargs):
        def jordan_product(x,y):
            P = x.parent()
            return P(tuple( xi*yi for (xi,yi) in zip(x,y) ))
        def inner_product(x,y):
            return x.inner_product(y)

        # Don't orthonormalize because our basis is already
        # orthonormal with respect to our inner-product.
        if not 'orthonormalize' in kwargs:
            kwargs['orthonormalize'] = False

        # But also don't pass check_field=False here, because the user
        # can pass in a field!
        standard_basis = FreeModule(ZZ, n).basis()
        super(HadamardEJA, self).__init__(standard_basis,
                                          jordan_product,
                                          inner_product,
                                          check_axioms=False,
                                          **kwargs)
        self.rank.set_cache(n)

        if n == 0:
            self.one.set_cache( self.zero() )
        else:
            self.one.set_cache( sum(self.gens()) )

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

    @classmethod
    def random_instance(cls, **kwargs):
        """
        Return a random instance of this type of algebra.
        """
        n = ZZ.random_element(cls._max_random_instance_size() + 1)
        return cls(n, **kwargs)


class BilinearFormEJA(ConcreteEuclideanJordanAlgebra):
    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, **kwargs):
        if not B.is_positive_definite():
            raise ValueError("bilinear form is not positive-definite")

        def inner_product(x,y):
            return (B*x).inner_product(y)

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

        n = B.nrows()
        standard_basis = FreeModule(ZZ, n).basis()
        super(BilinearFormEJA, self).__init__(standard_basis,
                                              jordan_product,
                                              inner_product,
                                              **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_size():
        r"""
        The maximum dimension of a random BilinearFormEJA.
        """
        return 5

    @classmethod
    def random_instance(cls, **kwargs):
        """
        Return a random instance of this algebra.
        """
        n = ZZ.random_element(cls._max_random_instance_size() + 1)
        if n.is_zero():
            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: e0,e1,e2,e3 = J.gens()
        sage: e0*e0
        e0
        sage: e0*e1
        e1
        sage: e0*e2
        e2
        sage: e0*e3
        e3
        sage: e1*e2
        0
        sage: e1*e3
        0
        sage: e2*e3
        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, **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 not 'orthonormalize' in kwargs:
            kwargs['orthonormalize'] = False

        # But also don't pass check_field=False here, because the user
        # can pass in a field!
        super(JordanSpinEJA, self).__init__(B,
                                            check_axioms=False,
                                            **kwargs)

    @staticmethod
    def _max_random_instance_size():
        r"""
        The maximum dimension of a random JordanSpinEJA.
        """
        return 5

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

        Needed here to override the implementation for ``BilinearFormEJA``.
        """
        n = ZZ.random_element(cls._max_random_instance_size() + 1)
        return cls(n, **kwargs)


class TrivialEJA(ConcreteEuclideanJordanAlgebra):
    """
    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, **kwargs):
        jordan_product = lambda x,y: x
        inner_product = lambda x,y: 0
        basis = ()
        super(TrivialEJA, self).__init__(basis,
                                         jordan_product,
                                         inner_product,
                                         **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, **kwargs):
        # We don't take a "size" argument so the superclass method is
        # inappropriate for us.
        return cls(**kwargs)

class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
    r"""
    The external (orthogonal) direct sum of two other Euclidean Jordan
    algebras. Essentially the Cartesian product of its two factors.
    Every Euclidean Jordan algebra decomposes into an orthogonal
    direct sum of simple Euclidean Jordan algebras, so no generality
    is lost by providing only this construction.

    SETUP::

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

    EXAMPLES::

        sage: J1 = HadamardEJA(2)
        sage: J2 = RealSymmetricEJA(3)
        sage: J = DirectSumEJA(J1,J2)
        sage: J.dimension()
        8
        sage: J.rank()
        5

    TESTS:

    The external direct sum construction is only valid when the two factors
    have the same base ring; an error is raised otherwise::

        sage: set_random_seed()
        sage: J1 = random_eja(field=AA)
        sage: J2 = random_eja(field=QQ,orthonormalize=False)
        sage: J = DirectSumEJA(J1,J2)
        Traceback (most recent call last):
        ...
        ValueError: algebras must share the same base field

    """
    def __init__(self, J1, J2, **kwargs):
        if J1.base_ring() != J2.base_ring():
            raise ValueError("algebras must share the same base field")
        field = J1.base_ring()

        self._factors = (J1, J2)
        n1 = J1.dimension()
        n2 = J2.dimension()
        n = n1+n2
        V = VectorSpace(field, n)
        mult_table = [ [ V.zero() for j in range(i+1) ]
                       for i in range(n) ]
        for i in range(n1):
            for j in range(i+1):
                p = (J1.monomial(i)*J1.monomial(j)).to_vector()
                mult_table[i][j] = V(p.list() + [field.zero()]*n2)

        for i in range(n2):
            for j in range(i+1):
                p = (J2.monomial(i)*J2.monomial(j)).to_vector()
                mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())

        # TODO: build the IP table here from the two constituent IP
        # matrices (it'll be block diagonal, I think).
        ip_table = [ [ field.zero() for j in range(i+1) ]
                       for i in range(n) ]
        super(DirectSumEJA, self).__init__(field,
                                           mult_table,
                                           ip_table,
                                           check_axioms=False,
                                           **kwargs)
        self.rank.set_cache(J1.rank() + J2.rank())


    def factors(self):
        r"""
        Return the pair of this algebra's factors.

        SETUP::

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

        EXAMPLES::

            sage: J1 = HadamardEJA(2, field=QQ)
            sage: J2 = JordanSpinEJA(3, field=QQ)
            sage: J = DirectSumEJA(J1,J2)
            sage: J.factors()
            (Euclidean Jordan algebra of dimension 2 over Rational Field,
             Euclidean Jordan algebra of dimension 3 over Rational Field)

        """
        return self._factors

    def projections(self):
        r"""
        Return a pair of projections onto this algebra's factors.

        SETUP::

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

        EXAMPLES::

            sage: J1 = JordanSpinEJA(2)
            sage: J2 = ComplexHermitianEJA(2)
            sage: J = DirectSumEJA(J1,J2)
            sage: (pi_left, pi_right) = J.projections()
            sage: J.one().to_vector()
            (1, 0, 1, 0, 0, 1)
            sage: pi_left(J.one()).to_vector()
            (1, 0)
            sage: pi_right(J.one()).to_vector()
            (1, 0, 0, 1)

        """
        (J1,J2) = self.factors()
        m = J1.dimension()
        n = J2.dimension()
        V_basis = self.vector_space().basis()
        # Need to specify the dimensions explicitly so that we don't
        # wind up with a zero-by-zero matrix when we want e.g. a
        # zero-by-two matrix (important for composing things).
        P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
        P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
        pi_left = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J1,P1)
        pi_right = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J2,P2)
        return (pi_left, pi_right)

    def inclusions(self):
        r"""
        Return the pair of inclusion maps from our factors into us.

        SETUP::

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

        EXAMPLES::

            sage: J1 = JordanSpinEJA(3)
            sage: J2 = RealSymmetricEJA(2)
            sage: J = DirectSumEJA(J1,J2)
            sage: (iota_left, iota_right) = J.inclusions()
            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:

        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 = DirectSumEJA(J1,J2)
            sage: (iota_left, iota_right) = J.inclusions()
            sage: (pi_left, pi_right) = J.projections()
            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

        """
        (J1,J2) = self.factors()
        m = J1.dimension()
        n = J2.dimension()
        V_basis = self.vector_space().basis()
        # Need to specify the dimensions explicitly so that we don't
        # wind up with a zero-by-zero matrix when we want e.g. a
        # two-by-zero matrix (important for composing things).
        I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m])
        I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:])
        iota_left = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,self,I1)
        iota_right = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,self,I2)
        return (iota_left, iota_right)

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

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

        SETUP::


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

        EXAMPLE::

            sage: J1 = HadamardEJA(3,field=QQ)
            sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
            sage: J = DirectSumEJA(J1,J2)
            sage: x1 = J1.one()
            sage: x2 = x1
            sage: y1 = J2.one()
            sage: y2 = y1
            sage: x1.inner_product(x2)
            3
            sage: y1.inner_product(y2)
            2
            sage: J.one().inner_product(J.one())
            5

        """
        (pi_left, pi_right) = self.projections()
        x1 = pi_left(x)
        x2 = pi_right(x)
        y1 = pi_left(y)
        y2 = pi_right(y)

        return (x1.inner_product(y1) + x2.inner_product(y2))



random_eja = ConcreteEuclideanJordanAlgebra.random_instance