]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: disable axiom checking for concrete algebras.
[sage.d.git] / mjo / eja / eja_algebra.py
index 14666c2f2a4d47b2f1a5cde1448dcb3b778aa72a..e075ed21ab287fcefcfbf7c9674a39281b48b69c 100644 (file)
@@ -119,10 +119,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         The ``field`` we're given must be real with ``check_field=True``::
 
-            sage: JordanSpinEJA(2,QQbar)
+            sage: JordanSpinEJA(2, field=QQbar)
             Traceback (most recent call last):
             ...
             ValueError: scalar field is not real
+            sage: JordanSpinEJA(2, field=QQbar, check_field=False)
+            Euclidean Jordan algebra of dimension 2 over Algebraic Field
 
         The multiplication table must be square with ``check_axioms=True``::
 
@@ -173,16 +175,27 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             if not all( len(l) == n for l in inner_product_table ):
                 raise ValueError(msg)
 
+            # Check commutativity of the Jordan product (symmetry of
+            # the multiplication table) and the commutativity of the
+            # inner-product (symmetry of the inner-product table)
+            # first if we're going to check them at all.. This has to
+            # be done before we define product_on_basis(), because
+            # that method assumes that self._multiplication_table is
+            # symmetric. And it has to be done before we build
+            # self._inner_product_matrix, because the process used to
+            # construct it assumes symmetry as well.
             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:
@@ -214,14 +227,18 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
                 elt = self.from_vector(multiplication_table[i][j])
                 self._multiplication_table[i][j] = 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': False}
+        ip_matrix_constructor = lambda i,j: inner_product_table[i][j] if j <= i else inner_product_table[j][i]
+        self._inner_product_matrix = matrix(field, n, ip_matrix_constructor)
+        self._inner_product_matrix._cache = {'hermitian': True}
+        self._inner_product_matrix.set_immutable()
 
         if check_axioms:
             if not self._is_jordanian():
@@ -276,6 +293,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: x = J.random_element()
             sage: J(x.to_vector().column()) == x
             True
+
         """
         msg = "not an element of this algebra"
         if elt == 0:
@@ -298,8 +316,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         # element's ring because the basis space might be an algebraic
         # closure whereas the base ring of the 3-by-3 identity matrix
         # could be QQ instead of QQbar.
+        #
+        # We pass check=False because the matrix basis is "guaranteed"
+        # to be linearly independent... right? Ha ha.
         V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols())
-        W = V.span_of_basis( _mat2vec(s) for s in self.matrix_basis() )
+        W = V.span_of_basis( (_mat2vec(s) for s in self.matrix_basis()),
+                             check=False)
 
         try:
             coords =  W.coordinate_vector(_mat2vec(elt))
@@ -378,7 +400,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         # Used to check whether or not something is zero in an inexact
         # ring. This number is sufficient to allow the construction of
-        # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
+        # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
         epsilon = 1e-16
 
         for i in range(self.dimension()):
@@ -475,7 +497,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,field=QQ,orthonormalize=False)
             sage: J.coordinate_polynomial_ring()
             Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
 
@@ -1117,7 +1139,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
     Element = FiniteDimensionalEuclideanJordanAlgebraElement
 
-class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlgebra):
+class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
     r"""
     New class for algebras whose supplied basis elements have all rational entries.
 
@@ -1140,16 +1162,26 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge
 
     """
     def __init__(self,
-                 field,
                  basis,
                  jordan_product,
                  inner_product,
+                 field=AA,
                  orthonormalize=True,
                  prefix='e',
                  category=None,
                  check_field=True,
                  check_axioms=True):
 
+        if check_field:
+            # Abuse the check_field parameter to check that the entries of
+            # out basis (in ambient coordinates) are in the field QQ.
+            if not all( all(b_i in QQ for b_i in b.list()) for b in basis ):
+                raise TypeError("basis not rational")
+
+        # Temporary(?) hack to ensure that the matrix and vector bases
+        # are over the same ring.
+        basis = tuple( b.change_ring(field) for b in basis )
+
         n = len(basis)
         vector_basis = basis
 
@@ -1160,6 +1192,7 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge
         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:
@@ -1173,27 +1206,97 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge
         # 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
+        self._rational_algebra = None
+
+        if orthonormalize:
+            if self.base_ring() is not QQ:
+                # There's no point in constructing the extra algebra if this
+                # one is already rational. If the original basis is rational
+                # but normalization would make it irrational, then this whole
+                # constructor will just fail anyway as it tries to stick an
+                # irrational number into a rational algebra.
+                #
+                # Note: the same Jordan and inner-products work here,
+                # because they are necessarily defined with respect to
+                # ambient coordinates and not any particular basis.
+                self._rational_algebra = RationalBasisEuclideanJordanAlgebra(
+                                          basis,
+                                          jordan_product,
+                                          inner_product,
+                                          field=QQ,
+                                          orthonormalize=False,
+                                          prefix=prefix,
+                                          category=category,
+                                          check_field=False,
+                                          check_axioms=False)
+
+            # Compute the deorthonormalized tables before we orthonormalize
+            # the given basis. The "check" parameter here guarantees that
+            # the basis is linearly-independent.
+            W = V.span_of_basis( vector_basis, check=check_axioms)
+
+            # 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)
+
+        # 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
-            vector_basis = gram_schmidt(vector_basis, inner_product)
-            W = V.span_of_basis( vector_basis )
+            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
 
             if basis_is_matrices:
-                from mjo.eja.eja_utils import _vec2mat
                 basis = tuple( map(_vec2mat,basis) )
 
-        W = V.span_of_basis( vector_basis )
+        W = V.span_of_basis( vector_basis, check=check_axioms)
 
-        mult_table = [ [0 for i in range(n)] for j in range(n) ]
-        ip_table = [ [0 for i in range(n)] for j in range(n) ]
+        # 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, check=check_axioms)
+        self._deortho_matrix = matrix( U.coordinate_vector(q)
+                                       for q in vector_basis )
 
-        # Note: the Jordan and inner- products are defined in terms
+        # If the superclass constructor is going to verify the
+        # symmetry of this table, it has better at least be
+        # square...
+        if check_axioms:
+            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) ]
+        else:
+            mult_table = [ [0 for j in range(i+1)] for i in range(n) ]
+            ip_table = [ [0 for j in range(i+1)] 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):
@@ -1216,9 +1319,12 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge
 
                 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 check_axioms:
+                    # The tables are square if we're verifying that they
+                    # are commutative.
+                    mult_table[j][i] = elt
+                    ip_table[j][i] = ip
 
         if basis_is_matrices:
             for m in basis:
@@ -1235,26 +1341,13 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge
                          check_field,
                          check_axioms)
 
-class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
-    r"""
-    Algebras whose basis consists of vectors with rational
-    entries. Equivalently, algebras whose multiplication tables
-    contain only rational coefficients.
-
-    When an EJA has a basis that can be made rational, we can speed up
-    the computation of its characteristic polynomial by doing it over
-    ``QQ``. All of the named EJA constructors that we provide fall
-    into this category.
-    """
     @cached_method
     def _charpoly_coefficients(self):
         r"""
-        Override the parent method with something that tries to compute
-        over a faster (non-extension) field.
-
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import JordanSpinEJA
+            sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
+            ....:                                  JordanSpinEJA)
 
         EXAMPLES:
 
@@ -1272,29 +1365,31 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr
             Algebraic Real Field
 
         """
-        if self.base_ring() is QQ:
+        if self.base_ring() is QQ or self._rational_algebra is None:
             # There's no need to construct *another* algebra over the
-            # rationals if this one is already over the rationals.
+            # rationals if this one is already over the
+            # rationals. Likewise, if we never orthonormalized our
+            # basis, we might as well just use the given one.
             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.
-        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))
+        # the same, because all we've done is a change of basis.
+        # Then, change back from QQ to our real base ring
+        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
+        # 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
 
-class ConcreteEuclideanJordanAlgebra:
+        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(RationalBasisEuclideanJordanAlgebra):
     r"""
     A class for the Euclidean Jordan algebras that we know by name.
 
@@ -1347,7 +1442,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.
 
@@ -1355,127 +1450,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 = [[field.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
-
-        super(MatrixEuclideanJordanAlgebra, self).__init__(field,
-                                                           mult_table,
-                                                           ip_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):
         """
@@ -1500,8 +1482,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()
@@ -1534,8 +1520,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
@@ -1558,9 +1544,9 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra,
 
     In theory, our "field" can be any subfield of the reals::
 
-        sage: RealSymmetricEJA(2, RDF)
+        sage: RealSymmetricEJA(2, field=RDF)
         Euclidean Jordan algebra of dimension 3 over Real Double Field
-        sage: RealSymmetricEJA(2, RR)
+        sage: RealSymmetricEJA(2, field=RR)
         Euclidean Jordan algebra of dimension 3 over Real Field with
         53 bits of precision
 
@@ -1601,7 +1587,7 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra,
 
     """
     @classmethod
-    def _denormalized_basis(cls, n, field):
+    def _denormalized_basis(cls, n):
         """
         Return a basis for the space of real symmetric n-by-n matrices.
 
@@ -1613,7 +1599,7 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra,
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
-            sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
+            sage: B = RealSymmetricEJA._denormalized_basis(n)
             sage: all( M.is_symmetric() for M in  B)
             True
 
@@ -1623,13 +1609,13 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra,
         S = []
         for i in range(n):
             for j in range(i+1):
-                Eij = matrix(field, n, lambda k,l: k==i and l==j)
+                Eij = matrix(ZZ, n, lambda k,l: k==i and l==j)
                 if i == j:
                     Sij = Eij
                 else:
                     Sij = Eij + Eij.transpose()
                 S.append(Sij)
-        return S
+        return tuple(S)
 
 
     @staticmethod
@@ -1637,20 +1623,24 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra,
         return 4 # Dimension 10
 
     @classmethod
-    def random_instance(cls, field=AA, **kwargs):
+    def random_instance(cls, **kwargs):
         """
         Return a random instance of this type of algebra.
         """
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
-        return cls(n, field, **kwargs)
+        return cls(n, **kwargs)
+
+    def __init__(self, n, **kwargs):
+        # 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
 
-    def __init__(self, n, field=AA, **kwargs):
-        basis = self._denormalized_basis(n, field)
-        super(RealSymmetricEJA, self).__init__(field,
-                                               basis,
-                                               check_axioms=False,
+        super(RealSymmetricEJA, self).__init__(self._denormalized_basis(n),
+                                               self.jordan_product,
+                                               self.trace_inner_product,
                                                **kwargs)
         self.rank.set_cache(n)
+        self.one.set_cache(self(matrix.identity(ZZ,n)))
 
 
 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
@@ -1781,7 +1771,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.
@@ -1803,16 +1793,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,
@@ -1827,9 +1817,9 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra,
 
     In theory, our "field" can be any subfield of the reals::
 
-        sage: ComplexHermitianEJA(2, RDF)
+        sage: ComplexHermitianEJA(2, field=RDF)
         Euclidean Jordan algebra of dimension 4 over Real Double Field
-        sage: ComplexHermitianEJA(2, RR)
+        sage: ComplexHermitianEJA(2, field=RR)
         Euclidean Jordan algebra of dimension 4 over Real Field with
         53 bits of precision
 
@@ -1871,7 +1861,7 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra,
     """
 
     @classmethod
-    def _denormalized_basis(cls, n, field):
+    def _denormalized_basis(cls, n):
         """
         Returns a basis for the space of complex Hermitian n-by-n matrices.
 
@@ -1890,15 +1880,16 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra,
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
             sage: field = QuadraticField(2, 'sqrt2')
-            sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
+            sage: B = ComplexHermitianEJA._denormalized_basis(n)
             sage: all( M.is_symmetric() for M in  B)
             True
 
         """
+        field = ZZ
         R = PolynomialRing(field, 'z')
         z = R.gen()
         F = field.extension(z**2 + 1, 'I')
-        I = F.gen()
+        I = F.gen(1)
 
         # This is like the symmetric case, but we need to be careful:
         #
@@ -1921,28 +1912,32 @@ 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)
+    def __init__(self, n, **kwargs):
+        # 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
+
+        super(ComplexHermitianEJA, self).__init__(self._denormalized_basis(n),
+                                                  self.jordan_product,
+                                                  self.trace_inner_product,
+                                                  **kwargs)
         self.rank.set_cache(n)
+        # TODO: pre-cache the identity!
 
     @staticmethod
     def _max_random_instance_size():
         return 3 # Dimension 9
 
     @classmethod
-    def random_instance(cls, field=AA, **kwargs):
+    def random_instance(cls, **kwargs):
         """
         Return a random instance of this type of algebra.
         """
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
-        return cls(n, field, **kwargs)
+        return cls(n, **kwargs)
 
 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
     @staticmethod
@@ -2076,7 +2071,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.
@@ -2098,16 +2093,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
@@ -2122,9 +2117,9 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
 
     In theory, our "field" can be any subfield of the reals::
 
-        sage: QuaternionHermitianEJA(2, RDF)
+        sage: QuaternionHermitianEJA(2, field=RDF)
         Euclidean Jordan algebra of dimension 6 over Real Double Field
-        sage: QuaternionHermitianEJA(2, RR)
+        sage: QuaternionHermitianEJA(2, field=RR)
         Euclidean Jordan algebra of dimension 6 over Real Field with
         53 bits of precision
 
@@ -2165,7 +2160,7 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
 
     """
     @classmethod
-    def _denormalized_basis(cls, n, field):
+    def _denormalized_basis(cls, n):
         """
         Returns a basis for the space of quaternion Hermitian n-by-n matrices.
 
@@ -2183,11 +2178,12 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
 
             sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
-            sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
+            sage: B = QuaternionHermitianEJA._denormalized_basis(n)
             sage: all( M.is_symmetric() for M in B )
             True
 
         """
+        field = ZZ
         Q = QuaternionAlgebra(QQ,-1,-1)
         I,J,K = Q.gens()
 
@@ -2217,16 +2213,20 @@ 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, **kwargs):
+        # 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
 
-    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__(self._denormalized_basis(n),
+                                                     self.jordan_product,
+                                                     self.trace_inner_product,
+                                                     **kwargs)
         self.rank.set_cache(n)
+        # TODO: cache one()!
 
     @staticmethod
     def _max_random_instance_size():
@@ -2236,16 +2236,15 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
         return 2 # Dimension 6
 
     @classmethod
-    def random_instance(cls, field=AA, **kwargs):
+    def random_instance(cls, **kwargs):
         """
         Return a random instance of this type of algebra.
         """
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
-        return cls(n, field, **kwargs)
+        return cls(n, **kwargs)
 
 
-class HadamardEJA(RationalBasisEuclideanJordanAlgebraNg,
-                  ConcreteEuclideanJordanAlgebra):
+class HadamardEJA(ConcreteEuclideanJordanAlgebra):
     """
     Return the Euclidean Jordan Algebra corresponding to the set
     `R^n` under the Hadamard product.
@@ -2285,17 +2284,25 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebraNg,
         (r0, r1, r2)
 
     """
-    def __init__(self, n, field=AA, **kwargs):
-        V = VectorSpace(field, n)
-        basis = V.basis()
-
+    def __init__(self, n, **kwargs):
         def jordan_product(x,y):
-            return V([ xi*yi for (xi,yi) in zip(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)
 
-        super(HadamardEJA, self).__init__(field,
-                                          basis,
+        # New defaults for keyword arguments. Don't orthonormalize
+        # because our basis is already orthonormal with respect to our
+        # inner-product. Don't check the axioms, because we know this
+        # is a valid EJA... but do double-check if the user passes
+        # check_axioms=True. Note: we DON'T override the "check_field"
+        # default here, because the user can pass in a field!
+        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)
@@ -2314,16 +2321,15 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebraNg,
         return 5
 
     @classmethod
-    def random_instance(cls, field=AA, **kwargs):
+    def random_instance(cls, **kwargs):
         """
         Return a random instance of this type of algebra.
         """
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
-        return cls(n, field, **kwargs)
+        return cls(n, **kwargs)
 
 
-class BilinearFormEJA(RationalBasisEuclideanJordanAlgebraNg,
-                      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 =
@@ -2404,27 +2410,30 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebraNg,
         sage: actual == expected
         True
     """
-    def __init__(self, B, field=AA, **kwargs):
+    def __init__(self, B, **kwargs):
         if not B.is_positive_definite():
             raise ValueError("bilinear form is not positive-definite")
 
-        n = B.nrows()
-        V = VectorSpace(field, n)
-
         def inner_product(x,y):
             return (B*x).inner_product(y)
 
         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)
             zbar = y0*xbar + x0*ybar
-            return V([z0] + zbar.list())
+            return P((z0,) + tuple(zbar))
 
-        super(BilinearFormEJA, self).__init__(field,
-                                              V.basis(),
+        # 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
+
+        n = B.nrows()
+        standard_basis = FreeModule(ZZ, n).basis()
+        super(BilinearFormEJA, self).__init__(standard_basis,
                                               jordan_product,
                                               inner_product,
                                               **kwargs)
@@ -2447,28 +2456,28 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebraNg,
         return 5
 
     @classmethod
-    def random_instance(cls, field=AA, **kwargs):
+    def random_instance(cls, **kwargs):
         """
         Return a random instance of this algebra.
         """
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
         if n.is_zero():
-            B = matrix.identity(field, n)
-            return cls(B, field, **kwargs)
+            B = matrix.identity(ZZ, n)
+            return cls(B, **kwargs)
 
-        B11 = matrix.identity(field,1)
-        M = matrix.random(field, n-1)
-        I = matrix.identity(field, n-1)
-        alpha = field.zero()
+        B11 = matrix.identity(ZZ, 1)
+        M = matrix.random(ZZ, n-1)
+        I = matrix.identity(ZZ, n-1)
+        alpha = ZZ.zero()
         while alpha.is_zero():
-            alpha = field.random_element().abs()
+            alpha = ZZ.random_element().abs()
         B22 = M.transpose()*M + alpha*I
 
         from sage.matrix.special import block_matrix
         B = block_matrix(2,2, [ [B11,   ZZ(0) ],
                                 [ZZ(0), B22 ] ])
 
-        return cls(B, field, **kwargs)
+        return cls(B, **kwargs)
 
 
 class JordanSpinEJA(BilinearFormEJA):
@@ -2521,11 +2530,18 @@ class JordanSpinEJA(BilinearFormEJA):
             True
 
     """
-    def __init__(self, n, field=AA, **kwargs):
-        # This is a special case of the BilinearFormEJA with the identity
-        # matrix as its bilinear form.
-        B = matrix.identity(field, n)
-        super(JordanSpinEJA, self).__init__(B, field, **kwargs)
+    def __init__(self, n, **kwargs):
+        # This is a special case of the BilinearFormEJA with the
+        # identity matrix as its bilinear form.
+        B = matrix.identity(ZZ, n)
+
+        # Don't orthonormalize because our basis is already
+        # orthonormal with respect to our inner-product.
+        if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
+
+        # But also don't pass check_field=False here, because the user
+        # can pass in a field!
+        super(JordanSpinEJA, self).__init__(B, **kwargs)
 
     @staticmethod
     def _max_random_instance_size():
@@ -2535,18 +2551,17 @@ class JordanSpinEJA(BilinearFormEJA):
         return 5
 
     @classmethod
-    def random_instance(cls, field=AA, **kwargs):
+    def random_instance(cls, **kwargs):
         """
         Return a random instance of this type of algebra.
 
         Needed here to override the implementation for ``BilinearFormEJA``.
         """
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
-        return cls(n, field, **kwargs)
+        return cls(n, **kwargs)
 
 
-class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra,
-                 ConcreteEuclideanJordanAlgebra):
+class TrivialEJA(ConcreteEuclideanJordanAlgebra):
     """
     The trivial Euclidean Jordan algebra consisting of only a zero element.
 
@@ -2575,13 +2590,18 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra,
         0
 
     """
-    def __init__(self, field=AA, **kwargs):
-        mult_table = []
-        ip_table = []
-        super(TrivialEJA, self).__init__(field,
-                                         mult_table,
-                                         ip_table,
-                                         check_axioms=False,
+    def __init__(self, **kwargs):
+        jordan_product = lambda x,y: x
+        inner_product = lambda x,y: 0
+        basis = ()
+
+        # New defaults for keyword arguments
+        if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
+        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
+
+        super(TrivialEJA, self).__init__(basis,
+                                         jordan_product,
+                                         inner_product,
                                          **kwargs)
         # The rank is zero using my definition, namely the dimension of the
         # largest subalgebra generated by any element.
@@ -2589,10 +2609,10 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra,
         self.one.set_cache( self.zero() )
 
     @classmethod
-    def random_instance(cls, field=AA, **kwargs):
+    def random_instance(cls, **kwargs):
         # We don't take a "size" argument so the superclass method is
         # inappropriate for us.
-        return cls(field, **kwargs)
+        return cls(**kwargs)
 
 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
     r"""
@@ -2625,8 +2645,8 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
     have the same base ring; an error is raised otherwise::
 
         sage: set_random_seed()
-        sage: J1 = random_eja(AA)
-        sage: J2 = random_eja(QQ)
+        sage: J1 = random_eja(field=AA)
+        sage: J2 = random_eja(field=QQ,orthonormalize=False)
         sage: J = DirectSumEJA(J1,J2)
         Traceback (most recent call last):
         ...
@@ -2643,21 +2663,22 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
         n2 = J2.dimension()
         n = n1+n2
         V = VectorSpace(field, n)
-        mult_table = [ [ V.zero() for j in range(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(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(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 = None
+        ip_table = [ [ field.zero() for j in range(i+1) ]
+                       for i in range(n) ]
         super(DirectSumEJA, self).__init__(field,
                                            mult_table,
                                            ip_table,
@@ -2678,8 +2699,8 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
 
         EXAMPLES::
 
-            sage: J1 = HadamardEJA(2,QQ)
-            sage: J2 = JordanSpinEJA(3,QQ)
+            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,
@@ -2808,8 +2829,8 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
 
         EXAMPLE::
 
-            sage: J1 = HadamardEJA(3,QQ)
-            sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
+            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