]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: drop a pointless "solve" in the EJA charpoly system.
[sage.d.git] / mjo / eja / eja_algebra.py
index ea62da8ca796031bd99382c5886bb1091d145b78..c35bf256453cfa9a3758c034a47a043745eb0c2b 100644 (file)
@@ -133,7 +133,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             deortho_vector_basis = tuple( V(b.list()) for b in basis )
 
             from mjo.eja.eja_utils import gram_schmidt
             deortho_vector_basis = tuple( V(b.list()) for b in basis )
 
             from mjo.eja.eja_utils import gram_schmidt
-            basis = gram_schmidt(basis, inner_product)
+            basis = tuple(gram_schmidt(basis, inner_product))
 
         # Save the (possibly orthonormalized) matrix basis for
         # later...
 
         # Save the (possibly orthonormalized) matrix basis for
         # later...
@@ -158,7 +158,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         # Now we actually compute the multiplication and inner-product
         # tables/matrices using the possibly-orthonormalized basis.
 
         # 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._inner_product_matrix = matrix.identity(field, n)
         self._multiplication_table = [ [0 for j in range(i+1)]
                                        for i in range(n) ]
 
         self._multiplication_table = [ [0 for j in range(i+1)]
                                        for i in range(n) ]
 
@@ -171,15 +171,20 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                 q_i = basis[i]
                 q_j = basis[j]
 
                 q_i = basis[i]
                 q_j = basis[j]
 
-                elt = jordan_product(q_i, q_j)
-                ip = inner_product(q_i, q_j)
-
                 # The jordan product returns a matrixy answer, so we
                 # have to convert it to the algebra coordinates.
                 # The jordan product returns a matrixy answer, so we
                 # have to convert it to the algebra coordinates.
+                elt = jordan_product(q_i, q_j)
                 elt = W.coordinate_vector(V(elt.list()))
                 self._multiplication_table[i][j] = self.from_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()
 
         self._inner_product_matrix._cache = {'hermitian': True}
         self._inner_product_matrix.set_immutable()
@@ -315,7 +320,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         This method should of course always return ``True``, unless
         this algebra was constructed with ``check_axioms=False`` and
 
         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
         """
 
         # Used to check whether or not something is zero in an inexact
@@ -732,7 +737,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         if self.is_trivial():
             return MatrixSpace(self.base_ring(), 0)
         else:
         if self.is_trivial():
             return MatrixSpace(self.base_ring(), 0)
         else:
-            return self._matrix_basis[0].matrix_space()
+            return self.matrix_basis()[0].parent()
 
 
     @cached_method
 
 
     @cached_method
@@ -745,23 +750,57 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: from mjo.eja.eja_algebra import (HadamardEJA,
             ....:                                  random_eja)
 
             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
 
 
             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:
 
         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: 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()
 
             sage: set_random_seed()
             sage: J = random_eja()
@@ -769,6 +808,27 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: expected = matrix.identity(J.base_ring(), J.dimension())
             sage: actual == expected
             True
             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::
 
         Ensure that the cached unit element (often precomputed by
         hand) agrees with the computed one::
@@ -780,6 +840,15 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: J.one() == cached
             True
 
             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.
         """
         # We can brute-force compute the matrices of the operators
         # that correspond to the basis elements of this algebra.
@@ -1032,6 +1101,21 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         r"""
         The `r` polynomial coefficients of the "characteristic polynomial
         of" function.
         r"""
         The `r` polynomial coefficients of the "characteristic polynomial
         of" function.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import random_eja
+
+        TESTS:
+
+        The theory shows that these are all homogeneous polynomials of
+        a known degree::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
+            True
+
         """
         n = self.dimension()
         R = self.coordinate_polynomial_ring()
         """
         n = self.dimension()
         R = self.coordinate_polynomial_ring()
@@ -1067,10 +1151,17 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         # The theory says that only the first "r" coefficients are
         # nonzero, and they actually live in the original polynomial
 
         # The theory says that only the first "r" coefficients are
         # nonzero, and they actually live in the original polynomial
-        # ring and not the fraction field. We negate them because
-        # in the actual characteristic polynomial, they get moved
-        # to the other side where x^r lives.
-        return -A_rref.solve_right(E*b).change_ring(R)[:r]
+        # ring and not the fraction field. We negate them because in
+        # the actual characteristic polynomial, they get moved to the
+        # other side where x^r lives. We don't bother to trim A_rref
+        # down to a square matrix and solve the resulting system,
+        # because the upper-left r-by-r portion of A_rref is
+        # guaranteed to be the identity matrix, so e.g.
+        #
+        #   A_rref.solve_right(Y)
+        #
+        # would just be returning Y.
+        return (-E*b)[:r].change_ring(R)
 
     @cached_method
     def rank(self):
 
     @cached_method
     def rank(self):
@@ -1131,7 +1222,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: set_random_seed()    # long time
             sage: J = random_eja()     # long time
 
             sage: set_random_seed()    # long time
             sage: J = random_eja()     # long time
-            sage: caches = J.rank()    # long time
+            sage: cached = J.rank()    # long time
             sage: J.rank.clear_cache() # long time
             sage: J.rank() == cached   # long time
             True
             sage: J.rank.clear_cache() # long time
             sage: J.rank() == cached   # long time
             True
@@ -1187,9 +1278,7 @@ class RationalBasisEJA(FiniteDimensionalEJA):
                  jordan_product,
                  inner_product,
                  field=AA,
                  jordan_product,
                  inner_product,
                  field=AA,
-                 orthonormalize=True,
                  check_field=True,
                  check_field=True,
-                 check_axioms=True,
                  **kwargs):
 
         if check_field:
                  **kwargs):
 
         if check_field:
@@ -1198,6 +1287,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")
 
             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.
         if field is not QQ:
             # There's no point in constructing the extra algebra if this
             # one is already rational.
@@ -1212,15 +1302,13 @@ class RationalBasisEJA(FiniteDimensionalEJA):
                                        field=QQ,
                                        orthonormalize=False,
                                        check_field=False,
                                        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,
 
         super().__init__(basis,
                          jordan_product,
                          inner_product,
                          field=field,
                          check_field=check_field,
-                         check_axioms=check_axioms,
                          **kwargs)
 
     @cached_method
                          **kwargs)
 
     @cached_method
@@ -1260,7 +1348,14 @@ class RationalBasisEJA(FiniteDimensionalEJA):
         a = ( a_i.change_ring(self.base_ring())
               for a_i in self._rational_algebra._charpoly_coefficients() )
 
         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
         # deorthonormalized ones.
         R = self.coordinate_polynomial_ring()
         from sage.modules.free_module_element import vector
@@ -1594,6 +1689,38 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
 
 
 class ComplexMatrixEJA(MatrixEJA):
 
 
 class ComplexMatrixEJA(MatrixEJA):
+    # A manual dictionary-cache for the complex_extension() method,
+    # since apparently @classmethods can't also be @cached_methods.
+    _complex_extension = {}
+
+    @classmethod
+    def complex_extension(cls,field):
+        r"""
+        The complex field that we embed/unembed, as an extension
+        of the given ``field``.
+        """
+        if field in cls._complex_extension:
+            return cls._complex_extension[field]
+
+        # Sage doesn't know how to adjoin the complex "i" (the root of
+        # x^2 + 1) to a field in a general way. Here, we just enumerate
+        # all of the cases that I have cared to support so far.
+        if field is AA:
+            # Sage doesn't know how to embed AA into QQbar, i.e. how
+            # to adjoin sqrt(-1) to AA.
+            F = QQbar
+        elif not field.is_exact():
+            # RDF or RR
+            F = field.complex_field()
+        else:
+            # Works for QQ and... maybe some other fields.
+            R = PolynomialRing(field, 'z')
+            z = R.gen()
+            F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
+
+        cls._complex_extension[field] = F
+        return F
+
     @staticmethod
     def dimension_over_reals():
         return 2
     @staticmethod
     def dimension_over_reals():
         return 2
@@ -1648,9 +1775,10 @@ class ComplexMatrixEJA(MatrixEJA):
 
         blocks = []
         for z in M.list():
 
         blocks = []
         for z in M.list():
-            a = z.list()[0] # real part, I guess
-            b = z.list()[1] # imag part, I guess
-            blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
+            a = z.real()
+            b = z.imag()
+            blocks.append(matrix(field, 2, [ [ a, b],
+                                             [-b, a] ]))
 
         return matrix.block(field, n, blocks)
 
 
         return matrix.block(field, n, blocks)
 
@@ -1689,26 +1817,7 @@ class ComplexMatrixEJA(MatrixEJA):
         super(ComplexMatrixEJA,cls).real_unembed(M)
         n = ZZ(M.nrows())
         d = cls.dimension_over_reals()
         super(ComplexMatrixEJA,cls).real_unembed(M)
         n = ZZ(M.nrows())
         d = cls.dimension_over_reals()
-
-        # If "M" was normalized, its base ring might have roots
-        # adjoined and they can stick around after unembedding.
-        field = M.base_ring()
-        R = PolynomialRing(field, 'z')
-        z = R.gen()
-
-        # Sage doesn't know how to adjoin the complex "i" (the root of
-        # x^2 + 1) to a field in a general way. Here, we just enumerate
-        # all of the cases that I have cared to support so far.
-        if field is AA:
-            # Sage doesn't know how to embed AA into QQbar, i.e. how
-            # to adjoin sqrt(-1) to AA.
-            F = QQbar
-        elif not field.is_exact():
-            # RDF or RR
-            F = field.complex_field()
-        else:
-            # Works for QQ and... maybe some other fields.
-            F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
+        F = cls.complex_extension(M.base_ring())
         i = F.gen()
 
         # Go top-left to bottom-right (reading order), converting every
         i = F.gen()
 
         # Go top-left to bottom-right (reading order), converting every
@@ -1804,7 +1913,6 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
-            sage: field = QuadraticField(2, 'sqrt2')
             sage: B = ComplexHermitianEJA._denormalized_basis(n)
             sage: all( M.is_symmetric() for M in  B)
             True
             sage: B = ComplexHermitianEJA._denormalized_basis(n)
             sage: all( M.is_symmetric() for M in  B)
             True
@@ -1822,18 +1930,27 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
         #   * The diagonal will (as a result) be real.
         #
         S = []
         #   * The diagonal will (as a result) be real.
         #
         S = []
+        Eij = matrix.zero(F,n)
         for i in range(n):
             for j in range(i+1):
         for i in range(n):
             for j in range(i+1):
-                Eij = matrix(F, n, lambda k,l: k==i and l==j)
+                # "build" E_ij
+                Eij[i,j] = 1
                 if i == j:
                     Sij = cls.real_embed(Eij)
                     S.append(Sij)
                 else:
                     # The second one has a minus because it's conjugated.
                 if i == j:
                     Sij = cls.real_embed(Eij)
                     S.append(Sij)
                 else:
                     # The second one has a minus because it's conjugated.
-                    Sij_real = cls.real_embed(Eij + Eij.transpose())
+                    Eij[j,i] = 1 # Eij = Eij + Eij.transpose()
+                    Sij_real = cls.real_embed(Eij)
                     S.append(Sij_real)
                     S.append(Sij_real)
-                    Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose())
+                    # Eij = I*Eij - I*Eij.transpose()
+                    Eij[i,j] = I
+                    Eij[j,i] = -I
+                    Sij_imag = cls.real_embed(Eij)
                     S.append(Sij_imag)
                     S.append(Sij_imag)
+                    Eij[j,i] = 0
+                # "erase" E_ij
+                Eij[i,j] = 0
 
         # Since we embedded these, we can drop back to the "field" that we
         # started with instead of the complex extension "F".
 
         # Since we embedded these, we can drop back to the "field" that we
         # started with instead of the complex extension "F".
@@ -1869,6 +1986,25 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
         return cls(n, **kwargs)
 
 class QuaternionMatrixEJA(MatrixEJA):
         return cls(n, **kwargs)
 
 class QuaternionMatrixEJA(MatrixEJA):
+
+    # A manual dictionary-cache for the quaternion_extension() method,
+    # since apparently @classmethods can't also be @cached_methods.
+    _quaternion_extension = {}
+
+    @classmethod
+    def quaternion_extension(cls,field):
+        r"""
+        The quaternion field that we embed/unembed, as an extension
+        of the given ``field``.
+        """
+        if field in cls._quaternion_extension:
+            return cls._quaternion_extension[field]
+
+        Q = QuaternionAlgebra(field,-1,-1)
+
+        cls._quaternion_extension[field] = Q
+        return Q
+
     @staticmethod
     def dimension_over_reals():
         return 4
     @staticmethod
     def dimension_over_reals():
         return 4
@@ -1973,8 +2109,7 @@ class QuaternionMatrixEJA(MatrixEJA):
 
         # Use the base ring of the matrix to ensure that its entries can be
         # multiplied by elements of the quaternion algebra.
 
         # Use the base ring of the matrix to ensure that its entries can be
         # multiplied by elements of the quaternion algebra.
-        field = M.base_ring()
-        Q = QuaternionAlgebra(field,-1,-1)
+        Q = cls.quaternion_extension(M.base_ring())
         i,j,k = Q.gens()
 
         # Go top-left to bottom-right (reading order), converting every
         i,j,k = Q.gens()
 
         # Go top-left to bottom-right (reading order), converting every
@@ -2089,23 +2224,39 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
         #   * The diagonal will (as a result) be real.
         #
         S = []
         #   * The diagonal will (as a result) be real.
         #
         S = []
+        Eij = matrix.zero(Q,n)
         for i in range(n):
             for j in range(i+1):
         for i in range(n):
             for j in range(i+1):
-                Eij = matrix(Q, n, lambda k,l: k==i and l==j)
+                # "build" E_ij
+                Eij[i,j] = 1
                 if i == j:
                     Sij = cls.real_embed(Eij)
                     S.append(Sij)
                 else:
                     # The second, third, and fourth ones have a minus
                     # because they're conjugated.
                 if i == j:
                     Sij = cls.real_embed(Eij)
                     S.append(Sij)
                 else:
                     # The second, third, and fourth ones have a minus
                     # because they're conjugated.
-                    Sij_real = cls.real_embed(Eij + Eij.transpose())
+                    # Eij = Eij + Eij.transpose()
+                    Eij[j,i] = 1
+                    Sij_real = cls.real_embed(Eij)
                     S.append(Sij_real)
                     S.append(Sij_real)
-                    Sij_I = cls.real_embed(I*Eij - I*Eij.transpose())
+                    # Eij = I*(Eij - Eij.transpose())
+                    Eij[i,j] = I
+                    Eij[j,i] = -I
+                    Sij_I = cls.real_embed(Eij)
                     S.append(Sij_I)
                     S.append(Sij_I)
-                    Sij_J = cls.real_embed(J*Eij - J*Eij.transpose())
+                    # Eij = J*(Eij - Eij.transpose())
+                    Eij[i,j] = J
+                    Eij[j,i] = -J
+                    Sij_J = cls.real_embed(Eij)
                     S.append(Sij_J)
                     S.append(Sij_J)
-                    Sij_K = cls.real_embed(K*Eij - K*Eij.transpose())
+                    # Eij = K*(Eij - Eij.transpose())
+                    Eij[i,j] = K
+                    Eij[j,i] = -K
+                    Sij_K = cls.real_embed(Eij)
                     S.append(Sij_K)
                     S.append(Sij_K)
+                    Eij[j,i] = 0
+                # "erase" E_ij
+                Eij[i,j] = 0
 
         # Since we embedded these, we can drop back to the "field" that we
         # started with instead of the quaternion algebra "Q".
 
         # Since we embedded these, we can drop back to the "field" that we
         # started with instead of the quaternion algebra "Q".
@@ -2311,31 +2462,38 @@ class BilinearFormEJA(ConcreteEJA):
         ....:              for j in range(n-1) ]
         sage: actual == expected
         True
         ....:              for j in range(n-1) ]
         sage: actual == expected
         True
+
     """
     def __init__(self, B, **kwargs):
     """
     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):
 
         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()
 
         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
             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()
 
         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)
                                               jordan_product,
                                               inner_product,
                                               **kwargs)
@@ -2516,7 +2674,8 @@ class TrivialEJA(ConcreteEJA):
         # inappropriate for us.
         return cls(**kwargs)
 
         # inappropriate for us.
         return cls(**kwargs)
 
-class DirectSumEJA(ConcreteEJA):
+
+class DirectSumEJA(FiniteDimensionalEJA):
     r"""
     The external (orthogonal) direct sum of two other Euclidean Jordan
     algebras. Essentially the Cartesian product of its two factors.
     r"""
     The external (orthogonal) direct sum of two other Euclidean Jordan
     algebras. Essentially the Cartesian product of its two factors.
@@ -2540,6 +2699,10 @@ class DirectSumEJA(ConcreteEJA):
         8
         sage: J.rank()
         5
         8
         sage: J.rank()
         5
+        sage: J.matrix_space()
+        The Cartesian product of (Full MatrixSpace of 2 by 1 dense matrices
+        over Algebraic Real Field, Full MatrixSpace of 3 by 3 dense matrices
+        over Algebraic Real Field)
 
     TESTS:
 
 
     TESTS:
 
@@ -2560,199 +2723,213 @@ class DirectSumEJA(ConcreteEJA):
             raise ValueError("algebras must share the same base field")
         field = J1.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::
+        M = J1.matrix_space().cartesian_product(J2.matrix_space())
+        self._cartprod_algebra = J1.cartesian_product(J2)
 
 
-            sage: from mjo.eja.eja_algebra import (HadamardEJA,
-            ....:                                  JordanSpinEJA,
-            ....:                                  DirectSumEJA)
+        self._matrix_basis = tuple( [M((a,0)) for a in J1.matrix_basis()] +
+                                    [M((0,b)) for b in J2.matrix_basis()] )
 
 
-        EXAMPLES::
+        n = len(self._matrix_basis)
+        self._sets = None
+        CombinatorialFreeModule.__init__(
+                         self,
+                         field,
+                         range(n),
+                         category=self._cartprod_algebra.category(),
+                         bracket=False,
+                         **kwargs)
+        self.rank.set_cache(J1.rank() + J2.rank())
 
 
-            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):
+    def product(self,x,y):
         r"""
         r"""
-        Return a pair of projections onto this algebra's factors.
-
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
             ....:                                  ComplexHermitianEJA,
             ....:                                  DirectSumEJA)
 
         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::
+        TESTS::
 
             sage: set_random_seed()
 
             sage: set_random_seed()
-            sage: J1 = random_eja()
-            sage: J2 = random_eja()
+            sage: J1 = JordanSpinEJA(3, field=QQ)
+            sage: J2 = ComplexHermitianEJA(2, field=QQ, orthonormalize=False)
             sage: J = DirectSumEJA(J1,J2)
             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()
+            sage: J.random_element()*J.random_element() in J
             True
 
         """
             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)
+        xv = self._cartprod_algebra.from_vector(x.to_vector())
+        yv = self._cartprod_algebra.from_vector(y.to_vector())
+        return self.from_vector((xv*yv).to_vector())
 
 
-    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.
+    def cartesian_factors(self):
+        r"""
+        Return the pair of this algebra's factors.
 
         SETUP::
 
 
         SETUP::
 
-
             sage: from mjo.eja.eja_algebra import (HadamardEJA,
             sage: from mjo.eja.eja_algebra import (HadamardEJA,
-            ....:                                  QuaternionHermitianEJA,
+            ....:                                  JordanSpinEJA,
             ....:                                  DirectSumEJA)
 
             ....:                                  DirectSumEJA)
 
-        EXAMPLE::
+        EXAMPLES::
 
 
-            sage: J1 = HadamardEJA(3,field=QQ)
-            sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: J1 = HadamardEJA(2, field=QQ)
+            sage: J2 = JordanSpinEJA(3, field=QQ)
             sage: J = DirectSumEJA(J1,J2)
             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
+            sage: J.cartesian_factors()
+            (Euclidean Jordan algebra of dimension 2 over Rational Field,
+             Euclidean Jordan algebra of dimension 3 over Rational Field)
 
         """
 
         """
-        (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))
+        return self._cartprod_algebra.cartesian_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))