]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: begin working on a new constructor.
[sage.d.git] / mjo / eja / eja_algebra.py
index 3f91ee3ae72c6281c9238804f146944d2c431f1e..44d83147bbeceb4acae52d573749545dec3c4b99 100644 (file)
@@ -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.
+
+
+SETUP::
+
+    sage: from mjo.eja.eja_algebra import random_eja
+
+EXAMPLES::
+
+    sage: random_eja()
+    Euclidean Jordan algebra of dimension...
+
 """
 
 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.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
-lazy_import('mjo.eja.eja_subalgebra',
-            'FiniteDimensionalEuclideanJordanSubalgebra')
+from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
 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):
             ...
-            ValueError: not a naturally-represented algebra element
+            ValueError: not an element of this algebra
 
         """
         return None
@@ -56,7 +64,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
                  mult_table,
                  prefix='e',
                  category=None,
-                 natural_basis=None,
+                 matrix_basis=None,
                  check_field=True,
                  check_axioms=True):
         """
@@ -107,7 +115,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             if not all( len(l) == n for l in mult_table ):
                 raise ValueError("multiplication table is not square")
 
-        self._natural_basis = natural_basis
+        self._matrix_basis = matrix_basis
 
         if category is None:
             category = MagmaticAlgebras(field).FiniteDimensional()
@@ -126,10 +134,15 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         # long run to have the multiplication table be in terms of
         # algebra elements. We do this after calling the superclass
         # constructor so that from_vector() knows what to do.
-        self._multiplication_table = [
-            list(map(lambda x: self.from_vector(x), ls))
-            for ls in mult_table
-        ]
+        self._multiplication_table = [ [ self.vector_space().zero()
+                                         for i in range(n) ]
+                                       for j in range(n) ]
+        # take advantage of symmetry
+        for i in range(n):
+            for j in range(i+1):
+                elt = self.from_vector(mult_table[i][j])
+                self._multiplication_table[i][j] = elt
+                self._multiplication_table[j][i] = elt
 
         if check_axioms:
             if not self._is_commutative():
@@ -141,7 +154,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
     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
@@ -169,13 +182,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             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
-        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()
@@ -186,9 +199,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             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.
@@ -200,40 +212,24 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             # 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
-        # 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.
-        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() )
 
-    @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.
+        try:
+            coords =  W.coordinate_vector(_mat2vec(elt))
+        except ArithmeticError:  # vector is not in free module
+            raise ValueError(msg)
 
-        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
+        return self.from_vector(coords)
 
     def _repr_(self):
         """
@@ -382,6 +378,28 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         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):
         """
@@ -393,7 +411,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import random_eja
+            sage: from mjo.eja.eja_algebra import (random_eja,
+            ....:                                  HadamardEJA,
+            ....:                                  BilinearFormEJA)
 
         EXAMPLES:
 
@@ -406,10 +426,34 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             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):
@@ -473,18 +517,33 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         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 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.)
+        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.
 
-        Note that this will always return a matrix. The standard basis
-        in `R^n` will be returned as `n`-by-`1` column matrices.
+        But, when we construct an algebra from a basis of matrices,
+        those matrix representations are lost in favor of coordinate
+        vectors *with respect to* that basis. We could eventually
+        convert back if we tried hard enough, but having the original
+        representations handy is valuable enough that we simply store
+        them and return them from this method.
+
+        Why implement this for non-matrix algebras? Avoiding special
+        cases for the :class:`BilinearFormEJA` pays with simplicity in
+        its own right. But mainly, we would like to be able to assume
+        that elements of a :class:`DirectSumEJA` can be displayed
+        nicely, without having to have special classes for direct sums
+        one of whose components was a matrix algebra.
 
         SETUP::
 
@@ -496,7 +555,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             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]
@@ -507,50 +566,38 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             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]
             )
-
         """
-        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 self._natural_basis
+            return self._matrix_basis
 
 
-    def natural_basis_space(self):
+    def matrix_space(self):
         """
-        Return the matrix space in which this algebra's natural basis
-        elements live.
+        Return the matrix space in which this algebra's elements live, if
+        we think of them as matrices (including column vectors of the
+        appropriate size).
 
         Generally this will be an `n`-by-`1` column-vector space,
         except when the algebra is trivial. There it's `n`-by-`n`
-        (where `n` is zero), to ensure that two elements of the
-        natural basis space (empty matrices) can be multiplied.
+        (where `n` is zero), to ensure that two elements of the matrix
+        space (empty matrices) can be multiplied.
+
+        Matrix algebras override this with something more useful.
         """
         if self.is_trivial():
             return MatrixSpace(self.base_ring(), 0)
-        elif self._natural_basis is None or len(self._natural_basis) == 0:
+        elif self._matrix_basis is None or len(self._matrix_basis) == 0:
             return MatrixSpace(self.base_ring(), self.dimension(), 1)
         else:
-            return self._natural_basis[0].matrix_space()
-
-
-    @staticmethod
-    def natural_inner_product(X,Y):
-        """
-        Compute the inner product of two naturally-represented elements.
-
-        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.
-        """
-        return (X.conjugate_transpose()*Y).trace()
+            return self._matrix_basis[0].matrix_space()
 
 
     @cached_method
@@ -588,6 +635,16 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             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.
@@ -599,19 +656,20 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         # 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.
-        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()) )
 
-        # 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):
@@ -666,22 +724,22 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             Vector space of degree 6 and dimension 2...
             sage: J1
             Euclidean Jordan algebra of dimension 3...
-            sage: J0.one().natural_representation()
+            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]).natural_representation()
+            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]).natural_representation()
+            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().natural_representation()
+            sage: J1.one().to_matrix()
             [1 0 0]
             [0 1 0]
             [0 0 0]
@@ -731,6 +789,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         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
@@ -831,21 +891,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return tuple( self.random_element(thorough)
                       for idx in range(count) )
 
-    @classmethod
-    def random_instance(cls, field=AA, **kwargs):
-        """
-        Return a random instance of this type of algebra.
-
-        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)
-
-        n = ZZ.random_element(cls._max_test_case_size() + 1)
-        return cls(n, field, **kwargs)
 
     @cached_method
     def _charpoly_coefficients(self):
@@ -854,8 +899,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         of" function.
         """
         n = self.dimension()
-        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()
         F = R.fraction_field()
 
@@ -950,38 +994,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         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: set_random_seed()    # long time
+            sage: J = random_eja()     # long time
+            sage: caches = J.rank()    # long time
+            sage: J.rank.clear_cache() # long time
+            sage: J.rank() == cached   # long time
+            True
 
-            sage: J = QuaternionHermitianEJA(2)
-            sage: J.rank.clear_cache()
-            sage: J.rank()
-            2
         """
         return len(self._charpoly_coefficients())
 
@@ -1006,32 +1025,75 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
     Element = FiniteDimensionalEuclideanJordanAlgebraElement
 
+class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlgebra):
+    def __init__(self,
+                 field,
+                 basis,
+                 jordan_product,
+                 inner_product,
+                 orthonormalize=True,
+                 prefix='e',
+                 category=None,
+                 check_field=True,
+                 check_axioms=True):
 
+        n = len(basis)
+        vector_basis = basis
 
-def random_eja(field=AA):
-    """
-    Return a "random" finite-dimensional Euclidean Jordan Algebra.
-
-    SETUP::
-
-        sage: from mjo.eja.eja_algebra import random_eja
-
-    TESTS::
+        from sage.matrix.matrix import is_Matrix
+        basis_is_matrices = False
 
-        sage: random_eja()
-        Euclidean Jordan algebra of dimension...
+        degree = 0
+        if n > 0:
+            if is_Matrix(basis[0]):
+                basis_is_matrices = True
+                vector_basis = tuple( map(_mat2vec,basis) )
+                degree = basis[0].nrows()**2
+            else:
+                degree = basis[0].degree()
 
-    """
-    classname = choice([TrivialEJA,
-                        HadamardEJA,
-                        JordanSpinEJA,
-                        RealSymmetricEJA,
-                        ComplexHermitianEJA,
-                        QuaternionHermitianEJA])
-    return classname.random_instance(field=field)
+        V = VectorSpace(field, degree)
 
+        self._deorthonormalization_matrix = matrix.identity(field,n)
+        if orthonormalize:
+            A = matrix(field, vector_basis)
+            # uh oh, this is only the "usual" inner product
+            Q,R = A.gram_schmidt(orthonormal=True)
+            self._deorthonormalization_matrix = R.inverse().transpose()
+            vector_basis = Q.rows()
+            W = V.span_of_basis( vector_basis )
+            if basis_is_matrices:
+                from mjo.eja.eja_utils import _vec2mat
+                basis = tuple( map(_vec2mat,vector_basis) )
 
+        mult_table = [ [0 for i in range(n)] for j in range(n) ]
+        ip_table = [ [0 for i in range(n)] for j in range(n) ]
 
+        for i in range(n):
+            for j in range(i+1):
+                # do another mat2vec because the multiplication
+                # table is in terms of vectors
+                elt = _mat2vec(jordan_product(basis[i],basis[j]))
+                elt = W.coordinate_vector(elt)
+                mult_table[i][j] = elt
+                mult_table[j][i] = elt
+                ip = inner_product(basis[i],basis[j])
+                ip_table[i][j] = ip
+                ip_table[j][i] = ip
+
+        self._inner_product_matrix = matrix(field,ip_table)
+
+        if basis_is_matrices:
+            for m in basis:
+                m.set_immutable()
+
+        super().__init__(field,
+                         mult_table,
+                         prefix,
+                         category,
+                         basis,
+                         check_field,
+                         check_axioms)
 
 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
     r"""
@@ -1077,7 +1139,7 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr
             return superclass._charpoly_coefficients()
 
         mult_table = tuple(
-            map(lambda x: x.to_vector(), ls)
+            tuple(map(lambda x: x.to_vector(), ls))
             for ls in self._multiplication_table
         )
 
@@ -1092,12 +1154,71 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr
         return tuple(map(lambda x: x.change_ring(self.base_ring()), a))
 
 
-class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
+class ConcreteEuclideanJordanAlgebra:
+    r"""
+    A class for the Euclidean Jordan algebras that we know by name.
+
+    These are the Jordan algebras whose basis, multiplication table,
+    rank, and so on are known a priori. More to the point, they are
+    the Euclidean Jordan algebras for which we are able to conjure up
+    a "random instance."
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
+
+    TESTS:
+
+    Our basis is normalized with respect to the algebra's inner
+    product, unless we specify otherwise::
+
+        sage: set_random_seed()
+        sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
+        sage: all( b.norm() == 1 for b in J.gens() )
+        True
+
+    Since our basis is orthonormal with respect to the algebra's inner
+    product, and since we know that this algebra is an EJA, any
+    left-multiplication operator's matrix will be symmetric because
+    natural->EJA basis representation is an isometry and within the
+    EJA the operator is self-adjoint by the Jordan axiom::
+
+        sage: set_random_seed()
+        sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
+        sage: x = J.random_element()
+        sage: x.operator().is_self_adjoint()
+        True
+    """
+
     @staticmethod
-    def _max_test_case_size():
-        # Play it safe, since this will be squared and the underlying
-        # field can have dimension 4 (quaternions) too.
-        return 2
+    def _max_random_instance_size():
+        """
+        Return an integer "size" that is an upper bound on the size of
+        this algebra when it is used in a random test
+        case. Unfortunately, the term "size" is ambiguous -- when
+        dealing with `R^n` under either the Hadamard or Jordan spin
+        product, the "size" refers to the dimension `n`. When dealing
+        with a matrix algebra (real symmetric or complex/quaternion
+        Hermitian), it refers to the size of the matrix, which is far
+        less than the dimension of the underlying vector space.
+
+        This method must be implemented in each subclass.
+        """
+        raise NotImplementedError
+
+    @classmethod
+    def random_instance(cls, 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):
 
     def __init__(self, field, basis, normalize_basis=True, **kwargs):
         """
@@ -1112,7 +1233,12 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         # time to ensure that it isn't a generator expression.
         basis = tuple(basis)
 
-        if len(basis) > 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.
@@ -1123,16 +1249,51 @@ 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(
-                ~(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))
 
-        Qs = self.multiplication_table_from_matrix_basis(basis)
+        # 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 )
+        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))
+
+                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,
-                                                           Qs,
-                                                           natural_basis=basis,
+                                                           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):
@@ -1145,7 +1306,7 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
             # 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.natural_basis(),
+        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.
@@ -1180,39 +1341,6 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         return tuple( a_i.subs(subs_dict) for a_i in a )
 
 
-    @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.
-        if len(basis) == 0:
-            return []
-
-        field = basis[0].base_ring()
-        dimension = basis[0].nrows()
-
-        V = VectorSpace(field, dimension**2)
-        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):
-                mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
-                mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
-
-        return mult_table
-
-
     @staticmethod
     def real_embed(M):
         """
@@ -1237,9 +1365,8 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         """
         raise NotImplementedError
 
-
     @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()
@@ -1272,7 +1399,8 @@ class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         return M
 
 
-class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra):
+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
@@ -1306,7 +1434,7 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra):
     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
@@ -1317,9 +1445,9 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra):
         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
@@ -1331,25 +1459,6 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra):
         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::
-
-        sage: set_random_seed()
-        sage: x = RealSymmetricEJA.random_instance().random_element()
-        sage: x.operator().matrix().is_symmetric()
-        True
-
     We can construct the (trivial) algebra of rank zero::
 
         sage: RealSymmetricEJA(0)
@@ -1389,9 +1498,16 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra):
 
 
     @staticmethod
-    def _max_test_case_size():
+    def _max_random_instance_size():
         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)
@@ -1435,8 +1551,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         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)
@@ -1531,9 +1646,9 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
 
 
     @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::
@@ -1548,20 +1663,21 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             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: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
+            sage: actual = ComplexHermitianEJA.matrix_inner_product(Xe,Ye)
             sage: actual == expected
             True
 
         """
-        return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/2
+        return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/2
 
 
-class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra):
+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,
@@ -1587,7 +1703,7 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra):
     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
@@ -1598,9 +1714,9 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra):
         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
@@ -1612,25 +1728,6 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra):
         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::
-
-        sage: set_random_seed()
-        sage: x = ComplexHermitianEJA.random_instance().random_element()
-        sage: x.operator().matrix().is_symmetric()
-        True
-
     We can construct the (trivial) algebra of rank zero::
 
         sage: ComplexHermitianEJA(0)
@@ -1700,6 +1797,17 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra):
                                                  **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
@@ -1731,8 +1839,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         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)
@@ -1834,9 +1941,9 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
 
 
     @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::
@@ -1851,21 +1958,22 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             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: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
+            sage: actual = QuaternionHermitianEJA.matrix_inner_product(Xe,Ye)
             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,
+                             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
@@ -1890,7 +1998,7 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra):
     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
@@ -1901,9 +2009,9 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra):
         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
@@ -1915,25 +2023,6 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra):
         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::
-
-        sage: set_random_seed()
-        sage: x = QuaternionHermitianEJA.random_instance().random_element()
-        sage: x.operator().matrix().is_symmetric()
-        True
-
     We can construct the (trivial) algebra of rank zero::
 
         sage: QuaternionHermitianEJA(0)
@@ -2004,8 +2093,24 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra):
                                                     **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):
+class HadamardEJA(RationalBasisEuclideanJordanAlgebra,
+                  ConcreteEuclideanJordanAlgebra):
     """
     Return the Euclidean Jordan Algebra corresponding to the set
     `R^n` under the Hadamard product.
@@ -2050,45 +2155,57 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra):
         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)
 
-    def inner_product(self, x, y):
-        """
-        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`::
+        if n == 0:
+            self.one.set_cache( self.zero() )
+        else:
+            self.one.set_cache( sum(self.gens()) )
 
-            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
+    @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.
         """
-        return x.to_vector().inner_product(y.to_vector())
+        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        return cls(n, field, **kwargs)
 
 
-class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
+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 =
-    (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::
 
@@ -2100,16 +2217,32 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
     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
 
+    An error is raised if the matrix `B` does not correspond to a
+    positive-definite bilinear form::
+
+        sage: B = matrix.random(QQ,2,3)
+        sage: J = BilinearFormEJA(B)
+        Traceback (most recent call last):
+        ...
+        ValueError: bilinear form is not positive-definite
+        sage: B = matrix.zero(QQ,3)
+        sage: J = BilinearFormEJA(B)
+        Traceback (most recent call last):
+        ...
+        ValueError: bilinear form is not positive-definite
+
     TESTS:
 
     We can create a zero-dimensional algebra::
 
-        sage: J = BilinearFormEJA(0)
+        sage: B = matrix.identity(AA,0)
+        sage: J = BilinearFormEJA(B)
         sage: J.basis()
         Finite family {}
 
@@ -2121,8 +2254,11 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
         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()))
@@ -2136,11 +2272,11 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
         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)]
@@ -2152,54 +2288,66 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
                 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
 
-        # 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).
+        # 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)
-        self.rank.set_cache(min(n,2))
 
-    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.
+        # The rank of this algebra is two, unless we're in a
+        # one-dimensional ambient space (because the rank is bounded
+        # by the ambient dimension).
+        self.rank.set_cache(min(n,2))
 
-        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):
@@ -2246,19 +2394,38 @@ class JordanSpinEJA(BilinearFormEJA):
             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.
-        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
 
-class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra):
+    @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,
+                 ConcreteEuclideanJordanAlgebra):
     """
     The trivial Euclidean Jordan algebra consisting of only a zero element.
 
@@ -2289,6 +2456,7 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra):
     """
     def __init__(self, field=AA, **kwargs):
         mult_table = []
+        self._inner_product_matrix = matrix(field,0)
         super(TrivialEJA, self).__init__(field,
                                          mult_table,
                                          check_axioms=False,
@@ -2296,7 +2464,13 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra):
         # The rank is zero using my definition, namely the dimension of the
         # largest subalgebra generated by any element.
         self.rank.set_cache(0)
+        self.one.set_cache( self.zero() )
 
+    @classmethod
+    def random_instance(cls, field=AA, **kwargs):
+        # We don't take a "size" argument so the superclass method is
+        # inappropriate for us.
+        return cls(field, **kwargs)
 
 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
     r"""
@@ -2308,7 +2482,8 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
 
     SETUP::
 
-        sage: from mjo.eja.eja_algebra import (HadamardEJA,
+        sage: from mjo.eja.eja_algebra import (random_eja,
+        ....:                                  HadamardEJA,
         ....:                                  RealSymmetricEJA,
         ....:                                  DirectSumEJA)
 
@@ -2322,8 +2497,25 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
         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, field=AA, **kwargs):
+    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()
@@ -2395,7 +2587,124 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
 
         """
         (J1,J2) = self.factors()
-        n = J1.dimension()
-        pi_left  = lambda x: J1.from_vector(x.to_vector()[:n])
-        pi_right = lambda x: J2.from_vector(x.to_vector()[n:])
+        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