eja: update todo, and rename "natural" to "matrix".

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

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index 34fc0c3bc7b1dfe250fd5f7ea5de4fbbc52c66e0..9d53bc5fdfd8fca177e43e49a187387fbe5a6085 100644 (file)
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -3,6 +3,17 @@ 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.
 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 itertools import repeat


@@ -13,16 +24,13 @@ 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.matrix.constructor import matrix
 from sage.matrix.matrix_space import MatrixSpace
 from sage.misc.cachefunc import cached_method

-from sage.misc.lazy_import import lazy_import
-from sage.misc.prandom import choice

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

-lazy_import('mjo.eja.eja_subalgebra',
-            'FiniteDimensionalEuclideanJordanSubalgebra')
+from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator

 from mjo.eja.eja_utils import _mat2vec
 
 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 from mjo.eja.eja_utils import _mat2vec
 
 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):


@@ -46,7 +54,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: J(1)
             Traceback (most recent call last):
             ...
             sage: J(1)
             Traceback (most recent call last):
             ...

-            ValueError: not a naturally-represented algebra element
+            ValueError: not an element of this algebra

 
         """
         return None
 
         """
         return None


@@ -54,15 +62,18 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
     def __init__(self,
                  field,
                  mult_table,
     def __init__(self,
                  field,
                  mult_table,

-                 rank,

                  prefix='e',
                  category=None,
                  prefix='e',
                  category=None,

-                 natural_basis=None,
-                 check=True):
+                 matrix_basis=None,
+                 check_field=True,
+                 check_axioms=True):

         """
         SETUP::
 
         """
         SETUP::
 

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

 
         EXAMPLES:
 
 
         EXAMPLES:
 


@@ -76,23 +87,35 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         TESTS:
 
 
         TESTS:
 

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

 
             sage: JordanSpinEJA(2,QQbar)
             Traceback (most recent call last):
             ...
 
             sage: JordanSpinEJA(2,QQbar)
             Traceback (most recent call last):
             ...

-            ValueError: field is not real
+            ValueError: scalar field is not real
+
+        The multiplication table must be square with ``check_axioms=True``::
+
+            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()))
+            Traceback (most recent call last):
+            ...
+            ValueError: multiplication table is not square

 
         """
 
         """

-        if check:
+        if check_field:

             if not field.is_subring(RR):
                 # Note: this does return true for the real algebraic
             if not field.is_subring(RR):
                 # Note: this does return true for the real algebraic

-                # field, and any quadratic field where we've specified
-                # a real embedding.
-                raise ValueError('field is not real')
+                # field, the rationals, and any quadratic field where
+                # we've specified a real embedding.
+                raise ValueError("scalar field is not real")

 
 

-        self._rank = rank
-        self._natural_basis = natural_basis
+        # The multiplication table had better be square
+        n = len(mult_table)
+        if check_axioms:
+            if not all( len(l) == n for l in mult_table ):
+                raise ValueError("multiplication table is not square")
+
+        self._matrix_basis = matrix_basis

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


@@ -100,7 +123,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
         fda.__init__(field,
 
         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
         fda.__init__(field,

-                     range(len(mult_table)),
+                     range(n),

                      prefix=prefix,
                      category=category)
         self.print_options(bracket='')
                      prefix=prefix,
                      category=category)
         self.print_options(bracket='')


@@ -116,10 +139,17 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             for ls in mult_table
         ]
 
             for ls in mult_table
         ]
 

+        if check_axioms:
+            if not self._is_commutative():
+                raise ValueError("algebra is not commutative")
+            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):
         """
 
     def _element_constructor_(self, elt):
         """

-        Construct an element of this algebra from its natural
+        Construct an element of this algebra from its vector or matrix

         representation.
 
         This gets called only after the parent element _call_ method
         representation.
 
         This gets called only after the parent element _call_ method


@@ -147,13 +177,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: J(A)
             Traceback (most recent call last):
             ...
             sage: J(A)
             Traceback (most recent call last):
             ...

-            ArithmeticError: vector is not in free module
+            ValueError: not an element of this algebra

 
         TESTS:
 
         Ensure that we can convert any element of the two non-matrix
 
         TESTS:
 
         Ensure that we can convert any element of the two non-matrix

-        simple algebras (whose natural representations are their usual
-        vector representations) back and forth faithfully::
+        simple algebras (whose matrix representations are columns)
+        back and forth faithfully::

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


@@ -164,9 +194,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: x = J.random_element()
             sage: J(x.to_vector().column()) == x
             True
             sage: x = J.random_element()
             sage: J(x.to_vector().column()) == x
             True

-

         """
         """

-        msg = "not a naturally-represented algebra element"
+        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.
         if elt == 0:
             # The superclass implementation of random_element()
             # needs to be able to coerce "0" into the algebra.


@@ -178,22 +207,24 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             # that the integer 3 belongs to the space of 2-by-2 matrices.
             raise ValueError(msg)
 
             # that the integer 3 belongs to the space of 2-by-2 matrices.
             raise ValueError(msg)
 

-        natural_basis = self.natural_basis()
-        basis_space = natural_basis[0].matrix_space()
-        if elt not in basis_space:
+        if elt not in self.matrix_space():

             raise ValueError(msg)
 
         # Thanks for nothing! Matrix spaces aren't vector spaces in
             raise ValueError(msg)
 
         # Thanks for nothing! Matrix spaces aren't vector spaces in

-        # Sage, so we have to figure out its natural-basis coordinates
+        # 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.
         # ourselves. We use the basis space's ring instead of the
         # element's ring because the basis space might be an algebraic
         # closure whereas the base ring of the 3-by-3 identity matrix
         # could be QQ instead of QQbar.

-        V = VectorSpace(basis_space.base_ring(), elt.nrows()*elt.ncols())
-        W = V.span_of_basis( _mat2vec(s) for s in natural_basis )
-        coords =  W.coordinate_vector(_mat2vec(elt))
-        return self.from_vector(coords)
+        V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols())
+        W = V.span_of_basis( _mat2vec(s) for s in self.matrix_basis() )

 
 

+        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):
         """
 
     def _repr_(self):
         """


@@ -219,168 +250,74 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
     def product_on_basis(self, i, j):
         return self._multiplication_table[i][j]
 
     def product_on_basis(self, i, j):
         return self._multiplication_table[i][j]
 

-    def _a_regular_element(self):
-        """
-        Guess a regular element. Needed to compute the basis for our
-        characteristic polynomial coefficients.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import random_eja
-
-        TESTS:
-
-        Ensure that this hacky method succeeds for every algebra that we
-        know how to construct::
-
-            sage: set_random_seed()
-            sage: J = random_eja()
-            sage: J._a_regular_element().is_regular()
-            True
+    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.

         """
         """

-        gs = self.gens()
-        z = self.sum( (i+1)*gs[i] for i in range(len(gs)) )
-        if not z.is_regular():
-            raise ValueError("don't know a regular element")
-        return z
-
-
-    @cached_method
-    def _charpoly_basis_space(self):
-        """
-        Return the vector space spanned by the basis used in our
-        characteristic polynomial coefficients. This is used not only to
-        compute those coefficients, but also any time we need to
-        evaluate the coefficients (like when we compute the trace or
-        determinant).
-        """
-        z = self._a_regular_element()
-        # Don't use the parent vector space directly here in case this
-        # happens to be a subalgebra. In that case, we would be e.g.
-        # two-dimensional but span_of_basis() would expect three
-        # coordinates.
-        V = VectorSpace(self.base_ring(), self.vector_space().dimension())
-        basis = [ (z**k).to_vector() for k in range(self.rank()) ]
-        V1 = V.span_of_basis( basis )
-        b =  (V1.basis() + V1.complement().basis())
-        return V.span_of_basis(b)
-
+        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)`.

 
 

-    @cached_method
-    def _charpoly_coeff(self, i):
-        """
-        Return the coefficient polynomial "a_{i}" of this algebra's
-        general characteristic polynomial.
-
-        Having this be a separate cached method lets us compute and
-        store the trace/determinant (a_{r-1} and a_{0} respectively)
-        separate from the entire characteristic polynomial.
-        """
-        (A_of_x, x, xr, detA) = self._charpoly_matrix_system()
-        R = A_of_x.base_ring()
-
-        if i == self.rank():
-            return R.one()
-        if i > self.rank():
-            # Guaranteed by theory
-            return R.zero()
-
-        # Danger: the in-place modification is done for performance
-        # reasons (reconstructing a matrix with huge polynomial
-        # entries is slow), but I don't know how cached_method works,
-        # so it's highly possible that we're modifying some global
-        # list variable by reference, here. In other words, you
-        # probably shouldn't call this method twice on the same
-        # algebra, at the same time, in two threads
-        Ai_orig = A_of_x.column(i)
-        A_of_x.set_column(i,xr)
-        numerator = A_of_x.det()
-        A_of_x.set_column(i,Ai_orig)
-
-        # We're relying on the theory here to ensure that each a_i is
-        # indeed back in R, and the added negative signs are to make
-        # the whole charpoly expression sum to zero.
-        return R(-numerator/detA)
+        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, RDF) with check_axioms=True.
+        epsilon = 1e-16

 
 

-    @cached_method
-    def _charpoly_matrix_system(self):
-        """
-        Compute the matrix whose entries A_ij are polynomials in
-        X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
-        corresponding to `x^r` and the determinent of the matrix A =
-        [A_ij]. In other words, all of the fixed (cachable) data needed
-        to compute the coefficients of the characteristic polynomial.
-        """
-        r = self.rank()
-        n = self.dimension()
+        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)

 
 

-        # Turn my vector space into a module so that "vectors" can
-        # have multivatiate polynomial entries.
-        names = tuple('X' + str(i) for i in range(1,n+1))
-        R = PolynomialRing(self.base_ring(), names)
-
-        # Using change_ring() on the parent's vector space doesn't work
-        # here because, in a subalgebra, that vector space has a basis
-        # and change_ring() tries to bring the basis along with it. And
-        # that doesn't work unless the new ring is a PID, which it usually
-        # won't be.
-        V = FreeModule(R,n)
-
-        # Now let x = (X1,X2,...,Xn) be the vector whose entries are
-        # indeterminates...
-        x = V(names)
-
-        # And figure out the "left multiplication by x" matrix in
-        # that setting.
-        lmbx_cols = []
-        monomial_matrices = [ self.monomial(i).operator().matrix()
-                              for i in range(n) ] # don't recompute these!
-        for k in range(n):
-            ek = self.monomial(k).to_vector()
-            lmbx_cols.append(
-              sum( x[i]*(monomial_matrices[i]*ek)
-                   for i in range(n) ) )
-        Lx = matrix.column(R, lmbx_cols)
-
-        # Now we can compute powers of x "symbolically"
-        x_powers = [self.one().to_vector(), x]
-        for d in range(2, r+1):
-            x_powers.append( Lx*(x_powers[-1]) )
-
-        idmat = matrix.identity(R, n)
-
-        W = self._charpoly_basis_space()
-        W = W.change_ring(R.fraction_field())
-
-        # Starting with the standard coordinates x = (X1,X2,...,Xn)
-        # and then converting the entries to W-coordinates allows us
-        # to pass in the standard coordinates to the charpoly and get
-        # back the right answer. Specifically, with x = (X1,X2,...,Xn),
-        # we have
-        #
-        #   W.coordinates(x^2) eval'd at (standard z-coords)
-        #     =
-        #   W-coords of (z^2)
-        #     =
-        #   W-coords of (standard coords of x^2 eval'd at std-coords of z)
-        #
-        # We want the middle equivalent thing in our matrix, but use
-        # the first equivalent thing instead so that we can pass in
-        # standard coordinates.
-        x_powers = [ W.coordinate_vector(xp) for xp in x_powers ]
-        l2 = [idmat.column(k-1) for k in range(r+1, n+1)]
-        A_of_x = matrix.column(R, n, (x_powers[:r] + l2))
-        return (A_of_x, x, x_powers[r], A_of_x.det())
+                    if self.base_ring().is_exact():
+                        if diff != 0:
+                            return False
+                    else:
+                        if diff.abs() > epsilon:
+                            return False

 
 

+        return True

 
     @cached_method
 
     @cached_method

-    def characteristic_polynomial(self):
+    def characteristic_polynomial_of(self):

         """
         """

-        Return a characteristic polynomial that works for all elements
-        of this algebra.
+        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 resulting polynomial has `n+1` variables, where `n` is the
         dimension of this algebra. The first `n` variables correspond to


@@ -400,7 +337,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         Alizadeh, Example 11.11::
 
             sage: J = JordanSpinEJA(3)
         Alizadeh, Example 11.11::
 
             sage: J = JordanSpinEJA(3)

-            sage: p = J.characteristic_polynomial(); p
+            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)
             X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
             sage: xvec = J.one().to_vector()
             sage: p(*xvec)


@@ -413,28 +350,51 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         any argument::
 
             sage: J = TrivialEJA()
         any argument::
 
             sage: J = TrivialEJA()

-            sage: J.characteristic_polynomial()
+            sage: J.characteristic_polynomial_of()

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

-        # The list of coefficient polynomials a_0, a_1, a_2, ..., a_n.
-        a = [ self._charpoly_coeff(i) for i in range(r+1) ]
+        # 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.
 
         # 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.

-        R = a[0].parent()

         S = PolynomialRing(self.base_ring(),'t')
         t = S.gen(0)
         S = PolynomialRing(self.base_ring(),'t')
         t = S.gen(0)

-        S = PolynomialRing(S, R.variable_names())
-        t = S(t)
+        if r > 0:
+            R = a[0].parent()
+            S = PolynomialRing(S, R.variable_names())
+            t = S(t)

 
 

-        return sum( a[k]*(t**k) for k in range(len(a)) )
+        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,QQ)
+            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):
         """
 
     def inner_product(self, x, y):
         """


@@ -446,7 +406,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         SETUP::
 
 
         SETUP::
 

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

 
         EXAMPLES:
 
 
         EXAMPLES:
 


@@ -459,10 +421,34 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: (x*y).inner_product(z) == y.inner_product(x*z)
             True
 
             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

         """
         """

-        X = x.natural_representation()
-        Y = y.natural_representation()
-        return self.natural_inner_product(X,Y)
+        B = self._inner_product_matrix
+        return (B*x.to_vector()).inner_product(y.to_vector())

 
 
     def is_trivial(self):
 
 
     def is_trivial(self):


@@ -526,18 +512,33 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return table(M, header_row=True, header_column=True, frame=True)
 
 
         return table(M, header_row=True, header_column=True, frame=True)
 
 

-    def natural_basis(self):
+    def matrix_basis(self):

         """
         """

-        Return a more-natural representation of this algebra's basis.
+        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.

 
 

-        Every finite-dimensional Euclidean Jordan Algebra is a direct
-        sum of five simple algebras, four of which comprise Hermitian
-        matrices. This method returns the original "natural" basis
-        for our underlying vector space. (Typically, the natural basis
-        is used to construct the multiplication table in the first place.)
+        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.

 
 

-        Note that this will always return a matrix. The standard basis
-        in `R^n` will be returned as `n`-by-`1` column matrices.
+        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::
 
 
         SETUP::
 


@@ -549,7 +550,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: J = RealSymmetricEJA(2)
             sage: J.basis()
             Finite family {0: e0, 1: e1, 2: e2}
             sage: J = RealSymmetricEJA(2)
             sage: J.basis()
             Finite family {0: e0, 1: e1, 2: e2}

-            sage: J.natural_basis()
+            sage: J.matrix_basis()

             (
             [1 0]  [                  0 0.7071067811865475?]  [0 0]
             [0 0], [0.7071067811865475?                   0], [0 1]
             (
             [1 0]  [                  0 0.7071067811865475?]  [0 0]
             [0 0], [0.7071067811865475?                   0], [0 1]


@@ -560,43 +561,38 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: J = JordanSpinEJA(2)
             sage: J.basis()
             Finite family {0: e0, 1: e1}
             sage: J = JordanSpinEJA(2)
             sage: J.basis()
             Finite family {0: e0, 1: e1}

-            sage: J.natural_basis()
+            sage: J.matrix_basis()

             (
             [1]  [0]
             [0], [1]
             )
             (
             [1]  [0]
             [0], [1]
             )

-

         """
         """

-        if self._natural_basis is None:
-            M = self.natural_basis_space()
+        if self._matrix_basis is None:
+            M = self.matrix_space()

             return tuple( M(b.to_vector()) for b in self.basis() )
         else:
             return tuple( M(b.to_vector()) for b in self.basis() )
         else:

-            return self._natural_basis
+            return self._matrix_basis

 
 
 
 

-    def natural_basis_space(self):
-        """
-        Return the matrix space in which this algebra's natural basis
-        elements live.
+    def matrix_space(self):

         """
         """

-        if self._natural_basis is None or len(self._natural_basis) == 0:
-            return MatrixSpace(self.base_ring(), self.dimension(), 1)
-        else:
-            return self._natural_basis[0].matrix_space()
-
+        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).

 
 

-    @staticmethod
-    def natural_inner_product(X,Y):
-        """
-        Compute the inner product of two naturally-represented elements.
+        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.

 
 

-        For example in the real symmetric matrix EJA, this will compute
-        the trace inner-product of two n-by-n symmetric matrices. The
-        default should work for the real cartesian product EJA, the
-        Jordan spin EJA, and the real symmetric matrices. The others
-        will have to be overridden.
+        Matrix algebras override this with something more useful.

         """
         """

-        return (X.conjugate_transpose()*Y).trace()
+        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
 
 
     @cached_method


@@ -634,6 +630,16 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: actual == expected
             True
 
             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.
         """
         # We can brute-force compute the matrices of the operators
         # that correspond to the basis elements of this algebra.


@@ -645,19 +651,20 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         # appeal to the "long vectors" isometry.
         oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]
 
         # appeal to the "long vectors" isometry.
         oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]
 

-        # Now we use basis linear algebra to find the coefficients,
+        # 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.
         # of the matrices-as-vectors-linear-combination, which should
         # work for the original algebra basis too.

-        A = matrix.column(self.base_ring(), oper_vecs)
+        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()) )
 
 
         # 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.
-        coeffs = A.solve_right(b)
-        return self.linear_combination(zip(self.gens(), coeffs))
+        # 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):
 
 
     def peirce_decomposition(self, c):


@@ -712,6 +719,25 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             Vector space of degree 6 and dimension 2...
             sage: J1
             Euclidean Jordan algebra of dimension 3...
             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:
 
 
         TESTS:
 


@@ -726,9 +752,10 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
             True
 
             sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
             True
 

-        The identity elements in the two subalgebras are the
-        projections onto their respective subspaces of the
-        superalgebra's identity element::
+        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: set_random_seed()
             sage: J = random_eja()


@@ -738,6 +765,16 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             ....:         x = J.random_element()
             sage: c = x.subalgebra_idempotent()
             sage: J0,J5,J1 = J.peirce_decomposition(c)
             ....:         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()
             sage: J1(c) == J1.one()
             True
             sage: J0(J.one() - c) == J0.one()


@@ -747,6 +784,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         if not c.is_idempotent():
             raise ValueError("element is not idempotent: %s" % c)
 
         if not c.is_idempotent():
             raise ValueError("element is not idempotent: %s" % c)
 

+        from mjo.eja.eja_subalgebra import FiniteDimensionalEuclideanJordanSubalgebra
+

         # Default these to what they should be if they turn out to be
         # trivial, because eigenspaces_left() won't return eigenvalues
         # corresponding to trivial spaces (e.g. it returns only the
         # 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


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

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

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


@@ -773,143 +814,89 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return (J0, J5, J1)
 
 
         return (J0, J5, J1)
 
 

-    def random_elements(self, count):
-        """
-        Return ``count`` random elements as a tuple.
+    def random_element(self, thorough=False):
+        r"""
+        Return a random element of this algebra.

 
 

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

 
 

-            sage: from mjo.eja.eja_algebra import JordanSpinEJA
+        INPUT:

 
 

-        EXAMPLES::
+        - ``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()

 
 

-            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 self.from_vector(V.coordinate_vector(v))

 
 

+    def random_elements(self, count, thorough=False):

         """
         """

-        return  tuple( self.random_element() for idx in range(count) )
+        Return ``count`` random elements as a tuple.

 
 

+        INPUT:

 
 

-    def _rank_computation1(self):
-        r"""
-        Compute the rank of this algebra using highly suspicious voodoo.
-
-        ALGORITHM:
-
-        We first compute the basis representation of the operator L_x
-        using polynomial indeterminates are placeholders for the
-        coordinates of "x", which is arbitrary. We then use that
-        matrix to compute the (polynomial) entries of x^0, x^1, ...,
-        x^d,... for increasing values of "d", starting at zero. The
-        idea is that. If we also add "coefficient variables" a_0,
-        a_1,...  to the ring, we can form the linear combination
-        a_0*x^0 + ... + a_d*x^d = 0, and ask what dimension the
-        solution space has as an affine variety. When "d" is smaller
-        than the rank, we expect that dimension to be the number of
-        coordinates of "x", since we can set *those* to whatever we
-        want, but linear independence forces the coefficients a_i to
-        be zero. Eventually, when "d" passes the rank, the dimension
-        of the solution space begins to grow, because we can *still*
-        set the coordinates of "x" arbitrarily, but now there are some
-        coefficients that make the sum zero as well. So, when the
-        dimension of the variety jumps, we return the corresponding
-        "d" as the rank of the algebra. This appears to work.
+        - ``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::
 
 
         SETUP::
 

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

 
         EXAMPLES::
 
 
         EXAMPLES::
 

-            sage: J = HadamardEJA(5)
-            sage: J._rank_computation() == J.rank()
-            True
-            sage: J = JordanSpinEJA(5)
-            sage: J._rank_computation() == J.rank()
-            True
-            sage: J = RealSymmetricEJA(4)
-            sage: J._rank_computation() == J.rank()
-            True
-            sage: J = ComplexHermitianEJA(3)
-            sage: J._rank_computation() == J.rank()
+            sage: J = JordanSpinEJA(3)
+            sage: x,y,z = J.random_elements(3)
+            sage: all( [ x in J, y in J, z in J ])

             True
             True

-            sage: J = QuaternionHermitianEJA(2)
-            sage: J._rank_computation() == J.rank()
+            sage: len( J.random_elements(10) ) == 10

             True
 
         """
             True
 
         """

-        n = self.dimension()
-        var_names = [ "X" + str(z) for z in range(1,n+1) ]
-        d = 0
-        ideal_dim = len(var_names)
-        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[d+k]*self.monomial(k).operator().matrix()[i,j]
-                        for k in range(n) )
+        return tuple( self.random_element(thorough)
+                      for idx in range(count) )

 
 

-        while ideal_dim == len(var_names):
-            coeff_names = [ "a" + str(z) for z in range(d) ]
-            R = PolynomialRing(self.base_ring(), coeff_names + var_names)
-            vars = R.gens()
-            L_x = matrix(R, n, n, L_x_i_j)
-            x_powers = [ vars[k]*(L_x**k)*self.one().to_vector()
-                         for k in range(d) ]
-            eqs = [ sum(x_powers[k][j] for k in range(d)) for j in range(n) ]
-            ideal_dim = R.ideal(eqs).dimension()
-            d += 1
-
-        # Subtract one because we increment one too many times, and
-        # subtract another one because "d" is one greater than the
-        # answer anyway; when d=3, we go up to x^2.
-        return d-2
-
-    def _rank_computation2(self):
-        r"""
-        Instead of using the dimension of an ideal, find the rank of a
-        matrix containing polynomials.
-        """
-        n = self.dimension()
-        var_names = [ "X" + str(z) for z in range(1,n+1) ]
-        R = PolynomialRing(self.base_ring(), var_names)
-        vars = R.gens()

 
 

-        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(R, n, n, L_x_i_j)
-        x_powers = [ (vars[k]*(L_x**k)*self.one().to_vector()).row()
-                     for k in range(n) ]
-
-        from sage.matrix.constructor import block_matrix
-        M = block_matrix(n,1,x_powers)
-        return M.rank()
-
-    def _rank_computation3(self):
+    @cached_method
+    def _charpoly_coefficients(self):

         r"""
         r"""

-        Similar to :meth:`_rank_computation2`, but it stops echelonizing
-        as soon as it hits the first zero row.
+        The `r` polynomial coefficients of the "characteristic polynomial
+        of" function.

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

-        if n == 0:
-            return 0
-        elif n == 1:
-            return 1
-
-        var_names = [ "X" + str(z) for z in range(1,n+1) ]
-        R = PolynomialRing(self.base_ring(), var_names)
+        R = self.coordinate_polynomial_ring()

         vars = R.gens()
         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
 
         def L_x_i_j(i,j):
             # From a result in my book, these are the entries of the


@@ -917,39 +904,48 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
                         for k in range(n) )
 
             return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
                         for k in range(n) )
 

-        L_x = matrix(R, n, n, L_x_i_j)
-        x_powers = [ vars[k]*(L_x**k)*self.one().to_vector()
-                     for k in range(n) ]
-
-        # Can assume n >= 2
-        rows = [x_powers[0]]
-        M = matrix(rows)
-        old_rank = 1
-
-        for d in range(1,n):
-            rows = M.rows() + [x_powers[d]]
-            M = matrix(rows)
-            M.echelonize()
-            new_rank = M.rank()
-            if new_rank == old_rank:
-                return new_rank
-            else:
-                old_rank = new_rank
+        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):
     def rank(self):

-        """
+        r"""

         Return the rank of this EJA.
 
         Return the rank of this EJA.
 

-        ALGORITHM:
-
-        The author knows of no algorithm to compute the rank of an EJA
-        where only the multiplication table is known. In lieu of one, we
-        require the rank to be specified when the algebra is created,
-        and simply pass along that number here.
+        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::
 
 
         SETUP::
 

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

             ....:                                  RealSymmetricEJA,
             ....:                                  ComplexHermitianEJA,
             ....:                                  QuaternionHermitianEJA,
             ....:                                  RealSymmetricEJA,
             ....:                                  ComplexHermitianEJA,
             ....:                                  QuaternionHermitianEJA,


@@ -990,8 +986,43 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: r > 0 or (r == 0 and J.is_trivial())
             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: J = HadamardEJA(4)
+            sage: J.rank.clear_cache()
+            sage: J.rank()
+            4
+
+        ::
+
+            sage: J = JordanSpinEJA(4)
+            sage: J.rank.clear_cache()
+            sage: J.rank()
+            2
+
+        ::
+
+            sage: J = RealSymmetricEJA(3)
+            sage: J.rank.clear_cache()
+            sage: J.rank()
+            3
+
+        ::
+
+            sage: J = ComplexHermitianEJA(2)
+            sage: J.rank.clear_cache()
+            sage: J.rank()
+            2
+
+        ::
+
+            sage: J = QuaternionHermitianEJA(2)
+            sage: J.rank.clear_cache()
+            sage: J.rank()
+            2

         """
         """

-        return self._rank
+        return len(self._charpoly_coefficients())

 
 
     def vector_space(self):
 
 
     def vector_space(self):


@@ -1015,179 +1046,150 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
     Element = FiniteDimensionalEuclideanJordanAlgebraElement
 
 
     Element = FiniteDimensionalEuclideanJordanAlgebraElement
 
 

-class KnownRankEJA(object):
-    """
-    A class for algebras that we actually know we can construct.  The
-    main issue is that, for most of our methods to make sense, we need
-    to know the rank of our algebra. Thus we can't simply generate a
-    "random" algebra, or even check that a given basis and product
-    satisfy the axioms; because even if everything looks OK, we wouldn't
-    know the rank we need to actuallty build the thing.
-
-    Not really a subclass of FDEJA because doing that causes method
-    resolution errors, e.g.
-
-      TypeError: Error when calling the metaclass bases
-      Cannot create a consistent method resolution
-      order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra,
-      KnownRankEJA
-
+class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
+    r"""
+    Algebras whose basis consists of vectors with rational
+    entries. Equivalently, algebras whose multiplication tables
+    contain only rational coefficients.
+
+    When an EJA has a basis that can be made rational, we can speed up
+    the computation of its characteristic polynomial by doing it over
+    ``QQ``. All of the named EJA constructors that we provide fall
+    into this category.

     """
     """

-    @staticmethod
-    def _max_test_case_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 quite vague -- 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.
-
-        We default to five in this class, which is safe in `R^n`. The
-        matrix algebra subclasses (or any class where the "size" is
-        interpreted to be far less than the dimension) should override
-        with a smaller number.
-        """
-        return 5
+    @cached_method
+    def _charpoly_coefficients(self):
+        r"""
+        Override the parent method with something that tries to compute
+        over a faster (non-extension) field.

 
 

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

 
 

-        Beware, this will crash for "most instances" because the
-        constructor below looks wrong.
-        """
-        if cls is TrivialEJA:
-            # The TrivialEJA class doesn't take an "n" argument because
-            # there's only one.
-            return cls(field)
+            sage: from mjo.eja.eja_algebra import JordanSpinEJA

 
 

-        n = ZZ.random_element(cls._max_test_case_size()) + 1
-        return cls(n, field, **kwargs)
+        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)::

 
 

-class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
-    """
-    Return the Euclidean Jordan Algebra corresponding to the set
-    `R^n` under the Hadamard product.
+            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:
+            # There's no need to construct *another* algebra over the
+            # rationals if this one is already over the rationals.
+            superclass = super(RationalBasisEuclideanJordanAlgebra, self)
+            return superclass._charpoly_coefficients()
+
+        mult_table = tuple(
+            map(lambda x: x.to_vector(), ls)
+            for ls in self._multiplication_table
+        )
+
+        # Do the computation over the rationals. The answer will be
+        # the same, because our basis coordinates are (essentially)
+        # rational.
+        J = FiniteDimensionalEuclideanJordanAlgebra(QQ,
+                                                    mult_table,
+                                                    check_field=False,
+                                                    check_axioms=False)
+        a = J._charpoly_coefficients()
+        return tuple(map(lambda x: x.change_ring(self.base_ring()), a))
+
+
+class ConcreteEuclideanJordanAlgebra:
+    r"""
+    A class for the Euclidean Jordan algebras that we know by name.

 
 

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

 
     TESTS:
 
 
     TESTS:
 

-    We can change the generator prefix::
+    Our basis is normalized with respect to the algebra's inner
+    product, unless we specify otherwise::

 
 

-        sage: HadamardEJA(3, prefix='r').gens()
-        (r0, r1, r2)
+        sage: set_random_seed()
+        sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
+        sage: all( b.norm() == 1 for b in J.gens() )
+        True

 
 

-    """
-    def __init__(self, n, field=AA, **kwargs):
-        V = VectorSpace(field, n)
-        mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
-                       for i in range(n) ]
+    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::

 
 

-        fdeja = super(HadamardEJA, self)
-        return fdeja.__init__(field, mult_table, rank=n, **kwargs)
+        sage: set_random_seed()
+        sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
+        sage: x = J.random_element()
+        sage: x.operator().is_self_adjoint()
+        True
+    """

 
 

-    def inner_product(self, x, y):
+    @staticmethod
+    def _max_random_instance_size():

         """
         """

-        Faster to reimplement than to use natural representations.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import HadamardEJA
-
-        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: X = x.natural_representation()
-            sage: Y = y.natural_representation()
-            sage: x.inner_product(y) == J.natural_inner_product(X,Y)
-            True
+        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.

         """
         """

-        return x.to_vector().inner_product(y.to_vector())
-
-
-def random_eja(field=AA, nontrivial=False):
-    """
-    Return a "random" finite-dimensional Euclidean Jordan Algebra.
-
-    SETUP::
-
-        sage: from mjo.eja.eja_algebra import random_eja
-
-    TESTS::
-
-        sage: random_eja()
-        Euclidean Jordan algebra of dimension...
-
-    """
-    eja_classes = KnownRankEJA.__subclasses__()
-    if nontrivial:
-        eja_classes.remove(TrivialEJA)
-    classname = choice(eja_classes)
-    return classname.random_instance(field=field)
-
-
+        raise NotImplementedError

 
 

+    @classmethod
+    def random_instance(cls, field=AA, **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__())
+        return eja_class.random_instance(field)

 
 
 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
 
 
 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):

-    @staticmethod
-    def _max_test_case_size():
-        # Play it safe, since this will be squared and the underlying
-        # field can have dimension 4 (quaternions) too.
-        return 2

 
 

-    def __init__(self, field, basis, rank, normalize_basis=True, **kwargs):
+    def __init__(self, field, basis, normalize_basis=True, **kwargs):

         """
         Compared to the superclass constructor, we take a basis instead of
         a multiplication table because the latter can be computed in terms
         of the former when the product is known (like it is here).
         """
         """
         Compared to the superclass constructor, we take a basis instead of
         a multiplication table because the latter can be computed in terms
         of the former when the product is known (like it is here).
         """

-        # Used in this class's fast _charpoly_coeff() override.
+        # Used in this class's fast _charpoly_coefficients() override.

         self._basis_normalizers = None
 
         # We're going to loop through this a few times, so now's a good
         # time to ensure that it isn't a generator expression.
         basis = tuple(basis)
 
         self._basis_normalizers = None
 
         # We're going to loop through this a few times, so now's a good
         # time to ensure that it isn't a generator expression.
         basis = tuple(basis)
 

-        if rank > 1 and normalize_basis:
+        algebra_dim = len(basis)
+        degree = 0 # size of the matrices
+        if algebra_dim > 0:
+            degree = basis[0].nrows()
+
+        if algebra_dim > 1 and normalize_basis:

             # We'll need sqrt(2) to normalize the basis, and this
             # winds up in the multiplication table, so the whole
             # algebra needs to be over the field extension.
             # We'll need sqrt(2) to normalize the basis, and this
             # winds up in the multiplication table, so the whole
             # algebra needs to be over the field extension.


@@ -1198,81 +1200,96 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
                 field = field.extension(p, 'sqrt2', embedding=RLF(2).sqrt())
                 basis = tuple( s.change_ring(field) for s in basis )
             self._basis_normalizers = tuple(
                 field = field.extension(p, 'sqrt2', embedding=RLF(2).sqrt())
                 basis = tuple( s.change_ring(field) for s in basis )
             self._basis_normalizers = tuple(

-                ~(self.natural_inner_product(s,s).sqrt()) for s in basis )
+                ~(self.matrix_inner_product(s,s).sqrt()) for s in basis )

             basis = tuple(s*c for (s,c) in zip(basis,self._basis_normalizers))
 
             basis = tuple(s*c for (s,c) in zip(basis,self._basis_normalizers))
 

-        Qs = self.multiplication_table_from_matrix_basis(basis)
-
-        fdeja = super(MatrixEuclideanJordanAlgebra, self)
-        return fdeja.__init__(field,
-                              Qs,
-                              rank=rank,
-                              natural_basis=basis,
-                              **kwargs)
-
-
-    @cached_method
-    def _charpoly_coeff(self, i):
-        """
-        Override the parent method with something that tries to compute
-        over a faster (non-extension) field.
-        """
-        if self._basis_normalizers is None:
-            # We didn't normalize, so assume that the basis we started
-            # with had entries in a nice field.
-            return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coeff(i)
-        else:
-            basis = ( (b/n) for (b,n) in zip(self.natural_basis(),
-                                             self._basis_normalizers) )
-
-            # Do this over the rationals and convert back at the end.
-            J = MatrixEuclideanJordanAlgebra(QQ,
-                                             basis,
-                                             self.rank(),
-                                             normalize_basis=False)
-            (_,x,_,_) = J._charpoly_matrix_system()
-            p = J._charpoly_coeff(i)
-            # p might be missing some vars, have to substitute "optionally"
-            pairs = zip(x.base_ring().gens(), self._basis_normalizers)
-            substitutions = { v: v*c for (v,c) in pairs }
-            result = p.subs(substitutions)
-
-            # The result of "subs" can be either a coefficient-ring
-            # element or a polynomial. Gotta handle both cases.
-            if result in QQ:
-                return self.base_ring()(result)
-            else:
-                return result.change_ring(self.base_ring())
-
-
-    @staticmethod
-    def multiplication_table_from_matrix_basis(basis):
-        """
-        At least three of the five simple Euclidean Jordan algebras have the
-        symmetric multiplication (A,B) |-> (AB + BA)/2, where the
-        multiplication on the right is matrix multiplication. Given a basis
-        for the underlying matrix space, this function returns a
-        multiplication table (obtained by looping through the basis
-        elements) for an algebra of those matrices.
-        """
-        # In S^2, for example, we nominally have four coordinates even
-        # though the space is of dimension three only. The vector space V
-        # is supposed to hold the entire long vector, and the subspace W
-        # of V will be spanned by the vectors that arise from symmetric
-        # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
-        field = basis[0].base_ring()
-        dimension = basis[0].nrows()
-
-        V = VectorSpace(field, dimension**2)
+        # Now compute the multiplication and inner product tables.
+        # We have to do this *after* normalizing the basis, because
+        # scaling affects the answers.
+        V = VectorSpace(field, degree**2)

         W = V.span_of_basis( _mat2vec(s) for s in basis )
         W = V.span_of_basis( _mat2vec(s) for s in basis )

-        n = len(basis)
-        mult_table = [[W.zero() for j in range(n)] for i in range(n)]
-        for i in range(n):
-            for j in range(n):
+        mult_table = [[W.zero() for j in range(algebra_dim)]
+                                for i in range(algebra_dim)]
+        ip_table = [[W.zero() for j in range(algebra_dim)]
+                              for i in range(algebra_dim)]
+        for i in range(algebra_dim):
+            for j in range(algebra_dim):

                 mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
                 mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
 
                 mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
                 mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
 

-        return mult_table
+                try:
+                    # HACK: ignore the error here if we don't need the
+                    # inner product (as is the case when we construct
+                    # a dummy QQ-algebra for fast charpoly coefficients.
+                    ip_table[i][j] = self.matrix_inner_product(basis[i],
+                                                                basis[j])
+                except:
+                    pass
+
+        try:
+            # HACK PART DEUX
+            self._inner_product_matrix = matrix(field,ip_table)
+        except:
+            pass
+
+        super(MatrixEuclideanJordanAlgebra, self).__init__(field,
+                                                           mult_table,
+                                                           matrix_basis=basis,
+                                                           **kwargs)
+
+        if algebra_dim == 0:
+            self.one.set_cache(self.zero())
+        else:
+            n = basis[0].nrows()
+            # The identity wrt (A,B) -> (AB + BA)/2 is independent of the
+            # details of this algebra.
+            self.one.set_cache(self(matrix.identity(field,n)))
+
+
+    @cached_method
+    def _charpoly_coefficients(self):
+        r"""
+        Override the parent method with something that tries to compute
+        over a faster (non-extension) field.
+        """
+        if self._basis_normalizers is None or self.base_ring() is QQ:
+            # We didn't normalize, or the basis we started with had
+            # entries in a nice field already. Just compute the thing.
+            return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coefficients()
+
+        basis = ( (b/n) for (b,n) in zip(self.matrix_basis(),
+                                         self._basis_normalizers) )
+
+        # Do this over the rationals and convert back at the end.
+        # Only works because we know the entries of the basis are
+        # integers. The argument ``check_axioms=False`` is required
+        # because the trace inner-product method for this
+        # class is a stub and can't actually be checked.
+        J = MatrixEuclideanJordanAlgebra(QQ,
+                                         basis,
+                                         normalize_basis=False,
+                                         check_field=False,
+                                         check_axioms=False)
+        a = J._charpoly_coefficients()
+
+        # Unfortunately, changing the basis does change the
+        # coefficients of the characteristic polynomial, but since
+        # these are really the coefficients of the "characteristic
+        # polynomial of" function, everything is still nice and
+        # unevaluated. It's therefore "obvious" how scaling the
+        # basis affects the coordinate variables X1, X2, et
+        # cetera. Scaling the first basis vector up by "n" adds a
+        # factor of 1/n into every "X1" term, for example. So here
+        # we simply undo the basis_normalizer scaling that we
+        # performed earlier.
+        #
+        # The a[0] access here is safe because trivial algebras
+        # won't have any basis normalizers and therefore won't
+        # make it to this "else" branch.
+        XS = a[0].parent().gens()
+        subs_dict = { XS[i]: self._basis_normalizers[i]*XS[i]
+                      for i in range(len(XS)) }
+        return tuple( a_i.subs(subs_dict) for a_i in a )

 
 
     @staticmethod
 
 
     @staticmethod


@@ -1299,23 +1316,17 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         """
         raise NotImplementedError
 
         """
         raise NotImplementedError
 

-

     @classmethod
     @classmethod

-    def natural_inner_product(cls,X,Y):
+    def matrix_inner_product(cls,X,Y):

         Xu = cls.real_unembed(X)
         Yu = cls.real_unembed(Y)
         tr = (Xu*Yu).trace()
 
         Xu = cls.real_unembed(X)
         Yu = cls.real_unembed(Y)
         tr = (Xu*Yu).trace()
 

-        if tr in RLF:
-            # It's real already.
-            return tr
-
-        # Otherwise, try the thing that works for complex numbers; and
-        # if that doesn't work, the thing that works for quaternions.

         try:
         try:

-            return tr.vector()[0] # real part, imag part is index 1
+            # Works in QQ, AA, RDF, et cetera.
+            return tr.real()

         except AttributeError:
         except AttributeError:

-            # A quaternions doesn't have a vector() method, but does
+            # 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]
             # have coefficient_tuple() method that returns the
             # coefficients of 1, i, j, and k -- in that order.
             return tr.coefficient_tuple()[0]


@@ -1339,7 +1350,8 @@ class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         return M
 
 
         return M
 
 

-class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
+class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra,
+                       ConcreteEuclideanJordanAlgebra):

     """
     The rank-n simple EJA consisting of real symmetric n-by-n
     matrices, the usual symmetric Jordan product, and the trace inner
     """
     The rank-n simple EJA consisting of real symmetric n-by-n
     matrices, the usual symmetric Jordan product, and the trace inner


@@ -1373,7 +1385,7 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
     The dimension of this algebra is `(n^2 + n) / 2`::
 
         sage: set_random_seed()
     The dimension of this algebra is `(n^2 + n) / 2`::
 
         sage: set_random_seed()

-        sage: n_max = RealSymmetricEJA._max_test_case_size()
+        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
         sage: n = ZZ.random_element(1, n_max)
         sage: J = RealSymmetricEJA(n)
         sage: J.dimension() == (n^2 + n)/2


@@ -1384,9 +1396,9 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
         sage: set_random_seed()
         sage: J = RealSymmetricEJA.random_instance()
         sage: x,y = J.random_elements(2)
         sage: set_random_seed()
         sage: J = RealSymmetricEJA.random_instance()
         sage: x,y = J.random_elements(2)

-        sage: actual = (x*y).natural_representation()
-        sage: X = x.natural_representation()
-        sage: Y = y.natural_representation()
+        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: expected = (X*Y + Y*X)/2
         sage: actual == expected
         True


@@ -1398,24 +1410,10 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
         sage: RealSymmetricEJA(3, prefix='q').gens()
         (q0, q1, q2, q3, q4, q5)
 
         sage: RealSymmetricEJA(3, prefix='q').gens()
         (q0, q1, q2, q3, q4, q5)
 

-    Our natural basis is normalized with respect to the natural inner
-    product unless we specify otherwise::
-
-        sage: set_random_seed()
-        sage: J = RealSymmetricEJA.random_instance()
-        sage: all( b.norm() == 1 for b in J.gens() )
-        True
-
-    Since our natural basis is normalized with respect to the natural
-    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::
+    We can construct the (trivial) algebra of rank zero::

 
 

-        sage: set_random_seed()
-        sage: x = RealSymmetricEJA.random_instance().random_element()
-        sage: x.operator().matrix().is_symmetric()
-        True
+        sage: RealSymmetricEJA(0)
+        Euclidean Jordan algebra of dimension 0 over Algebraic Real Field

 
     """
     @classmethod
 
     """
     @classmethod


@@ -1451,13 +1449,24 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
 
 
     @staticmethod
 
 
     @staticmethod

-    def _max_test_case_size():
+    def _max_random_instance_size():

         return 4 # Dimension 10
 
         return 4 # Dimension 10
 

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

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

-        super(RealSymmetricEJA, self).__init__(field, basis, n, **kwargs)
+        super(RealSymmetricEJA, self).__init__(field,
+                                               basis,
+                                               check_axioms=False,
+                                               **kwargs)
+        self.rank.set_cache(n)

 
 
 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
 
 
 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):


@@ -1493,8 +1502,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         Embedding is a homomorphism (isomorphism, in fact)::
 
             sage: set_random_seed()
         Embedding is a homomorphism (isomorphism, in fact)::
 
             sage: set_random_seed()

-            sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
-            sage: n = ZZ.random_element(n_max)
+            sage: n = ZZ.random_element(3)

             sage: F = QuadraticField(-1, 'I')
             sage: X = random_matrix(F, n)
             sage: Y = random_matrix(F, n)
             sage: F = QuadraticField(-1, 'I')
             sage: X = random_matrix(F, n)
             sage: Y = random_matrix(F, n)


@@ -1589,9 +1597,9 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
 
 
     @classmethod
 
 
     @classmethod

-    def natural_inner_product(cls,X,Y):
+    def matrix_inner_product(cls,X,Y):

         """
         """

-        Compute a natural inner product in this algebra directly from
+        Compute a matrix inner product in this algebra directly from

         its real embedding.
 
         SETUP::
         its real embedding.
 
         SETUP::


@@ -1606,20 +1614,21 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             sage: set_random_seed()
             sage: J = ComplexHermitianEJA.random_instance()
             sage: x,y = J.random_elements(2)
             sage: set_random_seed()
             sage: J = ComplexHermitianEJA.random_instance()
             sage: x,y = J.random_elements(2)

-            sage: Xe = x.natural_representation()
-            sage: Ye = y.natural_representation()
+            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: X = ComplexHermitianEJA.real_unembed(Xe)
             sage: Y = ComplexHermitianEJA.real_unembed(Ye)
             sage: expected = (X*Y).trace().real()

-            sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
+            sage: actual = ComplexHermitianEJA.matrix_inner_product(Xe,Ye)

             sage: actual == expected
             True
 
         """
             sage: actual == expected
             True
 
         """

-        return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/2
+        return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/2

 
 
 
 

-class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
+class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra,
+                          ConcreteEuclideanJordanAlgebra):

     """
     The rank-n simple EJA consisting of complex Hermitian n-by-n
     matrices over the real numbers, the usual symmetric Jordan product,
     """
     The rank-n simple EJA consisting of complex Hermitian n-by-n
     matrices over the real numbers, the usual symmetric Jordan product,


@@ -1645,7 +1654,7 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
     The dimension of this algebra is `n^2`::
 
         sage: set_random_seed()
     The dimension of this algebra is `n^2`::
 
         sage: set_random_seed()

-        sage: n_max = ComplexHermitianEJA._max_test_case_size()
+        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
         sage: n = ZZ.random_element(1, n_max)
         sage: J = ComplexHermitianEJA(n)
         sage: J.dimension() == n^2


@@ -1656,9 +1665,9 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
         sage: set_random_seed()
         sage: J = ComplexHermitianEJA.random_instance()
         sage: x,y = J.random_elements(2)
         sage: set_random_seed()
         sage: J = ComplexHermitianEJA.random_instance()
         sage: x,y = J.random_elements(2)

-        sage: actual = (x*y).natural_representation()
-        sage: X = x.natural_representation()
-        sage: Y = y.natural_representation()
+        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: expected = (X*Y + Y*X)/2
         sage: actual == expected
         True


@@ -1670,24 +1679,10 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
         sage: ComplexHermitianEJA(2, prefix='z').gens()
         (z0, z1, z2, z3)
 
         sage: ComplexHermitianEJA(2, prefix='z').gens()
         (z0, z1, z2, z3)
 

-    Our natural basis is normalized with respect to the natural inner
-    product unless we specify otherwise::
-
-        sage: set_random_seed()
-        sage: J = ComplexHermitianEJA.random_instance()
-        sage: all( b.norm() == 1 for b in J.gens() )
-        True
-
-    Since our natural basis is normalized with respect to the natural
-    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::
+    We can construct the (trivial) algebra of rank zero::

 
 

-        sage: set_random_seed()
-        sage: x = ComplexHermitianEJA.random_instance().random_element()
-        sage: x.operator().matrix().is_symmetric()
-        True
+        sage: ComplexHermitianEJA(0)
+        Euclidean Jordan algebra of dimension 0 over Algebraic Real Field

 
     """
 
 
     """
 


@@ -1747,8 +1742,23 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
 
     def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n,field)
 
     def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n,field)

-        super(ComplexHermitianEJA,self).__init__(field, basis, n, **kwargs)
+        super(ComplexHermitianEJA,self).__init__(field,
+                                                 basis,
+                                                 check_axioms=False,
+                                                 **kwargs)
+        self.rank.set_cache(n)

 
 

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

 
 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
     @staticmethod
 
 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
     @staticmethod


@@ -1780,8 +1790,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         Embedding is a homomorphism (isomorphism, in fact)::
 
             sage: set_random_seed()
         Embedding is a homomorphism (isomorphism, in fact)::
 
             sage: set_random_seed()

-            sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
-            sage: n = ZZ.random_element(n_max)
+            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: Q = QuaternionAlgebra(QQ,-1,-1)
             sage: X = random_matrix(Q, n)
             sage: Y = random_matrix(Q, n)


@@ -1883,9 +1892,9 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
 
 
     @classmethod
 
 
     @classmethod

-    def natural_inner_product(cls,X,Y):
+    def matrix_inner_product(cls,X,Y):

         """
         """

-        Compute a natural inner product in this algebra directly from
+        Compute a matrix inner product in this algebra directly from

         its real embedding.
 
         SETUP::
         its real embedding.
 
         SETUP::


@@ -1900,22 +1909,22 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             sage: set_random_seed()
             sage: J = QuaternionHermitianEJA.random_instance()
             sage: x,y = J.random_elements(2)
             sage: set_random_seed()
             sage: J = QuaternionHermitianEJA.random_instance()
             sage: x,y = J.random_elements(2)

-            sage: Xe = x.natural_representation()
-            sage: Ye = y.natural_representation()
+            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: X = QuaternionHermitianEJA.real_unembed(Xe)
             sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
             sage: expected = (X*Y).trace().coefficient_tuple()[0]

-            sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
+            sage: actual = QuaternionHermitianEJA.matrix_inner_product(Xe,Ye)

             sage: actual == expected
             True
 
         """
             sage: actual == expected
             True
 
         """

-        return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/4
+        return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/4

 
 
 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
 
 
 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,

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


@@ -1940,7 +1949,7 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
     The dimension of this algebra is `2*n^2 - n`::
 
         sage: set_random_seed()
     The dimension of this algebra is `2*n^2 - n`::
 
         sage: set_random_seed()

-        sage: n_max = QuaternionHermitianEJA._max_test_case_size()
+        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
         sage: n = ZZ.random_element(1, n_max)
         sage: J = QuaternionHermitianEJA(n)
         sage: J.dimension() == 2*(n^2) - n


@@ -1951,9 +1960,9 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
         sage: set_random_seed()
         sage: J = QuaternionHermitianEJA.random_instance()
         sage: x,y = J.random_elements(2)
         sage: set_random_seed()
         sage: J = QuaternionHermitianEJA.random_instance()
         sage: x,y = J.random_elements(2)

-        sage: actual = (x*y).natural_representation()
-        sage: X = x.natural_representation()
-        sage: Y = y.natural_representation()
+        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: expected = (X*Y + Y*X)/2
         sage: actual == expected
         True


@@ -1965,24 +1974,10 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
         sage: QuaternionHermitianEJA(2, prefix='a').gens()
         (a0, a1, a2, a3, a4, a5)
 
         sage: QuaternionHermitianEJA(2, prefix='a').gens()
         (a0, a1, a2, a3, a4, a5)
 

-    Our natural basis is normalized with respect to the natural inner
-    product unless we specify otherwise::
-
-        sage: set_random_seed()
-        sage: J = QuaternionHermitianEJA.random_instance()
-        sage: all( b.norm() == 1 for b in J.gens() )
-        True
-
-    Since our natural basis is normalized with respect to the natural
-    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::
+    We can construct the (trivial) algebra of rank zero::

 
 

-        sage: set_random_seed()
-        sage: x = QuaternionHermitianEJA.random_instance().random_element()
-        sage: x.operator().matrix().is_symmetric()
-        True
+        sage: QuaternionHermitianEJA(0)
+        Euclidean Jordan algebra of dimension 0 over Algebraic Real Field

 
     """
     @classmethod
 
     """
     @classmethod


@@ -2043,17 +2038,125 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
 
     def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n,field)
 
     def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n,field)

-        super(QuaternionHermitianEJA,self).__init__(field, basis, n, **kwargs)
+        super(QuaternionHermitianEJA,self).__init__(field,
+                                                    basis,
+                                                    check_axioms=False,
+                                                    **kwargs)
+        self.rank.set_cache(n)
+
+    @staticmethod
+    def _max_random_instance_size():
+        r"""
+        The maximum rank of a random QuaternionHermitianEJA.
+        """
+        return 2 # Dimension 6
+
+    @classmethod
+    def random_instance(cls, field=AA, **kwargs):
+        """
+        Return a random instance of this type of algebra.
+        """
+        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        return cls(n, field, **kwargs)
+
+
+class HadamardEJA(RationalBasisEuclideanJordanAlgebra,
+                  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, field=AA, **kwargs):
+        V = VectorSpace(field, n)
+        mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
+                       for i in range(n) ]
+
+        # Inner products are real numbers and not algebra
+        # elements, so once we turn the algebra element
+        # into a vector in inner_product(), we never go
+        # back. As a result -- contrary to what we do with
+        # self._multiplication_table -- we store the inner
+        # product table as a plain old matrix and not as
+        # an algebra operator.
+        ip_table = matrix.identity(field,n)
+        self._inner_product_matrix = ip_table
+
+        super(HadamardEJA, self).__init__(field,
+                                          mult_table,
+                                          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, field=AA, **kwargs):
+        """
+        Return a random instance of this type of algebra.
+        """
+        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        return cls(n, field, **kwargs)

 
 
 
 

-class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
+class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
+                      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 =
     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 =

-    (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
-    symmetric positive-definite "bilinear form" matrix. It has
-    dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
-    when ``B`` is the identity matrix of order ``n-1``.
+    (<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::
 
 
     SETUP::
 


@@ -2065,16 +2168,32 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
     When no bilinear form is specified, the identity matrix is used,
     and the resulting algebra is the Jordan spin algebra::
 
     When no bilinear form is specified, the identity matrix is used,
     and the resulting algebra is the Jordan spin algebra::
 

-        sage: J0 = BilinearFormEJA(3)
+        sage: B = matrix.identity(AA,3)
+        sage: J0 = BilinearFormEJA(B)

         sage: J1 = JordanSpinEJA(3)
         sage: J0.multiplication_table() == J0.multiplication_table()
         True
 
         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::
 
     TESTS:
 
     We can create a zero-dimensional algebra::
 

-        sage: J = BilinearFormEJA(0)
+        sage: B = matrix.identity(AA,0)
+        sage: J = BilinearFormEJA(B)

         sage: J.basis()
         Finite family {}
 
         sage: J.basis()
         Finite family {}
 


@@ -2086,8 +2205,11 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
         sage: set_random_seed()
         sage: n = ZZ.random_element(5)
         sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
         sage: set_random_seed()
         sage: n = ZZ.random_element(5)
         sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')

-        sage: B = M.transpose()*M
-        sage: J = BilinearFormEJA(n, B=B)
+        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)

         sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
         sage: V = J.vector_space()
         sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
         sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
         sage: V = J.vector_space()
         sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))


@@ -2101,11 +2223,11 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
         sage: actual == expected
         True
     """
         sage: actual == expected
         True
     """

-    def __init__(self, n, field=AA, B=None, **kwargs):
-        if B is None:
-            self._B = matrix.identity(field, max(0,n-1))
-        else:
-            self._B = B
+    def __init__(self, B, field=AA, **kwargs):
+        n = B.nrows()
+
+        if not B.is_positive_definite():
+            raise ValueError("bilinear form is not positive-definite")

 
         V = VectorSpace(field, n)
         mult_table = [[V.zero() for j in range(n)] for i in range(n)]
 
         V = VectorSpace(field, n)
         mult_table = [[V.zero() for j in range(n)] for i in range(n)]


@@ -2117,51 +2239,66 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
                 xbar = x[1:]
                 y0 = y[0]
                 ybar = y[1:]
                 xbar = x[1:]
                 y0 = y[0]
                 ybar = y[1:]

-                z0 = x0*y0 + (self._B*xbar).inner_product(ybar)
+                z0 = (B*x).inner_product(y)

                 zbar = y0*xbar + x0*ybar
                 z = V([z0] + zbar.list())
                 mult_table[i][j] = z
 
                 zbar = y0*xbar + x0*ybar
                 z = V([z0] + zbar.list())
                 mult_table[i][j] = z
 

+        # Inner products are real numbers and not algebra
+        # elements, so once we turn the algebra element
+        # into a vector in inner_product(), we never go
+        # back. As a result -- contrary to what we do with
+        # self._multiplication_table -- we store the inner
+        # product table as a plain old matrix and not as
+        # an algebra operator.
+        ip_table = B
+        self._inner_product_matrix = ip_table
+
+        super(BilinearFormEJA, self).__init__(field,
+                                              mult_table,
+                                              check_axioms=False,
+                                              **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).
         # 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).

-        fdeja = super(BilinearFormEJA, self)
-        return fdeja.__init__(field, mult_table, rank=min(n,2), **kwargs)
-
-    def inner_product(self, x, y):
-        r"""
-        Half of the trace inner product.
-
-        This is defined so that the special case of the Jordan spin
-        algebra gets the usual inner product.
+        self.rank.set_cache(min(n,2))

 
 

-        SETUP::
+        if n == 0:
+            self.one.set_cache( self.zero() )
+        else:
+            self.one.set_cache( self.monomial(0) )

 
 

-            sage: from mjo.eja.eja_algebra import BilinearFormEJA
+    @staticmethod
+    def _max_random_instance_size():
+        r"""
+        The maximum dimension of a random BilinearFormEJA.
+        """
+        return 5

 
 

-        TESTS:
+    @classmethod
+    def random_instance(cls, field=AA, **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(field, n)
+            return cls(B, field, **kwargs)

 
 

-        Ensure that this is one-half of the trace inner-product when
-        the algebra isn't just the reals (when ``n`` isn't one). This
-        is in Faraut and Koranyi, and also my "On the symmetry..."
-        paper::
+        B11 = matrix.identity(field,1)
+        M = matrix.random(field, n-1)
+        I = matrix.identity(field, n-1)
+        alpha = field.zero()
+        while alpha.is_zero():
+            alpha = field.random_element().abs()
+        B22 = M.transpose()*M + alpha*I

 
 

-            sage: set_random_seed()
-            sage: n = ZZ.random_element(2,5)
-            sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
-            sage: B = M.transpose()*M
-            sage: J = BilinearFormEJA(n, B=B)
-            sage: x = J.random_element()
-            sage: y = J.random_element()
-            sage: x.inner_product(y) == (x*y).trace()/2
-            True
+        from sage.matrix.special import block_matrix
+        B = block_matrix(2,2, [ [B11,   ZZ(0) ],
+                                [ZZ(0), B22 ] ])

 
 

-        """
-        xvec = x.to_vector()
-        xbar = xvec[1:]
-        yvec = y.to_vector()
-        ybar = yvec[1:]
-        return x[0]*y[0] + (self._B*xbar).inner_product(ybar)
+        return cls(B, field, **kwargs)

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


@@ -2208,19 +2345,38 @@ class JordanSpinEJA(BilinearFormEJA):
             sage: set_random_seed()
             sage: J = JordanSpinEJA.random_instance()
             sage: x,y = J.random_elements(2)
             sage: set_random_seed()
             sage: J = JordanSpinEJA.random_instance()
             sage: x,y = J.random_elements(2)

-            sage: X = x.natural_representation()
-            sage: Y = y.natural_representation()
-            sage: x.inner_product(y) == J.natural_inner_product(X,Y)
+            sage: actual = x.inner_product(y)
+            sage: expected = x.to_vector().inner_product(y.to_vector())
+            sage: actual == expected

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

-        return super(JordanSpinEJA, self).__init__(n, field, **kwargs)
+        B = matrix.identity(field, n)
+        super(JordanSpinEJA, self).__init__(B, field, **kwargs)
+
+    @staticmethod
+    def _max_random_instance_size():
+        r"""
+        The maximum dimension of a random JordanSpinEJA.
+        """
+        return 5
+
+    @classmethod
+    def random_instance(cls, field=AA, **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, field, **kwargs)

 
 
 
 

-class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
+class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra,
+                 ConcreteEuclideanJordanAlgebra):

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


@@ -2251,7 +2407,255 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
     """
     def __init__(self, field=AA, **kwargs):
         mult_table = []
     """
     def __init__(self, field=AA, **kwargs):
         mult_table = []

-        fdeja = super(TrivialEJA, self)
+        self._inner_product_matrix = matrix(field,0)
+        super(TrivialEJA, self).__init__(field,
+                                         mult_table,
+                                         check_axioms=False,
+                                         **kwargs)

         # The rank is zero using my definition, namely the dimension of the
         # largest subalgebra generated by any element.
         # The rank is zero using my definition, namely the dimension of the
         # largest subalgebra generated by any element.

-        return fdeja.__init__(field, mult_table, rank=0, **kwargs)
+        self.rank.set_cache(0)
+        self.one.set_cache( self.zero() )
+
+    @classmethod
+    def random_instance(cls, field=AA, **kwargs):
+        # We don't take a "size" argument so the superclass method is
+        # inappropriate for us.
+        return cls(field, **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(AA)
+        sage: J2 = random_eja(QQ)
+        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(n) ]
+                       for i in range(n) ]
+        for i in range(n1):
+            for j in range(n1):
+                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(n2):
+                p = (J2.monomial(i)*J2.monomial(j)).to_vector()
+                mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
+
+        super(DirectSumEJA, self).__init__(field,
+                                           mult_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,QQ)
+            sage: J2 = JordanSpinEJA(3,QQ)
+            sage: J = DirectSumEJA(J1,J2)
+            sage: J.factors()
+            (Euclidean Jordan algebra of dimension 2 over Rational Field,
+             Euclidean Jordan algebra of dimension 3 over Rational Field)
+
+        """
+        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,QQ)
+            sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=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