]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: eliminate the special element subalgebra class.
[sage.d.git] / mjo / eja / eja_algebra.py
index 1250fbda7981aba5474af0a3d2615c969bc9d1e1..83ec50e5a215417811fe09d72c61251081772160 100644 (file)
@@ -113,60 +113,36 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # we see in things like x = 1*e1 + 2*e2.
         vector_basis = basis
 
-        from sage.structure.element import is_Matrix
-        basis_is_matrices = False
-
         degree = 0
         if n > 0:
-            if is_Matrix(basis[0]):
-                if basis[0].is_square():
-                    # TODO: this ugly is_square() hack works around the problem
-                    # of passing to_matrix()ed vectors in as the basis from a
-                    # subalgebra. They aren't REALLY matrices, at least not of
-                    # the type that we assume here... Ugh.
-                    basis_is_matrices = True
-                    from mjo.eja.eja_utils import _vec2mat
-                    vector_basis = tuple( map(_mat2vec,basis) )
-                    degree = basis[0].nrows()**2
-                else:
-                    # convert from column matrices to vectors, yuck
-                    basis = tuple( map(_mat2vec,basis) )
-                    vector_basis = basis
-                    degree = basis[0].degree()
-            else:
-                degree = basis[0].degree()
+            # Works on both column and square matrices...
+            degree = len(basis[0].list())
 
-        # Build an ambient space that fits...
+        # Build an ambient space that fits our matrix basis when
+        # written out as "long vectors."
         V = VectorSpace(field, degree)
 
-        # We overwrite the name "vector_basis" in a second, but never modify it
-        # in place, to this effectively makes a copy of it.
-        deortho_vector_basis = vector_basis
+        # The matrix that will hole the orthonormal -> unorthonormal
+        # coordinate transformation.
         self._deortho_matrix = None
 
         if orthonormalize:
-            from mjo.eja.eja_utils import gram_schmidt
-            if basis_is_matrices:
-                vector_ip = lambda x,y: inner_product(_vec2mat(x), _vec2mat(y))
-                vector_basis = gram_schmidt(vector_basis, vector_ip)
-            else:
-                vector_basis = gram_schmidt(vector_basis, inner_product)
+            # Save a copy of the un-orthonormalized basis for later.
+            # Convert it to ambient V (vector) coordinates while we're
+            # at it, because we'd have to do it later anyway.
+            deortho_vector_basis = tuple( V(b.list()) for b in basis )
 
-            # Normalize the "matrix" basis, too!
-            basis = vector_basis
-
-            if basis_is_matrices:
-                basis = tuple( map(_vec2mat,basis) )
+            from mjo.eja.eja_utils import gram_schmidt
+            basis = tuple(gram_schmidt(basis, inner_product))
 
-        # Save the matrix "basis" for later... this is the last time we'll
-        # reference it in this constructor.
-        if basis_is_matrices:
-            self._matrix_basis = basis
-        else:
-            MS = MatrixSpace(self.base_ring(), degree, 1)
-            self._matrix_basis = tuple( MS(b) for b in basis )
+        # Save the (possibly orthonormalized) matrix basis for
+        # later...
+        self._matrix_basis = basis
 
-        # Now create the vector space for the algebra...
+        # Now create the vector space for the algebra, which will have
+        # its own set of non-ambient coordinates (in terms of the
+        # supplied basis).
+        vector_basis = tuple( V(b.list()) for b in basis )
         W = V.span_of_basis( vector_basis, check=check_axioms)
 
         if orthonormalize:
@@ -182,40 +158,33 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         # Now we actually compute the multiplication and inner-product
         # tables/matrices using the possibly-orthonormalized basis.
-        self._inner_product_matrix = matrix.zero(field, n)
-        self._multiplication_table = [ [0 for j in range(i+1)] for i in range(n) ]
+        self._inner_product_matrix = matrix.identity(field, n)
+        self._multiplication_table = [ [0 for j in range(i+1)]
+                                       for i in range(n) ]
 
-        print("vector_basis:")
-        print(vector_basis)
         # Note: the Jordan and inner-products are defined in terms
         # of the ambient basis. It's important that their arguments
         # are in ambient coordinates as well.
         for i in range(n):
             for j in range(i+1):
                 # ortho basis w.r.t. ambient coords
-                q_i = vector_basis[i]
-                q_j = vector_basis[j]
-
-                if basis_is_matrices:
-                    q_i = _vec2mat(q_i)
-                    q_j = _vec2mat(q_j)
+                q_i = basis[i]
+                q_j = basis[j]
 
+                # The jordan product returns a matrixy answer, so we
+                # have to convert it to the algebra coordinates.
                 elt = jordan_product(q_i, q_j)
-                ip = inner_product(q_i, q_j)
-
-                if basis_is_matrices:
-                    # do another mat2vec because the multiplication
-                    # table is in terms of vectors
-                    elt = _mat2vec(elt)
-
-                # TODO: the jordan product turns things back into
-                # matrices here even if they're supposed to be
-                # vectors. ugh. Can we get rid of vectors all together
-                # please?
-                elt = W.coordinate_vector(elt)
+                elt = W.coordinate_vector(V(elt.list()))
                 self._multiplication_table[i][j] = self.from_vector(elt)
-                self._inner_product_matrix[i,j] = ip
-                self._inner_product_matrix[j,i] = ip
+
+                if not orthonormalize:
+                    # If we're orthonormalizing the basis with respect
+                    # to an inner-product, then the inner-product
+                    # matrix with respect to the resulting basis is
+                    # just going to be the identity.
+                    ip = inner_product(q_i, q_j)
+                    self._inner_product_matrix[i,j] = ip
+                    self._inner_product_matrix[j,i] = ip
 
         self._inner_product_matrix._cache = {'hermitian': True}
         self._inner_product_matrix.set_immutable()
@@ -351,7 +320,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         This method should of course always return ``True``, unless
         this algebra was constructed with ``check_axioms=False`` and
-        passed an invalid multiplication table.
+        passed an invalid Jordan or inner-product.
         """
 
         # Used to check whether or not something is zero in an inexact
@@ -781,23 +750,57 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: from mjo.eja.eja_algebra import (HadamardEJA,
             ....:                                  random_eja)
 
-        EXAMPLES::
+        EXAMPLES:
+
+        We can compute unit element in the Hadamard EJA::
 
             sage: J = HadamardEJA(5)
             sage: J.one()
             e0 + e1 + e2 + e3 + e4
 
+        The unit element in the Hadamard EJA is inherited in the
+        subalgebras generated by its elements::
+
+            sage: J = HadamardEJA(5)
+            sage: J.one()
+            e0 + e1 + e2 + e3 + e4
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
+            sage: A.one()
+            f0
+            sage: A.one().superalgebra_element()
+            e0 + e1 + e2 + e3 + e4
+
         TESTS:
 
-        The identity element acts like the identity::
+        The identity element acts like the identity, regardless of
+        whether or not we orthonormalize::
 
             sage: set_random_seed()
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: J.one()*x == x and x*J.one() == x
             True
+            sage: A = x.subalgebra_generated_by()
+            sage: y = A.random_element()
+            sage: A.one()*y == y and y*A.one() == y
+            True
+
+        ::
+
+            sage: set_random_seed()
+            sage: J = random_eja(field=QQ, orthonormalize=False)
+            sage: x = J.random_element()
+            sage: J.one()*x == x and x*J.one() == x
+            True
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
+            sage: y = A.random_element()
+            sage: A.one()*y == y and y*A.one() == y
+            True
 
-        The matrix of the unit element's operator is the identity::
+        The matrix of the unit element's operator is the identity,
+        regardless of the base field and whether or not we
+        orthonormalize::
 
             sage: set_random_seed()
             sage: J = random_eja()
@@ -805,6 +808,27 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: expected = matrix.identity(J.base_ring(), J.dimension())
             sage: actual == expected
             True
+            sage: x = J.random_element()
+            sage: A = x.subalgebra_generated_by()
+            sage: actual = A.one().operator().matrix()
+            sage: expected = matrix.identity(A.base_ring(), A.dimension())
+            sage: actual == expected
+            True
+
+        ::
+
+            sage: set_random_seed()
+            sage: J = random_eja(field=QQ, orthonormalize=False)
+            sage: actual = J.one().operator().matrix()
+            sage: expected = matrix.identity(J.base_ring(), J.dimension())
+            sage: actual == expected
+            True
+            sage: x = J.random_element()
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
+            sage: actual = A.one().operator().matrix()
+            sage: expected = matrix.identity(A.base_ring(), A.dimension())
+            sage: actual == expected
+            True
 
         Ensure that the cached unit element (often precomputed by
         hand) agrees with the computed one::
@@ -816,6 +840,15 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: J.one() == cached
             True
 
+        ::
+
+            sage: set_random_seed()
+            sage: J = random_eja(field=QQ, orthonormalize=False)
+            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.
@@ -1223,9 +1256,7 @@ class RationalBasisEJA(FiniteDimensionalEJA):
                  jordan_product,
                  inner_product,
                  field=AA,
-                 orthonormalize=True,
                  check_field=True,
-                 check_axioms=True,
                  **kwargs):
 
         if check_field:
@@ -1234,6 +1265,7 @@ class RationalBasisEJA(FiniteDimensionalEJA):
             if not all( all(b_i in QQ for b_i in b.list()) for b in basis ):
                 raise TypeError("basis not rational")
 
+        self._rational_algebra = None
         if field is not QQ:
             # There's no point in constructing the extra algebra if this
             # one is already rational.
@@ -1248,15 +1280,13 @@ class RationalBasisEJA(FiniteDimensionalEJA):
                                        field=QQ,
                                        orthonormalize=False,
                                        check_field=False,
-                                       check_axioms=False,
-                                       **kwargs)
+                                       check_axioms=False)
 
         super().__init__(basis,
                          jordan_product,
                          inner_product,
                          field=field,
                          check_field=check_field,
-                         check_axioms=check_axioms,
                          **kwargs)
 
     @cached_method
@@ -1296,7 +1326,14 @@ class RationalBasisEJA(FiniteDimensionalEJA):
         a = ( a_i.change_ring(self.base_ring())
               for a_i in self._rational_algebra._charpoly_coefficients() )
 
-        # Now convert the coordinate variables back to the
+        if self._deortho_matrix is None:
+            # This can happen if our base ring was, say, AA and we
+            # chose not to (or didn't need to) orthonormalize. It's
+            # still faster to do the computations over QQ even if
+            # the numbers in the boxes stay the same.
+            return tuple(a)
+
+        # Otherwise, convert the coordinate variables back to the
         # deorthonormalized ones.
         R = self.coordinate_polynomial_ring()
         from sage.modules.free_module_element import vector
@@ -2222,11 +2259,16 @@ class HadamardEJA(ConcreteEJA):
 
     """
     def __init__(self, n, **kwargs):
-        def jordan_product(x,y):
-            P = x.parent()
-            return P(tuple( xi*yi for (xi,yi) in zip(x,y) ))
-        def inner_product(x,y):
-            return x.inner_product(y)
+        if n == 0:
+            jordan_product = lambda x,y: x
+            inner_product = lambda x,y: x
+        else:
+            def jordan_product(x,y):
+                P = x.parent()
+                return P( xi*yi for (xi,yi) in zip(x,y) )
+
+            def inner_product(x,y):
+                return (x.T*y)[0,0]
 
         # New defaults for keyword arguments. Don't orthonormalize
         # because our basis is already orthonormal with respect to our
@@ -2237,12 +2279,8 @@ class HadamardEJA(ConcreteEJA):
         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
-
-        standard_basis = FreeModule(ZZ, n).basis()
-        super(HadamardEJA, self).__init__(standard_basis,
-                                          jordan_product,
-                                          inner_product,
-                                          **kwargs)
+        column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )
+        super().__init__(column_basis, jordan_product, inner_product, **kwargs)
         self.rank.set_cache(n)
 
         if n == 0:
@@ -2346,31 +2384,38 @@ class BilinearFormEJA(ConcreteEJA):
         ....:              for j in range(n-1) ]
         sage: actual == expected
         True
+
     """
     def __init__(self, B, **kwargs):
-        if not B.is_positive_definite():
-            raise ValueError("bilinear form is not positive-definite")
+        # The matrix "B" is supplied by the user in most cases,
+        # so it makes sense to check whether or not its positive-
+        # definite unless we are specifically asked not to...
+        if ("check_axioms" not in kwargs) or kwargs["check_axioms"]:
+            if not B.is_positive_definite():
+                raise ValueError("bilinear form is not positive-definite")
+
+        # However, all of the other data for this EJA is computed
+        # by us in manner that guarantees the axioms are
+        # satisfied. So, again, unless we are specifically asked to
+        # verify things, we'll skip the rest of the checks.
+        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
         def inner_product(x,y):
-            return (B*x).inner_product(y)
+            return (y.T*B*x)[0,0]
 
         def jordan_product(x,y):
             P = x.parent()
-            x0 = x[0]
-            xbar = x[1:]
-            y0 = y[0]
-            ybar = y[1:]
-            z0 = inner_product(x,y)
+            x0 = x[0,0]
+            xbar = x[1:,0]
+            y0 = y[0,0]
+            ybar = y[1:,0]
+            z0 = inner_product(y,x)
             zbar = y0*xbar + x0*ybar
-            return P((z0,) + tuple(zbar))
-
-        # We know this is a valid EJA, but will double-check
-        # if the user passes check_axioms=True.
-        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
+            return P([z0] + zbar.list())
 
         n = B.nrows()
-        standard_basis = FreeModule(ZZ, n).basis()
-        super(BilinearFormEJA, self).__init__(standard_basis,
+        column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )
+        super(BilinearFormEJA, self).__init__(column_basis,
                                               jordan_product,
                                               inner_product,
                                               **kwargs)
@@ -2551,243 +2596,243 @@ class TrivialEJA(ConcreteEJA):
         # inappropriate for us.
         return cls(**kwargs)
 
-class DirectSumEJA(ConcreteEJA):
-    r"""
-    The external (orthogonal) direct sum of two other Euclidean Jordan
-    algebras. Essentially the Cartesian product of its two factors.
-    Every Euclidean Jordan algebra decomposes into an orthogonal
-    direct sum of simple Euclidean Jordan algebras, so no generality
-    is lost by providing only this construction.
-
-    SETUP::
-
-        sage: from mjo.eja.eja_algebra import (random_eja,
-        ....:                                  HadamardEJA,
-        ....:                                  RealSymmetricEJA,
-        ....:                                  DirectSumEJA)
-
-    EXAMPLES::
-
-        sage: J1 = HadamardEJA(2)
-        sage: J2 = RealSymmetricEJA(3)
-        sage: J = DirectSumEJA(J1,J2)
-        sage: J.dimension()
-        8
-        sage: J.rank()
-        5
-
-    TESTS:
-
-    The external direct sum construction is only valid when the two factors
-    have the same base ring; an error is raised otherwise::
-
-        sage: set_random_seed()
-        sage: J1 = random_eja(field=AA)
-        sage: J2 = random_eja(field=QQ,orthonormalize=False)
-        sage: J = DirectSumEJA(J1,J2)
-        Traceback (most recent call last):
-        ...
-        ValueError: algebras must share the same base field
-
-    """
-    def __init__(self, J1, J2, **kwargs):
-        if J1.base_ring() != J2.base_ring():
-            raise ValueError("algebras must share the same base field")
-        field = J1.base_ring()
-
-        self._factors = (J1, J2)
-        n1 = J1.dimension()
-        n2 = J2.dimension()
-        n = n1+n2
-        V = VectorSpace(field, n)
-        mult_table = [ [ V.zero() for j in range(i+1) ]
-                       for i in range(n) ]
-        for i in range(n1):
-            for j in range(i+1):
-                p = (J1.monomial(i)*J1.monomial(j)).to_vector()
-                mult_table[i][j] = V(p.list() + [field.zero()]*n2)
-
-        for i in range(n2):
-            for j in range(i+1):
-                p = (J2.monomial(i)*J2.monomial(j)).to_vector()
-                mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
-
-        # TODO: build the IP table here from the two constituent IP
-        # matrices (it'll be block diagonal, I think).
-        ip_table = [ [ field.zero() for j in range(i+1) ]
-                       for i in range(n) ]
-        super(DirectSumEJA, self).__init__(field,
-                                           mult_table,
-                                           ip_table,
-                                           check_axioms=False,
-                                           **kwargs)
-        self.rank.set_cache(J1.rank() + J2.rank())
-
-
-    def factors(self):
-        r"""
-        Return the pair of this algebra's factors.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import (HadamardEJA,
-            ....:                                  JordanSpinEJA,
-            ....:                                  DirectSumEJA)
-
-        EXAMPLES::
-
-            sage: J1 = HadamardEJA(2, field=QQ)
-            sage: J2 = JordanSpinEJA(3, field=QQ)
-            sage: J = DirectSumEJA(J1,J2)
-            sage: J.factors()
-            (Euclidean Jordan algebra of dimension 2 over Rational Field,
-             Euclidean Jordan algebra of dimension 3 over Rational Field)
-
-        """
-        return self._factors
-
-    def projections(self):
-        r"""
-        Return a pair of projections onto this algebra's factors.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
-            ....:                                  ComplexHermitianEJA,
-            ....:                                  DirectSumEJA)
-
-        EXAMPLES::
-
-            sage: J1 = JordanSpinEJA(2)
-            sage: J2 = ComplexHermitianEJA(2)
-            sage: J = DirectSumEJA(J1,J2)
-            sage: (pi_left, pi_right) = J.projections()
-            sage: J.one().to_vector()
-            (1, 0, 1, 0, 0, 1)
-            sage: pi_left(J.one()).to_vector()
-            (1, 0)
-            sage: pi_right(J.one()).to_vector()
-            (1, 0, 0, 1)
-
-        """
-        (J1,J2) = self.factors()
-        m = J1.dimension()
-        n = J2.dimension()
-        V_basis = self.vector_space().basis()
-        # Need to specify the dimensions explicitly so that we don't
-        # wind up with a zero-by-zero matrix when we want e.g. a
-        # zero-by-two matrix (important for composing things).
-        P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
-        P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
-        pi_left = FiniteDimensionalEJAOperator(self,J1,P1)
-        pi_right = FiniteDimensionalEJAOperator(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 = FiniteDimensionalEJAOperator(J1,self,I1)
-        iota_right = FiniteDimensionalEJAOperator(J2,self,I2)
-        return (iota_left, iota_right)
-
-    def inner_product(self, x, y):
-        r"""
-        The standard Cartesian inner-product.
-
-        We project ``x`` and ``y`` onto our factors, and add up the
-        inner-products from the subalgebras.
-
-        SETUP::
-
-
-            sage: from mjo.eja.eja_algebra import (HadamardEJA,
-            ....:                                  QuaternionHermitianEJA,
-            ....:                                  DirectSumEJA)
-
-        EXAMPLE::
-
-            sage: J1 = HadamardEJA(3,field=QQ)
-            sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
-            sage: J = DirectSumEJA(J1,J2)
-            sage: x1 = J1.one()
-            sage: x2 = x1
-            sage: y1 = J2.one()
-            sage: y2 = y1
-            sage: x1.inner_product(x2)
-            3
-            sage: y1.inner_product(y2)
-            2
-            sage: J.one().inner_product(J.one())
-            5
-
-        """
-        (pi_left, pi_right) = self.projections()
-        x1 = pi_left(x)
-        x2 = pi_right(x)
-        y1 = pi_left(y)
-        y2 = pi_right(y)
-
-        return (x1.inner_product(y1) + x2.inner_product(y2))
+class DirectSumEJA(ConcreteEJA):
+    r"""
+    The external (orthogonal) direct sum of two other Euclidean Jordan
+    algebras. Essentially the Cartesian product of its two factors.
+    Every Euclidean Jordan algebra decomposes into an orthogonal
+    direct sum of simple Euclidean Jordan algebras, so no generality
+    is lost by providing only this construction.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import (random_eja,
+        ....:                                  HadamardEJA,
+        ....:                                  RealSymmetricEJA,
+        ....:                                  DirectSumEJA)
+
+    EXAMPLES::
+
+        sage: J1 = HadamardEJA(2)
+        sage: J2 = RealSymmetricEJA(3)
+        sage: J = DirectSumEJA(J1,J2)
+        sage: J.dimension()
+        8
+        sage: J.rank()
+        5
+
+    TESTS:
+
+    The external direct sum construction is only valid when the two factors
+    have the same base ring; an error is raised otherwise::
+
+        sage: set_random_seed()
+        sage: J1 = random_eja(field=AA)
+        sage: J2 = random_eja(field=QQ,orthonormalize=False)
+        sage: J = DirectSumEJA(J1,J2)
+        Traceback (most recent call last):
+        ...
+        ValueError: algebras must share the same base field
+
+    """
+    def __init__(self, J1, J2, **kwargs):
+        if J1.base_ring() != J2.base_ring():
+            raise ValueError("algebras must share the same base field")
+        field = J1.base_ring()
+
+        self._factors = (J1, J2)
+        n1 = J1.dimension()
+        n2 = J2.dimension()
+        n = n1+n2
+        V = VectorSpace(field, n)
+        mult_table = [ [ V.zero() for j in range(i+1) ]
+                       for i in range(n) ]
+        for i in range(n1):
+            for j in range(i+1):
+                p = (J1.monomial(i)*J1.monomial(j)).to_vector()
+                mult_table[i][j] = V(p.list() + [field.zero()]*n2)
+
+        for i in range(n2):
+            for j in range(i+1):
+                p = (J2.monomial(i)*J2.monomial(j)).to_vector()
+                mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
+
+        # TODO: build the IP table here from the two constituent IP
+        # matrices (it'll be block diagonal, I think).
+        ip_table = [ [ field.zero() for j in range(i+1) ]
+                       for i in range(n) ]
+        super(DirectSumEJA, self).__init__(field,
+                                           mult_table,
+                                           ip_table,
+                                           check_axioms=False,
+                                           **kwargs)
+        self.rank.set_cache(J1.rank() + J2.rank())
+
+
+    def factors(self):
+        r"""
+        Return the pair of this algebra's factors.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            ....:                                  JordanSpinEJA,
+            ....:                                  DirectSumEJA)
+
+        EXAMPLES::
+
+            sage: J1 = HadamardEJA(2, field=QQ)
+            sage: J2 = JordanSpinEJA(3, field=QQ)
+            sage: J = DirectSumEJA(J1,J2)
+            sage: J.factors()
+            (Euclidean Jordan algebra of dimension 2 over Rational Field,
+             Euclidean Jordan algebra of dimension 3 over Rational Field)
+
+        """
+        return self._factors
+
+    def projections(self):
+        r"""
+        Return a pair of projections onto this algebra's factors.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            ....:                                  ComplexHermitianEJA,
+            ....:                                  DirectSumEJA)
+
+        EXAMPLES::
+
+            sage: J1 = JordanSpinEJA(2)
+            sage: J2 = ComplexHermitianEJA(2)
+            sage: J = DirectSumEJA(J1,J2)
+            sage: (pi_left, pi_right) = J.projections()
+            sage: J.one().to_vector()
+            (1, 0, 1, 0, 0, 1)
+            sage: pi_left(J.one()).to_vector()
+            (1, 0)
+            sage: pi_right(J.one()).to_vector()
+            (1, 0, 0, 1)
+
+        """
+        (J1,J2) = self.factors()
+        m = J1.dimension()
+        n = J2.dimension()
+        V_basis = self.vector_space().basis()
+        # Need to specify the dimensions explicitly so that we don't
+        # wind up with a zero-by-zero matrix when we want e.g. a
+        # zero-by-two matrix (important for composing things).
+        P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
+        P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
+        pi_left = FiniteDimensionalEJAOperator(self,J1,P1)
+        pi_right = FiniteDimensionalEJAOperator(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 = FiniteDimensionalEJAOperator(J1,self,I1)
+        iota_right = FiniteDimensionalEJAOperator(J2,self,I2)
+        return (iota_left, iota_right)
+
+    def inner_product(self, x, y):
+        r"""
+        The standard Cartesian inner-product.
+
+        We project ``x`` and ``y`` onto our factors, and add up the
+        inner-products from the subalgebras.
+
+        SETUP::
+
+
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            ....:                                  QuaternionHermitianEJA,
+            ....:                                  DirectSumEJA)
+
+        EXAMPLE::
+
+            sage: J1 = HadamardEJA(3,field=QQ)
+            sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: J = DirectSumEJA(J1,J2)
+            sage: x1 = J1.one()
+            sage: x2 = x1
+            sage: y1 = J2.one()
+            sage: y2 = y1
+            sage: x1.inner_product(x2)
+            3
+            sage: y1.inner_product(y2)
+            2
+            sage: J.one().inner_product(J.one())
+            5
+
+        """
+        (pi_left, pi_right) = self.projections()
+        x1 = pi_left(x)
+        x2 = pi_right(x)
+        y1 = pi_left(y)
+        y2 = pi_right(y)
+
+        return (x1.inner_product(y1) + x2.inner_product(y2))