]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: convert matrix algebas to the new constructor, fix all tests.
[sage.d.git] / mjo / eja / eja_algebra.py
index adcf75503d2a55add4c32dc56380396abd1f8d9d..efa10e2a856f90cbae752a0fcb19eefe9a637ded 100644 (file)
@@ -34,7 +34,9 @@ from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
 from mjo.eja.eja_utils import _mat2vec
 
 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
-
+    r"""
+    The lowest-level class for representing a Euclidean Jordan algebra.
+    """
     def _coerce_map_from_base_ring(self):
         """
         Disable the map from the base ring into the algebra.
@@ -61,13 +63,23 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
     def __init__(self,
                  field,
-                 mult_table,
+                 multiplication_table,
+                 inner_product_table,
                  prefix='e',
                  category=None,
                  matrix_basis=None,
                  check_field=True,
                  check_axioms=True):
         """
+        INPUT:
+
+          * field -- the scalar field for this algebra (must be real)
+
+          * multiplication_table -- the multiplication table for this
+            algebra's implicit basis. Only the lower-triangular portion
+            of the table is used, since the multiplication is assumed
+            to be commutative.
+
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (
@@ -85,6 +97,24 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: x*y == y*x
             True
 
+        An error is raised if the Jordan product is not commutative::
+
+            sage: JP = ((1,2),(0,0))
+            sage: IP = ((1,0),(0,1))
+            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
+            Traceback (most recent call last):
+            ...
+            ValueError: Jordan product is not commutative
+
+        An error is raised if the inner-product is not commutative::
+
+            sage: JP = ((1,0),(0,1))
+            sage: IP = ((1,2),(0,0))
+            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
+            Traceback (most recent call last):
+            ...
+            ValueError: inner-product is not commutative
+
         TESTS:
 
         The ``field`` we're given must be real with ``check_field=True``::
@@ -96,11 +126,26 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         The multiplication table must be square with ``check_axioms=True``::
 
-            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()))
+            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()),((1,),))
             Traceback (most recent call last):
             ...
             ValueError: multiplication table is not square
 
+        The multiplication and inner-product tables must be the same
+        size (and in particular, the inner-product table must also be
+        square) with ``check_axioms=True``::
+
+            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),(()))
+            Traceback (most recent call last):
+            ...
+            ValueError: multiplication and inner-product tables are
+            different sizes
+            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),((1,2),))
+            Traceback (most recent call last):
+            ...
+            ValueError: multiplication and inner-product tables are
+            different sizes
+
         """
         if check_field:
             if not field.is_subring(RR):
@@ -109,12 +154,35 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
                 # we've specified a real embedding.
                 raise ValueError("scalar field is not real")
 
-        # The multiplication table had better be square
-        n = len(mult_table)
+
+        # The multiplication and inner-product tables should be square
+        # if the user wants us to verify them. And we verify them as
+        # soon as possible, because we want to exploit their symmetry.
+        n = len(multiplication_table)
         if check_axioms:
-            if not all( len(l) == n for l in mult_table ):
+            if not all( len(l) == n for l in multiplication_table ):
                 raise ValueError("multiplication table is not square")
 
+            # If the multiplication table is square, we can check if
+            # the inner-product table is square by comparing it to the
+            # multiplication table's dimensions.
+            msg = "multiplication and inner-product tables are different sizes"
+            if not len(inner_product_table) == n:
+                raise ValueError(msg)
+
+            if not all( len(l) == n for l in inner_product_table ):
+                raise ValueError(msg)
+
+            if not all(    multiplication_table[j][i]
+                        == multiplication_table[i][j]
+                        for i in range(n)
+                        for j in range(i+1) ):
+                raise ValueError("Jordan product is not commutative")
+            if not all( inner_product_table[j][i]
+                        == inner_product_table[i][j]
+                        for i in range(n)
+                        for j in range(i+1) ):
+                raise ValueError("inner-product is not commutative")
         self._matrix_basis = matrix_basis
 
         if category is None:
@@ -134,19 +202,31 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         # long run to have the multiplication table be in terms of
         # algebra elements. We do this after calling the superclass
         # constructor so that from_vector() knows what to do.
+        #
+        # Note: we take advantage of symmetry here, and only store
+        # the lower-triangular portion of the table.
         self._multiplication_table = [ [ self.vector_space().zero()
-                                         for i in range(n) ]
-                                       for j in range(n) ]
-        # take advantage of symmetry
+                                         for j in range(i+1) ]
+                                       for i in range(n) ]
+
         for i in range(n):
             for j in range(i+1):
-                elt = self.from_vector(mult_table[i][j])
+                elt = self.from_vector(multiplication_table[i][j])
                 self._multiplication_table[i][j] = elt
-                self._multiplication_table[j][i] = elt
+
+        self._multiplication_table = tuple(map(tuple, self._multiplication_table))
+
+        # Save our inner product as a matrix, since the efficiency of
+        # matrix multiplication will usually outweigh the fact that we
+        # have to store a redundant upper- or lower-triangular part.
+        # Pre-cache the fact that these are Hermitian (real symmetric,
+        # in fact) in case some e.g. matrix multiplication routine can
+        # take advantage of it.
+        self._inner_product_matrix = matrix(field, inner_product_table)
+        self._inner_product_matrix._cache = {'hermitian': True}
+        self._inner_product_matrix.set_immutable()
 
         if check_axioms:
-            if not self._is_commutative():
-                raise ValueError("algebra is not commutative")
             if not self._is_jordanian():
                 raise ValueError("Jordan identity does not hold")
             if not self._inner_product_is_associative():
@@ -253,7 +333,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return fmt.format(self.dimension(), self.base_ring())
 
     def product_on_basis(self, i, j):
-        return self._multiplication_table[i][j]
+        # We only stored the lower-triangular portion of the
+        # multiplication table.
+        if j <= i:
+            return self._multiplication_table[i][j]
+        else:
+            return self._multiplication_table[j][i]
 
     def _is_commutative(self):
         r"""
@@ -393,7 +478,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: J = HadamardEJA(2)
             sage: J.coordinate_polynomial_ring()
             Multivariate Polynomial Ring in X1, X2...
-            sage: J = RealSymmetricEJA(3,QQ)
+            sage: J = RealSymmetricEJA(3,QQ,orthonormalize=False)
             sage: J.coordinate_polynomial_ring()
             Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
 
@@ -509,10 +594,20 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             +----++----+----+----+----+
 
         """
-        M = list(self._multiplication_table) # copy
-        for i in range(len(M)):
-            # M had better be "square"
+        n = self.dimension()
+        M = [ [ self.zero() for j in range(n) ]
+              for i in range(n) ]
+        for i in range(n):
+            for j in range(i+1):
+                M[i][j] = self._multiplication_table[i][j]
+                M[j][i] = M[i][j]
+
+        for i in range(n):
+            # Prepend the left "header" column entry Can't do this in
+            # the loop because it messes up the symmetry.
             M[i] = [self.monomial(i)] + M[i]
+
+        # Prepend the header row.
         M = [["*"] + list(self.gens())] + M
         return table(M, header_row=True, header_column=True, frame=True)
 
@@ -1025,27 +1120,196 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
     Element = FiniteDimensionalEuclideanJordanAlgebraElement
 
-
 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
     r"""
-    Algebras whose basis consists of vectors with rational
-    entries. Equivalently, algebras whose multiplication tables
-    contain only rational coefficients.
-
-    When an EJA has a basis that can be made rational, we can speed up
-    the computation of its characteristic polynomial by doing it over
-    ``QQ``. All of the named EJA constructors that we provide fall
-    into this category.
+    New class for algebras whose supplied basis elements have all rational entries.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import BilinearFormEJA
+
+    EXAMPLES:
+
+    The supplied basis is orthonormalized by default::
+
+        sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
+        sage: J = BilinearFormEJA(B)
+        sage: J.matrix_basis()
+        (
+        [1]  [  0]  [   0]
+        [0]  [1/5]  [32/5]
+        [0], [  0], [   5]
+        )
+
     """
+    def __init__(self,
+                 field,
+                 basis,
+                 jordan_product,
+                 inner_product,
+                 orthonormalize=True,
+                 prefix='e',
+                 category=None,
+                 check_field=True,
+                 check_axioms=True):
+
+        n = len(basis)
+        vector_basis = basis
+
+        from sage.structure.element import is_Matrix
+        basis_is_matrices = False
+
+        degree = 0
+        if n > 0:
+            if is_Matrix(basis[0]):
+                basis_is_matrices = True
+                from mjo.eja.eja_utils import _vec2mat
+                vector_basis = tuple( map(_mat2vec,basis) )
+                degree = basis[0].nrows()**2
+            else:
+                degree = basis[0].degree()
+
+        V = VectorSpace(field, degree)
+
+        # If we were asked to orthonormalize, and if the orthonormal
+        # basis is different from the given one, then we also want to
+        # compute multiplication and inner-product tables for the
+        # deorthonormalized basis. These can be used later to
+        # construct a deorthonormalized copy of this algebra over QQ
+        # in which several operations are much faster.
+        self._deortho_multiplication_table = None
+        self._deortho_inner_product_table = None
+
+        if orthonormalize:
+            # Compute the deorthonormalized tables before we orthonormalize
+            # the given basis.
+            W = V.span_of_basis( vector_basis )
+
+            # TODO: use symmetry
+            self._deortho_multiplication_table = [ [0 for j in range(n)]
+                                                   for i in range(n) ]
+            self._deortho_inner_product_table = [ [0 for j in range(n)]
+                                                  for i in range(n) ]
+
+            # 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):
+                    # given 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)
+
+                    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: use symmetry
+                    elt = W.coordinate_vector(elt)
+                    self._deortho_multiplication_table[i][j] = elt
+                    self._deortho_multiplication_table[j][i] = elt
+                    self._deortho_inner_product_table[i][j] = ip
+                    self._deortho_inner_product_table[j][i] = ip
+
+        if self._deortho_multiplication_table is not None:
+            self._deortho_multiplication_table = tuple(map(tuple, self._deortho_multiplication_table))
+        if self._deortho_inner_product_table is not None:
+            self._deortho_inner_product_table = tuple(map(tuple, self._deortho_inner_product_table))
+
+        # 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
+        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)
+
+            W = V.span_of_basis( vector_basis )
+
+            # Normalize the "matrix" basis, too!
+            basis = vector_basis
+
+            if basis_is_matrices:
+                basis = tuple( map(_vec2mat,basis) )
+
+        W = V.span_of_basis( vector_basis )
+
+        # Now "W" is the vector space of our algebra coordinates. The
+        # variables "X1", "X2",...  refer to the entries of vectors in
+        # W. Thus to convert back and forth between the orthonormal
+        # coordinates and the given ones, we need to stick the original
+        # basis in W.
+        U = V.span_of_basis( deortho_vector_basis )
+        self._deortho_matrix = matrix( U.coordinate_vector(q)
+                                       for q in vector_basis )
+
+        # TODO: use symmetry
+        mult_table = [ [0 for j in range(n)] for i in range(n) ]
+        ip_table = [ [0 for j in range(n)] for i in range(n) ]
+
+        # 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)
+
+                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: use symmetry
+                elt = W.coordinate_vector(elt)
+                mult_table[i][j] = elt
+                mult_table[j][i] = elt
+                ip_table[i][j] = ip
+                ip_table[j][i] = ip
+
+        if basis_is_matrices:
+            for m in basis:
+                m.set_immutable()
+        else:
+            basis = tuple( x.column() for x in basis )
+
+        super().__init__(field,
+                         mult_table,
+                         ip_table,
+                         prefix,
+                         category,
+                         basis, # matrix basis
+                         check_field,
+                         check_axioms)
+
     @cached_method
     def _charpoly_coefficients(self):
         r"""
-        Override the parent method with something that tries to compute
-        over a faster (non-extension) field.
-
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import JordanSpinEJA
+            sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
+            ....:                                  JordanSpinEJA)
 
         EXAMPLES:
 
@@ -1069,23 +1333,27 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr
             superclass = super(RationalBasisEuclideanJordanAlgebra, self)
             return superclass._charpoly_coefficients()
 
-        mult_table = tuple(
-            tuple(map(lambda x: x.to_vector(), ls))
-            for ls in self._multiplication_table
-        )
-
         # Do the computation over the rationals. The answer will be
-        # the same, because our basis coordinates are (essentially)
-        # rational.
+        # the same, because all we've done is a change of basis.
         J = FiniteDimensionalEuclideanJordanAlgebra(QQ,
-                                                    mult_table,
-                                                    check_field=False,
-                                                    check_axioms=False)
-        a = J._charpoly_coefficients()
-        return tuple(map(lambda x: x.change_ring(self.base_ring()), a))
+                                                    self._deortho_multiplication_table,
+                                                    self._deortho_inner_product_table)
+
+        # Change back from QQ to our real base ring
+        a = ( a_i.change_ring(self.base_ring())
+              for a_i in J._charpoly_coefficients() )
 
+        # Now convert the coordinate variables back to the
+        # deorthonormalized ones.
+        R = self.coordinate_polynomial_ring()
+        from sage.modules.free_module_element import vector
+        X = vector(R, R.gens())
+        BX = self._deortho_matrix*X
+
+        subs_dict = { X[i]: BX[i] for i in range(len(X)) }
+        return tuple( a_i.subs(subs_dict) for a_i in a )
 
-class ConcreteEuclideanJordanAlgebra:
+class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra):
     r"""
     A class for the Euclidean Jordan algebras that we know by name.
 
@@ -1138,7 +1406,7 @@ class ConcreteEuclideanJordanAlgebra:
         raise NotImplementedError
 
     @classmethod
-    def random_instance(cls, field=AA, **kwargs):
+    def random_instance(cls, *args, **kwargs):
         """
         Return a random instance of this type of algebra.
 
@@ -1146,132 +1414,14 @@ class ConcreteEuclideanJordanAlgebra:
         """
         from sage.misc.prandom import choice
         eja_class = choice(cls.__subclasses__())
-        return eja_class.random_instance(field)
-
-
-class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
-
-    def __init__(self, field, basis, normalize_basis=True, **kwargs):
-        """
-        Compared to the superclass constructor, we take a basis instead of
-        a multiplication table because the latter can be computed in terms
-        of the former when the product is known (like it is here).
-        """
-        # Used in this class's fast _charpoly_coefficients() override.
-        self._basis_normalizers = None
-
-        # We're going to loop through this a few times, so now's a good
-        # time to ensure that it isn't a generator expression.
-        basis = tuple(basis)
-
-        algebra_dim = len(basis)
-        degree = 0 # size of the matrices
-        if algebra_dim > 0:
-            degree = basis[0].nrows()
-
-        if algebra_dim > 1 and normalize_basis:
-            # We'll need sqrt(2) to normalize the basis, and this
-            # winds up in the multiplication table, so the whole
-            # algebra needs to be over the field extension.
-            R = PolynomialRing(field, 'z')
-            z = R.gen()
-            p = z**2 - 2
-            if p.is_irreducible():
-                field = field.extension(p, 'sqrt2', embedding=RLF(2).sqrt())
-                basis = tuple( s.change_ring(field) for s in basis )
-            self._basis_normalizers = tuple(
-                ~(self.matrix_inner_product(s,s).sqrt()) for s in basis )
-            basis = tuple(s*c for (s,c) in zip(basis,self._basis_normalizers))
-
-        # Now compute the multiplication and inner product tables.
-        # We have to do this *after* normalizing the basis, because
-        # scaling affects the answers.
-        V = VectorSpace(field, degree**2)
-        W = V.span_of_basis( _mat2vec(s) for s in basis )
-        mult_table = [[W.zero() for j in range(algebra_dim)]
-                                for i in range(algebra_dim)]
-        ip_table = [[W.zero() for j in range(algebra_dim)]
-                              for i in range(algebra_dim)]
-        for i in range(algebra_dim):
-            for j in range(algebra_dim):
-                mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
-                mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
-
-                try:
-                    # HACK: ignore the error here if we don't need the
-                    # inner product (as is the case when we construct
-                    # a dummy QQ-algebra for fast charpoly coefficients.
-                    ip_table[i][j] = self.matrix_inner_product(basis[i],
-                                                                basis[j])
-                except:
-                    pass
-
-        try:
-            # HACK PART DEUX
-            self._inner_product_matrix = matrix(field,ip_table)
-        except:
-            pass
-
-        super(MatrixEuclideanJordanAlgebra, self).__init__(field,
-                                                           mult_table,
-                                                           matrix_basis=basis,
-                                                           **kwargs)
-
-        if algebra_dim == 0:
-            self.one.set_cache(self.zero())
-        else:
-            n = basis[0].nrows()
-            # The identity wrt (A,B) -> (AB + BA)/2 is independent of the
-            # details of this algebra.
-            self.one.set_cache(self(matrix.identity(field,n)))
-
 
-    @cached_method
-    def _charpoly_coefficients(self):
-        r"""
-        Override the parent method with something that tries to compute
-        over a faster (non-extension) field.
-        """
-        if self._basis_normalizers is None or self.base_ring() is QQ:
-            # We didn't normalize, or the basis we started with had
-            # entries in a nice field already. Just compute the thing.
-            return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coefficients()
-
-        basis = ( (b/n) for (b,n) in zip(self.matrix_basis(),
-                                         self._basis_normalizers) )
-
-        # Do this over the rationals and convert back at the end.
-        # Only works because we know the entries of the basis are
-        # integers. The argument ``check_axioms=False`` is required
-        # because the trace inner-product method for this
-        # class is a stub and can't actually be checked.
-        J = MatrixEuclideanJordanAlgebra(QQ,
-                                         basis,
-                                         normalize_basis=False,
-                                         check_field=False,
-                                         check_axioms=False)
-        a = J._charpoly_coefficients()
-
-        # Unfortunately, changing the basis does change the
-        # coefficients of the characteristic polynomial, but since
-        # these are really the coefficients of the "characteristic
-        # polynomial of" function, everything is still nice and
-        # unevaluated. It's therefore "obvious" how scaling the
-        # basis affects the coordinate variables X1, X2, et
-        # cetera. Scaling the first basis vector up by "n" adds a
-        # factor of 1/n into every "X1" term, for example. So here
-        # we simply undo the basis_normalizer scaling that we
-        # performed earlier.
-        #
-        # The a[0] access here is safe because trivial algebras
-        # won't have any basis normalizers and therefore won't
-        # make it to this "else" branch.
-        XS = a[0].parent().gens()
-        subs_dict = { XS[i]: self._basis_normalizers[i]*XS[i]
-                      for i in range(len(XS)) }
-        return tuple( a_i.subs(subs_dict) for a_i in a )
+        # These all bubble up to the RationalBasisEuclideanJordanAlgebra
+        # superclass constructor, so any (kw)args valid there are also
+        # valid here.
+        return eja_class.random_instance(*args, **kwargs)
 
 
+class MatrixEuclideanJordanAlgebra:
     @staticmethod
     def real_embed(M):
         """
@@ -1296,8 +1446,12 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         """
         raise NotImplementedError
 
+    @staticmethod
+    def jordan_product(X,Y):
+        return (X*Y + Y*X)/2
+
     @classmethod
-    def matrix_inner_product(cls,X,Y):
+    def trace_inner_product(cls,X,Y):
         Xu = cls.real_unembed(X)
         Yu = cls.real_unembed(Y)
         tr = (Xu*Yu).trace()
@@ -1330,8 +1484,8 @@ class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         return M
 
 
-class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra,
-                       ConcreteEuclideanJordanAlgebra):
+class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra,
+                       RealMatrixEuclideanJordanAlgebra):
     """
     The rank-n simple EJA consisting of real symmetric n-by-n
     matrices, the usual symmetric Jordan product, and the trace inner
@@ -1444,9 +1598,11 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra,
         basis = self._denormalized_basis(n, field)
         super(RealSymmetricEJA, self).__init__(field,
                                                basis,
-                                               check_axioms=False,
+                                               self.jordan_product,
+                                               self.trace_inner_product,
                                                **kwargs)
         self.rank.set_cache(n)
+        self.one.set_cache(self(matrix.identity(field,n)))
 
 
 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
@@ -1577,7 +1733,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
 
 
     @classmethod
-    def matrix_inner_product(cls,X,Y):
+    def trace_inner_product(cls,X,Y):
         """
         Compute a matrix inner product in this algebra directly from
         its real embedding.
@@ -1599,16 +1755,16 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             sage: X = ComplexHermitianEJA.real_unembed(Xe)
             sage: Y = ComplexHermitianEJA.real_unembed(Ye)
             sage: expected = (X*Y).trace().real()
-            sage: actual = ComplexHermitianEJA.matrix_inner_product(Xe,Ye)
+            sage: actual = ComplexHermitianEJA.trace_inner_product(Xe,Ye)
             sage: actual == expected
             True
 
         """
-        return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/2
+        return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/2
 
 
-class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra,
-                          ConcreteEuclideanJordanAlgebra):
+class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra,
+                          ComplexMatrixEuclideanJordanAlgebra):
     """
     The rank-n simple EJA consisting of complex Hermitian n-by-n
     matrices over the real numbers, the usual symmetric Jordan product,
@@ -1717,16 +1873,18 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra,
 
         # Since we embedded these, we can drop back to the "field" that we
         # started with instead of the complex extension "F".
-        return ( s.change_ring(field) for s in S )
+        return tuple( s.change_ring(field) for s in S )
 
 
     def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n,field)
-        super(ComplexHermitianEJA,self).__init__(field,
-                                                 basis,
-                                                 check_axioms=False,
-                                                 **kwargs)
+        super(ComplexHermitianEJA, self).__init__(field,
+                                               basis,
+                                               self.jordan_product,
+                                               self.trace_inner_product,
+                                               **kwargs)
         self.rank.set_cache(n)
+        # TODO: pre-cache the identity!
 
     @staticmethod
     def _max_random_instance_size():
@@ -1872,7 +2030,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
 
 
     @classmethod
-    def matrix_inner_product(cls,X,Y):
+    def trace_inner_product(cls,X,Y):
         """
         Compute a matrix inner product in this algebra directly from
         its real embedding.
@@ -1894,16 +2052,16 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             sage: X = QuaternionHermitianEJA.real_unembed(Xe)
             sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
             sage: expected = (X*Y).trace().coefficient_tuple()[0]
-            sage: actual = QuaternionHermitianEJA.matrix_inner_product(Xe,Ye)
+            sage: actual = QuaternionHermitianEJA.trace_inner_product(Xe,Ye)
             sage: actual == expected
             True
 
         """
-        return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/4
+        return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/4
 
 
-class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
-                             ConcreteEuclideanJordanAlgebra):
+class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra,
+                             QuaternionMatrixEuclideanJordanAlgebra):
     r"""
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the
@@ -2013,16 +2171,18 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
 
         # Since we embedded these, we can drop back to the "field" that we
         # started with instead of the quaternion algebra "Q".
-        return ( s.change_ring(field) for s in S )
+        return tuple( s.change_ring(field) for s in S )
 
 
     def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n,field)
-        super(QuaternionHermitianEJA,self).__init__(field,
-                                                    basis,
-                                                    check_axioms=False,
-                                                    **kwargs)
+        super(QuaternionHermitianEJA, self).__init__(field,
+                                                     basis,
+                                                     self.jordan_product,
+                                                     self.trace_inner_product,
+                                                     **kwargs)
         self.rank.set_cache(n)
+        # TODO: cache one()!
 
     @staticmethod
     def _max_random_instance_size():
@@ -2040,8 +2200,7 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
         return cls(n, field, **kwargs)
 
 
-class HadamardEJA(RationalBasisEuclideanJordanAlgebra,
-                  ConcreteEuclideanJordanAlgebra):
+class HadamardEJA(ConcreteEuclideanJordanAlgebra):
     """
     Return the Euclidean Jordan Algebra corresponding to the set
     `R^n` under the Hadamard product.
@@ -2083,22 +2242,17 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra,
     """
     def __init__(self, n, field=AA, **kwargs):
         V = VectorSpace(field, n)
-        mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
-                       for i in range(n) ]
+        basis = V.basis()
 
-        # Inner products are real numbers and not algebra
-        # elements, so once we turn the algebra element
-        # into a vector in inner_product(), we never go
-        # back. As a result -- contrary to what we do with
-        # self._multiplication_table -- we store the inner
-        # product table as a plain old matrix and not as
-        # an algebra operator.
-        ip_table = matrix.identity(field,n)
-        self._inner_product_matrix = ip_table
+        def jordan_product(x,y):
+            return V([ xi*yi for (xi,yi) in zip(x,y) ])
+        def inner_product(x,y):
+            return x.inner_product(y)
 
         super(HadamardEJA, self).__init__(field,
-                                          mult_table,
-                                          check_axioms=False,
+                                          basis,
+                                          jordan_product,
+                                          inner_product,
                                           **kwargs)
         self.rank.set_cache(n)
 
@@ -2123,8 +2277,7 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra,
         return cls(n, field, **kwargs)
 
 
-class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
-                      ConcreteEuclideanJordanAlgebra):
+class BilinearFormEJA(ConcreteEuclideanJordanAlgebra):
     r"""
     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
     with the half-trace inner product and jordan product ``x*y =
@@ -2180,7 +2333,9 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
     We can check the multiplication condition given in the Jordan, von
     Neumann, and Wigner paper (and also discussed on my "On the
     symmetry..." paper). Note that this relies heavily on the standard
-    choice of basis, as does anything utilizing the bilinear form matrix::
+    choice of basis, as does anything utilizing the bilinear form
+    matrix.  We opt not to orthonormalize the basis, because if we
+    did, we would have to normalize the `s_{i}` in a similar manner::
 
         sage: set_random_seed()
         sage: n = ZZ.random_element(5)
@@ -2189,10 +2344,10 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
         sage: B22 = M.transpose()*M
         sage: B = block_matrix(2,2,[ [B11,0  ],
         ....:                        [0, B22 ] ])
-        sage: J = BilinearFormEJA(B)
+        sage: J = BilinearFormEJA(B, orthonormalize=False)
         sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
         sage: V = J.vector_space()
-        sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
+        sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
         ....:         for ei in eis ]
         sage: actual = [ sis[i]*sis[j]
         ....:            for i in range(n-1)
@@ -2204,39 +2359,28 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
         True
     """
     def __init__(self, B, field=AA, **kwargs):
-        n = B.nrows()
-
         if not B.is_positive_definite():
             raise ValueError("bilinear form is not positive-definite")
 
+        n = B.nrows()
         V = VectorSpace(field, n)
-        mult_table = [[V.zero() for j in range(n)] for i in range(n)]
-        for i in range(n):
-            for j in range(n):
-                x = V.gen(i)
-                y = V.gen(j)
-                x0 = x[0]
-                xbar = x[1:]
-                y0 = y[0]
-                ybar = y[1:]
-                z0 = (B*x).inner_product(y)
-                zbar = y0*xbar + x0*ybar
-                z = V([z0] + zbar.list())
-                mult_table[i][j] = z
-
-        # Inner products are real numbers and not algebra
-        # elements, so once we turn the algebra element
-        # into a vector in inner_product(), we never go
-        # back. As a result -- contrary to what we do with
-        # self._multiplication_table -- we store the inner
-        # product table as a plain old matrix and not as
-        # an algebra operator.
-        ip_table = B
-        self._inner_product_matrix = ip_table
+
+        def inner_product(x,y):
+            return (B*x).inner_product(y)
+
+        def jordan_product(x,y):
+            x0 = x[0]
+            xbar = x[1:]
+            y0 = y[0]
+            ybar = y[1:]
+            z0 = inner_product(x,y)
+            zbar = y0*xbar + x0*ybar
+            return V([z0] + zbar.list())
 
         super(BilinearFormEJA, self).__init__(field,
-                                              mult_table,
-                                              check_axioms=False,
+                                              V.basis(),
+                                              jordan_product,
+                                              inner_product,
                                               **kwargs)
 
         # The rank of this algebra is two, unless we're in a
@@ -2355,8 +2499,7 @@ class JordanSpinEJA(BilinearFormEJA):
         return cls(n, field, **kwargs)
 
 
-class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra,
-                 ConcreteEuclideanJordanAlgebra):
+class TrivialEJA(ConcreteEuclideanJordanAlgebra):
     """
     The trivial Euclidean Jordan algebra consisting of only a zero element.
 
@@ -2386,11 +2529,13 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra,
 
     """
     def __init__(self, field=AA, **kwargs):
-        mult_table = []
-        self._inner_product_matrix = matrix(field,0)
+        jordan_product = lambda x,y: x
+        inner_product = lambda x,y: field(0)
+        basis = ()
         super(TrivialEJA, self).__init__(field,
-                                         mult_table,
-                                         check_axioms=False,
+                                         basis,
+                                         jordan_product,
+                                         inner_product,
                                          **kwargs)
         # The rank is zero using my definition, namely the dimension of the
         # largest subalgebra generated by any element.
@@ -2464,8 +2609,12 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
                 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 = None
         super(DirectSumEJA, self).__init__(field,
                                            mult_table,
+                                           ip_table,
                                            check_axioms=False,
                                            **kwargs)
         self.rank.set_cache(J1.rank() + J2.rank())
@@ -2614,7 +2763,7 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
         EXAMPLE::
 
             sage: J1 = HadamardEJA(3,QQ)
-            sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
+            sage: J2 = QuaternionHermitianEJA(2,QQ,orthonormalize=False)
             sage: J = DirectSumEJA(J1,J2)
             sage: x1 = J1.one()
             sage: x2 = x1