]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: add random_instance() for BilinearFormEJA.
[sage.d.git] / mjo / eja / eja_algebra.py
index ed2e2587282e91d2672cbcbc851d37fcddf7ae74..8ca501d6ba944a539901d384a78a5c8a0f45a7ca 100644 (file)
@@ -57,11 +57,15 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
                  prefix='e',
                  category=None,
                  natural_basis=None,
-                 check=True):
+                 check_field=True,
+                 check_axioms=True):
         """
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import (JordanSpinEJA, random_eja)
+            sage: from mjo.eja.eja_algebra import (
+            ....:   FiniteDimensionalEuclideanJordanAlgebra,
+            ....:   JordanSpinEJA,
+            ....:   random_eja)
 
         EXAMPLES:
 
@@ -75,20 +79,33 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         TESTS:
 
-        The ``field`` we're given must be real::
+        The ``field`` we're given must be real with ``check_field=True``::
 
             sage: JordanSpinEJA(2,QQbar)
             Traceback (most recent call last):
             ...
-            ValueError: field is not real
+            ValueError: scalar field is not real
+
+        The multiplication table must be square with ``check_axioms=True``::
+
+            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()))
+            Traceback (most recent call last):
+            ...
+            ValueError: multiplication table is not square
 
         """
-        if check:
+        if check_field:
             if not field.is_subring(RR):
                 # Note: this does return true for the real algebraic
-                # field, and any quadratic field where we've specified
-                # a real embedding.
-                raise ValueError('field is not real')
+                # field, the rationals, and any quadratic field where
+                # we've specified a real embedding.
+                raise ValueError("scalar field is not real")
+
+        # The multiplication table had better be square
+        n = len(mult_table)
+        if check_axioms:
+            if not all( len(l) == n for l in mult_table ):
+                raise ValueError("multiplication table is not square")
 
         self._natural_basis = natural_basis
 
@@ -98,7 +115,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
         fda.__init__(field,
-                     range(len(mult_table)),
+                     range(n),
                      prefix=prefix,
                      category=category)
         self.print_options(bracket='')
@@ -114,6 +131,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             for ls in mult_table
         ]
 
+        if check_axioms:
+            if not self._is_commutative():
+                raise ValueError("algebra is not commutative")
+            if not self._is_jordanian():
+                raise ValueError("Jordan identity does not hold")
+            if not self._inner_product_is_associative():
+                raise ValueError("inner product is not associative")
 
     def _element_constructor_(self, elt):
         """
@@ -193,7 +217,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return self.from_vector(coords)
 
     @staticmethod
-    def _max_test_case_size():
+    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
@@ -209,7 +233,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         interpreted to be far less than the dimension) should override
         with a smaller number.
         """
-        return 5
+        raise NotImplementedError
 
     def _repr_(self):
         """
@@ -235,6 +259,67 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
     def product_on_basis(self, i, j):
         return self._multiplication_table[i][j]
 
+    def _is_commutative(self):
+        r"""
+        Whether or not this algebra's multiplication table is commutative.
+
+        This method should of course always return ``True``, unless
+        this algebra was constructed with ``check_axioms=False`` and
+        passed an invalid multiplication table.
+        """
+        return all( self.product_on_basis(i,j) == self.product_on_basis(i,j)
+                    for i in range(self.dimension())
+                    for j in range(self.dimension()) )
+
+    def _is_jordanian(self):
+        r"""
+        Whether or not this algebra's multiplication table respects the
+        Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
+
+        We only check one arrangement of `x` and `y`, so for a
+        ``True`` result to be truly true, you should also check
+        :meth:`_is_commutative`. This method should of course always
+        return ``True``, unless this algebra was constructed with
+        ``check_axioms=False`` and passed an invalid multiplication table.
+        """
+        return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
+                    ==
+                    (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j))
+                    for i in range(self.dimension())
+                    for j in range(self.dimension()) )
+
+    def _inner_product_is_associative(self):
+        r"""
+        Return whether or not this algebra's inner product `B` is
+        associative; that is, whether or not `B(xy,z) = B(x,yz)`.
+
+        This method should of course always return ``True``, unless
+        this algebra was constructed with ``check_axioms=False`` and
+        passed an invalid multiplication table.
+        """
+
+        # Used to check whether or not something is zero in an inexact
+        # ring. This number is sufficient to allow the construction of
+        # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
+        epsilon = 1e-16
+
+        for i in range(self.dimension()):
+            for j in range(self.dimension()):
+                for k in range(self.dimension()):
+                    x = self.monomial(i)
+                    y = self.monomial(j)
+                    z = self.monomial(k)
+                    diff = (x*y).inner_product(z) - x.inner_product(y*z)
+
+                    if self.base_ring().is_exact():
+                        if diff != 0:
+                            return False
+                    else:
+                        if diff.abs() > epsilon:
+                            return False
+
+        return True
+
     @cached_method
     def characteristic_polynomial_of(self):
         """
@@ -514,19 +599,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):
@@ -661,7 +747,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
                 J5 = eigspace
             else:
                 gens = tuple( self.from_vector(b) for b in eigspace.basis() )
-                subalg = FiniteDimensionalEuclideanJordanSubalgebra(self, gens)
+                subalg = FiniteDimensionalEuclideanJordanSubalgebra(self,
+                                                                    gens,
+                                                                    check_axioms=False)
                 if eigval == 0:
                     J0 = subalg
                 elif eigval == 1:
@@ -672,10 +760,61 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return (J0, J5, J1)
 
 
-    def random_elements(self, count):
+    def random_element(self, thorough=False):
+        r"""
+        Return a random element of this algebra.
+
+        Our algebra superclass method only returns a linear
+        combination of at most two basis elements. We instead
+        want the vector space "random element" method that
+        returns a more diverse selection.
+
+        INPUT:
+
+        - ``thorough`` -- (boolean; default False) whether or not we
+          should generate irrational coefficients for the random
+          element when our base ring is irrational; this slows the
+          algebra operations to a crawl, but any truly random method
+          should include them
+
+        """
+        # For a general base ring... maybe we can trust this to do the
+        # right thing? Unlikely, but.
+        V = self.vector_space()
+        v = V.random_element()
+
+        if self.base_ring() is AA:
+            # The "random element" method of the algebraic reals is
+            # stupid at the moment, and only returns integers between
+            # -2 and 2, inclusive:
+            #
+            #   https://trac.sagemath.org/ticket/30875
+            #
+            # Instead, we implement our own "random vector" method,
+            # and then coerce that into the algebra. We use the vector
+            # space degree here instead of the dimension because a
+            # subalgebra could (for example) be spanned by only two
+            # vectors, each with five coordinates.  We need to
+            # generate all five coordinates.
+            if thorough:
+                v *= QQbar.random_element().real()
+            else:
+                v *= QQ.random_element()
+
+        return self.from_vector(V.coordinate_vector(v))
+
+    def random_elements(self, count, thorough=False):
         """
         Return ``count`` random elements as a tuple.
 
+        INPUT:
+
+        - ``thorough`` -- (boolean; default False) whether or not we
+          should generate irrational coefficients for the random
+          elements when our base ring is irrational; this slows the
+          algebra operations to a crawl, but any truly random method
+          should include them
+
         SETUP::
 
             sage: from mjo.eja.eja_algebra import JordanSpinEJA
@@ -690,7 +829,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             True
 
         """
-        return tuple( self.random_element() for idx in range(count) )
+        return tuple( self.random_element(thorough)
+                      for idx in range(count) )
 
     @classmethod
     def random_instance(cls, field=AA, **kwargs):
@@ -700,12 +840,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         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)
+        n = ZZ.random_element(cls._max_random_instance_size() + 1)
         return cls(n, field, **kwargs)
 
     @cached_method
@@ -868,79 +1003,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
     Element = FiniteDimensionalEuclideanJordanAlgebraElement
 
 
-class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra):
-    """
-    Return the Euclidean Jordan Algebra corresponding to the set
-    `R^n` under the Hadamard product.
-
-    Note: this is nothing more than the Cartesian product of ``n``
-    copies of the spin algebra. Once Cartesian product algebras
-    are implemented, this can go.
-
-    SETUP::
-
-        sage: from mjo.eja.eja_algebra import HadamardEJA
-
-    EXAMPLES:
-
-    This multiplication table can be verified by hand::
-
-        sage: J = HadamardEJA(3)
-        sage: e0,e1,e2 = J.gens()
-        sage: e0*e0
-        e0
-        sage: e0*e1
-        0
-        sage: e0*e2
-        0
-        sage: e1*e1
-        e1
-        sage: e1*e2
-        0
-        sage: e2*e2
-        e2
-
-    TESTS:
-
-    We can change the generator prefix::
-
-        sage: HadamardEJA(3, prefix='r').gens()
-        (r0, r1, r2)
-
-    """
-    def __init__(self, n, field=AA, **kwargs):
-        V = VectorSpace(field, n)
-        mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
-                       for i in range(n) ]
-
-        fdeja = super(HadamardEJA, self)
-        fdeja.__init__(field, mult_table, **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`::
-
-            sage: set_random_seed()
-            sage: J = HadamardEJA.random_instance()
-            sage: x,y = J.random_elements(2)
-            sage: X = x.natural_representation()
-            sage: Y = y.natural_representation()
-            sage: x.inner_product(y) == J.natural_inner_product(X,Y)
-            True
-
-        """
-        return x.to_vector().inner_product(y.to_vector())
-
 
 def random_eja(field=AA):
     """
@@ -967,9 +1029,68 @@ def random_eja(field=AA):
 
 
 
+class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
+    r"""
+    Algebras whose basis consists of vectors with rational
+    entries. Equivalently, algebras whose multiplication tables
+    contain only rational coefficients.
+
+    When an EJA has a basis that can be made rational, we can speed up
+    the computation of its characteristic polynomial by doing it over
+    ``QQ``. All of the named EJA constructors that we provide fall
+    into this category.
+    """
+    @cached_method
+    def _charpoly_coefficients(self):
+        r"""
+        Override the parent method with something that tries to compute
+        over a faster (non-extension) field.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import JordanSpinEJA
+
+        EXAMPLES:
+
+        The base ring of the resulting polynomial coefficients is what
+        it should be, and not the rationals (unless the algebra was
+        already over the rationals)::
+
+            sage: J = JordanSpinEJA(3)
+            sage: J._charpoly_coefficients()
+            (X1^2 - X2^2 - X3^2, -2*X1)
+            sage: a0 = J._charpoly_coefficients()[0]
+            sage: J.base_ring()
+            Algebraic Real Field
+            sage: a0.base_ring()
+            Algebraic Real Field
+
+        """
+        if self.base_ring() is QQ:
+            # There's no need to construct *another* algebra over the
+            # rationals if this one is already over the rationals.
+            superclass = super(RationalBasisEuclideanJordanAlgebra, self)
+            return superclass._charpoly_coefficients()
+
+        mult_table = tuple(
+            map(lambda x: x.to_vector(), ls)
+            for ls in self._multiplication_table
+        )
+
+        # Do the computation over the rationals. The answer will be
+        # the same, because our basis coordinates are (essentially)
+        # rational.
+        J = FiniteDimensionalEuclideanJordanAlgebra(QQ,
+                                                    mult_table,
+                                                    check_field=False,
+                                                    check_axioms=False)
+        a = J._charpoly_coefficients()
+        return tuple(map(lambda x: x.change_ring(self.base_ring()), a))
+
+
 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
     @staticmethod
-    def _max_test_case_size():
+    def _max_random_instance_size():
         # Play it safe, since this will be squared and the underlying
         # field can have dimension 4 (quaternions) too.
         return 2
@@ -1003,9 +1124,10 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
 
         Qs = self.multiplication_table_from_matrix_basis(basis)
 
-        fdeja = super(MatrixEuclideanJordanAlgebra, self)
-        fdeja.__init__(field, Qs, natural_basis=basis, **kwargs)
-        return
+        super(MatrixEuclideanJordanAlgebra, self).__init__(field,
+                                                           Qs,
+                                                           natural_basis=basis,
+                                                           **kwargs)
 
 
     @cached_method
@@ -1014,40 +1136,44 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         Override the parent method with something that tries to compute
         over a faster (non-extension) field.
         """
-        if self._basis_normalizers is None:
-            # We didn't normalize, so assume that the basis we started
-            # with had entries in a nice field.
+        if self._basis_normalizers is None or self.base_ring() is QQ:
+            # We didn't normalize, or the basis we started with had
+            # entries in a nice field already. Just compute the thing.
             return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coefficients()
-        else:
-            basis = ( (b/n) for (b,n) in zip(self.natural_basis(),
-                                             self._basis_normalizers) )
-
-            # Do this over the rationals and convert back at the end.
-            # Only works because we know the entries of the basis are
-            # integers.
-            J = MatrixEuclideanJordanAlgebra(QQ,
-                                             basis,
-                                             normalize_basis=False)
-            a = J._charpoly_coefficients()
-
-            # Unfortunately, changing the basis does change the
-            # coefficients of the characteristic polynomial, but since
-            # these are really the coefficients of the "characteristic
-            # polynomial of" function, everything is still nice and
-            # unevaluated. It's therefore "obvious" how scaling the
-            # basis affects the coordinate variables X1, X2, et
-            # cetera. Scaling the first basis vector up by "n" adds a
-            # factor of 1/n into every "X1" term, for example. So here
-            # we simply undo the basis_normalizer scaling that we
-            # performed earlier.
-            #
-            # The a[0] access here is safe because trivial algebras
-            # won't have any basis normalizers and therefore won't
-            # make it to this "else" branch.
-            XS = a[0].parent().gens()
-            subs_dict = { XS[i]: self._basis_normalizers[i]*XS[i]
-                          for i in range(len(XS)) }
-            return tuple( a_i.subs(subs_dict) for a_i in a )
+
+        basis = ( (b/n) for (b,n) in zip(self.natural_basis(),
+                                         self._basis_normalizers) )
+
+        # Do this over the rationals and convert back at the end.
+        # Only works because we know the entries of the basis are
+        # integers. The argument ``check_axioms=False`` is required
+        # because the trace inner-product method for this
+        # class is a stub and can't actually be checked.
+        J = MatrixEuclideanJordanAlgebra(QQ,
+                                         basis,
+                                         normalize_basis=False,
+                                         check_field=False,
+                                         check_axioms=False)
+        a = J._charpoly_coefficients()
+
+        # Unfortunately, changing the basis does change the
+        # coefficients of the characteristic polynomial, but since
+        # these are really the coefficients of the "characteristic
+        # polynomial of" function, everything is still nice and
+        # unevaluated. It's therefore "obvious" how scaling the
+        # basis affects the coordinate variables X1, X2, et
+        # cetera. Scaling the first basis vector up by "n" adds a
+        # factor of 1/n into every "X1" term, for example. So here
+        # we simply undo the basis_normalizer scaling that we
+        # performed earlier.
+        #
+        # The a[0] access here is safe because trivial algebras
+        # won't have any basis normalizers and therefore won't
+        # make it to this "else" branch.
+        XS = a[0].parent().gens()
+        subs_dict = { XS[i]: self._basis_normalizers[i]*XS[i]
+                      for i in range(len(XS)) }
+        return tuple( a_i.subs(subs_dict) for a_i in a )
 
 
     @staticmethod
@@ -1114,16 +1240,11 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         Yu = cls.real_unembed(Y)
         tr = (Xu*Yu).trace()
 
-        if tr in RLF:
-            # It's real already.
-            return tr
-
-        # Otherwise, try the thing that works for complex numbers; and
-        # if that doesn't work, the thing that works for quaternions.
         try:
-            return tr.vector()[0] # real part, imag part is index 1
+            # Works in QQ, AA, RDF, et cetera.
+            return tr.real()
         except AttributeError:
-            # A quaternions doesn't have a vector() method, but does
+            # A quaternion doesn't have a real() method, but does
             # have coefficient_tuple() method that returns the
             # coefficients of 1, i, j, and k -- in that order.
             return tr.coefficient_tuple()[0]
@@ -1181,7 +1302,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
@@ -1264,13 +1385,16 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra):
 
 
     @staticmethod
-    def _max_test_case_size():
+    def _max_random_instance_size():
         return 4 # Dimension 10
 
 
     def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n, field)
-        super(RealSymmetricEJA, self).__init__(field, basis, **kwargs)
+        super(RealSymmetricEJA, self).__init__(field,
+                                               basis,
+                                               check_axioms=False,
+                                               **kwargs)
         self.rank.set_cache(n)
 
 
@@ -1307,7 +1431,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_max = ComplexMatrixEuclideanJordanAlgebra._max_random_instance_size()
             sage: n = ZZ.random_element(n_max)
             sage: F = QuadraticField(-1, 'I')
             sage: X = random_matrix(F, n)
@@ -1459,7 +1583,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
@@ -1566,7 +1690,10 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra):
 
     def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n,field)
-        super(ComplexHermitianEJA,self).__init__(field, basis, **kwargs)
+        super(ComplexHermitianEJA,self).__init__(field,
+                                                 basis,
+                                                 check_axioms=False,
+                                                 **kwargs)
         self.rank.set_cache(n)
 
 
@@ -1600,7 +1727,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_max = QuaternionMatrixEuclideanJordanAlgebra._max_random_instance_size()
             sage: n = ZZ.random_element(n_max)
             sage: Q = QuaternionAlgebra(QQ,-1,-1)
             sage: X = random_matrix(Q, n)
@@ -1759,7 +1886,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
@@ -1867,18 +1994,106 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra):
 
     def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n,field)
-        super(QuaternionHermitianEJA,self).__init__(field, basis, **kwargs)
+        super(QuaternionHermitianEJA,self).__init__(field,
+                                                    basis,
+                                                    check_axioms=False,
+                                                    **kwargs)
+        self.rank.set_cache(n)
+
+
+class HadamardEJA(RationalBasisEuclideanJordanAlgebra):
+    """
+    Return the Euclidean Jordan Algebra corresponding to the set
+    `R^n` under the Hadamard product.
+
+    Note: this is nothing more than the Cartesian product of ``n``
+    copies of the spin algebra. Once Cartesian product algebras
+    are implemented, this can go.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import HadamardEJA
+
+    EXAMPLES:
+
+    This multiplication table can be verified by hand::
+
+        sage: J = HadamardEJA(3)
+        sage: e0,e1,e2 = J.gens()
+        sage: e0*e0
+        e0
+        sage: e0*e1
+        0
+        sage: e0*e2
+        0
+        sage: e1*e1
+        e1
+        sage: e1*e2
+        0
+        sage: e2*e2
+        e2
+
+    TESTS:
+
+    We can change the generator prefix::
+
+        sage: HadamardEJA(3, prefix='r').gens()
+        (r0, r1, r2)
+
+    """
+    def __init__(self, n, field=AA, **kwargs):
+        V = VectorSpace(field, n)
+        mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
+                       for i in range(n) ]
+
+        super(HadamardEJA, self).__init__(field,
+                                          mult_table,
+                                          check_axioms=False,
+                                          **kwargs)
         self.rank.set_cache(n)
 
+    @staticmethod
+    def _max_random_instance_size():
+        return 5
+
+    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`::
+
+            sage: set_random_seed()
+            sage: J = HadamardEJA.random_instance()
+            sage: x,y = J.random_elements(2)
+            sage: X = x.natural_representation()
+            sage: Y = y.natural_representation()
+            sage: x.inner_product(y) == J.natural_inner_product(X,Y)
+            True
+
+        """
+        return x.to_vector().inner_product(y.to_vector())
+
 
-class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra):
+class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
     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::
 
@@ -1890,7 +2105,8 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra):
     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
@@ -1899,7 +2115,8 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra):
 
     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 {}
 
@@ -1911,8 +2128,11 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra):
         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()))
@@ -1926,11 +2146,12 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra):
         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 TypeError("matrix B is not positive-definite")
 
         V = VectorSpace(field, n)
         mult_table = [[V.zero() for j in range(n)] for i in range(n)]
@@ -1942,7 +2163,7 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra):
                 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
@@ -1950,10 +2171,41 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra):
         # The rank of this algebra is two, unless we're in a
         # one-dimensional ambient space (because the rank is bounded
         # by the ambient dimension).
-        fdeja = super(BilinearFormEJA, self)
-        fdeja.__init__(field, mult_table, **kwargs)
+        super(BilinearFormEJA, self).__init__(field,
+                                              mult_table,
+                                              check_axioms=False,
+                                              **kwargs)
         self.rank.set_cache(min(n,2))
 
+    @staticmethod
+    def _max_random_instance_size():
+        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 == 0:
+            # Special case needed since we use (n-1) below.
+            B = matrix.identity(field, 0)
+            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.
@@ -1973,21 +2225,15 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra):
         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):
@@ -2043,7 +2289,8 @@ 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.
-        return super(JordanSpinEJA, self).__init__(n, field, **kwargs)
+        B = matrix.identity(field, n)
+        super(JordanSpinEJA, self).__init__(B, field, **kwargs)
 
 
 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra):
@@ -2077,8 +2324,198 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra):
     """
     def __init__(self, field=AA, **kwargs):
         mult_table = []
-        fdeja = super(TrivialEJA, self)
+        super(TrivialEJA, self).__init__(field,
+                                         mult_table,
+                                         check_axioms=False,
+                                         **kwargs)
         # The rank is zero using my definition, namely the dimension of the
         # largest subalgebra generated by any element.
-        fdeja.__init__(field, mult_table, **kwargs)
         self.rank.set_cache(0)
+
+    @classmethod
+    def random_instance(cls, field=AA, **kwargs):
+        # We don't take a "size" argument so the superclass method is
+        # inappropriate for us.
+        return cls(field, **kwargs)
+
+class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
+    r"""
+    The external (orthogonal) direct sum of two other Euclidean Jordan
+    algebras. Essentially the Cartesian product of its two factors.
+    Every Euclidean Jordan algebra decomposes into an orthogonal
+    direct sum of simple Euclidean Jordan algebras, so no generality
+    is lost by providing only this construction.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import (HadamardEJA,
+        ....:                                  RealSymmetricEJA,
+        ....:                                  DirectSumEJA)
+
+    EXAMPLES::
+
+        sage: J1 = HadamardEJA(2)
+        sage: J2 = RealSymmetricEJA(3)
+        sage: J = DirectSumEJA(J1,J2)
+        sage: J.dimension()
+        8
+        sage: J.rank()
+        5
+
+    """
+    def __init__(self, J1, J2, field=AA, **kwargs):
+        self._factors = (J1, J2)
+        n1 = J1.dimension()
+        n2 = J2.dimension()
+        n = n1+n2
+        V = VectorSpace(field, n)
+        mult_table = [ [ V.zero() for j in range(n) ]
+                       for i in range(n) ]
+        for i in range(n1):
+            for j in range(n1):
+                p = (J1.monomial(i)*J1.monomial(j)).to_vector()
+                mult_table[i][j] = V(p.list() + [field.zero()]*n2)
+
+        for i in range(n2):
+            for j in range(n2):
+                p = (J2.monomial(i)*J2.monomial(j)).to_vector()
+                mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
+
+        super(DirectSumEJA, self).__init__(field,
+                                           mult_table,
+                                           check_axioms=False,
+                                           **kwargs)
+        self.rank.set_cache(J1.rank() + J2.rank())
+
+
+    def factors(self):
+        r"""
+        Return the pair of this algebra's factors.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            ....:                                  JordanSpinEJA,
+            ....:                                  DirectSumEJA)
+
+        EXAMPLES::
+
+            sage: J1 = HadamardEJA(2,QQ)
+            sage: J2 = JordanSpinEJA(3,QQ)
+            sage: J = DirectSumEJA(J1,J2)
+            sage: J.factors()
+            (Euclidean Jordan algebra of dimension 2 over Rational Field,
+             Euclidean Jordan algebra of dimension 3 over Rational Field)
+
+        """
+        return self._factors
+
+    def projections(self):
+        r"""
+        Return a pair of projections onto this algebra's factors.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            ....:                                  ComplexHermitianEJA,
+            ....:                                  DirectSumEJA)
+
+        EXAMPLES::
+
+            sage: J1 = JordanSpinEJA(2)
+            sage: J2 = ComplexHermitianEJA(2)
+            sage: J = DirectSumEJA(J1,J2)
+            sage: (pi_left, pi_right) = J.projections()
+            sage: J.one().to_vector()
+            (1, 0, 1, 0, 0, 1)
+            sage: pi_left(J.one()).to_vector()
+            (1, 0)
+            sage: pi_right(J.one()).to_vector()
+            (1, 0, 0, 1)
+
+        """
+        (J1,J2) = self.factors()
+        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:])
+        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 (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)
+
+        """
+        (J1,J2) = self.factors()
+        n = J1.dimension()
+        V_basis = self.vector_space().basis()
+        I1 = matrix.column(self.base_ring(), V_basis[:n])
+        I2 = matrix.column(self.base_ring(), V_basis[n:])
+        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))