]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: factor out some common tests.
[sage.d.git] / mjo / eja / eja_algebra.py
index 183e4cb978c233dbac31e321e855978309b27687..3b5828fed55098dfd3f4c95bf2354ac4e38efc87 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
@@ -14,7 +25,6 @@ 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,
@@ -216,25 +226,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         coords =  W.coordinate_vector(_mat2vec(elt))
         return self.from_vector(coords)
 
-    @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
-
     def _repr_(self):
         """
         Return a string representation of ``self``.
@@ -588,6 +579,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 +600,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):
@@ -831,21 +833,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):
@@ -1007,32 +994,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
     Element = FiniteDimensionalEuclideanJordanAlgebraElement
 
 
-
-def random_eja(field=AA):
-    """
-    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...
-
-    """
-    classname = choice([TrivialEJA,
-                        HadamardEJA,
-                        JordanSpinEJA,
-                        RealSymmetricEJA,
-                        ComplexHermitianEJA,
-                        QuaternionHermitianEJA])
-    return classname.random_instance(field=field)
-
-
-
-
 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
     r"""
     Algebras whose basis consists of vectors with rational
@@ -1092,12 +1053,72 @@ 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 natural basis is normalized with respect to the natural 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 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: J = ConcreteEuclideanJordanAlgebra.random_instance()
+        sage: x = J.random_element()
+        sage: x.operator().matrix().is_symmetric()
+        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 +1133,8 @@ 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)
+        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.
@@ -1133,6 +1155,14 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
                                                            natural_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):
@@ -1272,7 +1302,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 +1337,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
@@ -1331,25 +1362,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 +1401,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 +1454,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)
@@ -1561,7 +1579,8 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         return RealMatrixEuclideanJordanAlgebra.natural_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 +1606,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
@@ -1612,25 +1631,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 +1700,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 +1742,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)
@@ -1864,8 +1874,9 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         return RealMatrixEuclideanJordanAlgebra.natural_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 +1901,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
@@ -1915,25 +1926,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 +1996,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.
@@ -2056,6 +2064,27 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra):
                                           **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)
+
+
     def inner_product(self, x, y):
         """
         Faster to reimplement than to use natural representations.
@@ -2081,14 +2110,20 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra):
         return x.to_vector().inner_product(y.to_vector())
 
 
-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 +2135,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 +2172,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 +2190,12 @@ 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):
+        self._B = B
+        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,7 +2207,7 @@ 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
@@ -2166,6 +2221,42 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
                                               **kwargs)
         self.rank.set_cache(min(n,2))
 
+        if n == 0:
+            self.one.set_cache( self.zero() )
+        else:
+            self.one.set_cache( self.monomial(0) )
+
+    @staticmethod
+    def _max_random_instance_size():
+        r"""
+        The maximum dimension of a random BilinearFormEJA.
+        """
+        return 5
+
+    @classmethod
+    def random_instance(cls, 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)
+
+        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
+
+        from sage.matrix.special import block_matrix
+        B = block_matrix(2,2, [ [B11,   ZZ(0) ],
+                                [ZZ(0), B22 ] ])
+
+        return cls(B, field, **kwargs)
+
     def inner_product(self, x, y):
         r"""
         Half of the trace inner product.
@@ -2185,21 +2276,15 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
         paper::
 
             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: J = BilinearFormEJA.random_instance()
+            sage: n = J.dimension()
             sage: x = J.random_element()
             sage: y = J.random_element()
-            sage: x.inner_product(y) == (x*y).trace()/2
+            sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
             True
 
         """
-        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 (self._B*x.to_vector()).inner_product(y.to_vector())
 
 
 class JordanSpinEJA(BilinearFormEJA):
@@ -2255,10 +2340,29 @@ class JordanSpinEJA(BilinearFormEJA):
     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.
 
@@ -2296,7 +2400,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"""
@@ -2440,3 +2550,46 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
         iota_left = lambda x: self.from_vector(I1*x.to_vector())
         iota_right = lambda x: self.from_vector(I2*+x.to_vector())
         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)
+            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