eja: drop a pointless "solve" in the EJA charpoly system.

[sage.d.git] / mjo / eja / eja_algebra.py

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index 1250fbda7981aba5474af0a3d2615c969bc9d1e1..c35bf256453cfa9a3758c034a47a043745eb0c2b 100644 (file)
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -113,60 +113,36 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # we see in things like x = 1*e1 + 2*e2.
         vector_basis = basis
 
         # 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:
         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)
 
         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:
         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)
-
-            # Normalize the "matrix" basis, too!
-            basis = vector_basis
+            # 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 )

 
 

-            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:
         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.
 
         # 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
         # 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)
                 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._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()


@@ -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
 
         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


@@ -768,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


@@ -781,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
+
+        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: 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()


@@ -805,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::


@@ -816,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.


@@ -1068,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()


@@ -1103,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):


@@ -1167,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


@@ -1223,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:


@@ -1234,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.


@@ -1248,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


@@ -1296,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


@@ -1630,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


@@ -1684,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)
 


@@ -1725,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


@@ -1840,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


@@ -1858,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".


@@ -1905,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


@@ -2009,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


@@ -2125,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".


@@ -2222,11 +2337,16 @@ class HadamardEJA(ConcreteEJA):
 
     """
     def __init__(self, n, **kwargs):
 
     """
     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
 
         # New defaults for keyword arguments. Don't orthonormalize
         # because our basis is already orthonormal with respect to our


@@ -2237,12 +2357,8 @@ class HadamardEJA(ConcreteEJA):
         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
         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:
         self.rank.set_cache(n)
 
         if n == 0:


@@ -2346,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)


@@ -2551,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.


@@ -2575,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:
 


@@ -2595,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))