]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: reimplement multiplication_table() without member variables.
[sage.d.git] / mjo / eja / eja_algebra.py
index 4609ca2a4df7a4ff762e3eef00435a70b85e8cbb..3b13ce10bbf864ceb667f3ba4fc1f8b520327292 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,22 +97,57 @@ 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``::
 
-            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``::
 
-            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 +156,46 @@ 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)
+
+            # 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:
@@ -134,19 +215,32 @@ 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.
+        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_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():
@@ -199,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:
@@ -221,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))
@@ -253,7 +352,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"""
@@ -296,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()):
@@ -393,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...
 
@@ -509,11 +613,16 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             +----++----+----+----+----+
 
         """
-        M = list(self._multiplication_table) # copy
-        for i in range(len(M)):
-            # M had better be "square"
-            M[i] = [self.monomial(i)] + M[i]
-        M = [["*"] + list(self.gens())] + M
+        n = self.dimension()
+        # Prepend the header row.
+        M = [["*"] + list(self.gens())]
+
+        # And to each subsequent row, prepend an entry that belongs to
+        # the left-side "header column."
+        M += [ [self.monomial(i)] + [ self.product_on_basis(i,j)
+                                      for j in range(n) ]
+               for i in range(n) ]
+
         return table(M, header_row=True, header_column=True, frame=True)
 
 
@@ -1025,18 +1134,49 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
     Element = FiniteDimensionalEuclideanJordanAlgebraElement
 
-class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlgebra):
+class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
+    r"""
+    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,
+                 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
 
@@ -1047,6 +1187,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:
@@ -1054,38 +1195,124 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge
 
         V = VectorSpace(field, degree)
 
-        # Compute this from "Q" (obtained from Gram-Schmidt) below as
-        # R = Q.solve_right(A), where the rows of "Q" are the
-        # orthonormalized vector_basis and and the rows of "A" are the
-        # original vector_basis.
-        self._deorthonormalization_matrix = None
+        # Save a copy of an algebra with the original, rational basis
+        # and over QQ where computations are fast.
+        self._rational_algebra = None
+
+        if field is not QQ:
+            # There's no point in constructing the extra algebra if this
+            # one is already rational.
+            #
+            # 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)
+
+        if orthonormalize:
+            # 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:
-                from mjo.eja.eja_utils import _vec2mat
-                basis = tuple( map(_vec2mat,vector_basis) )
+                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
 
-        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) ]
+            if basis_is_matrices:
+                basis = tuple( map(_vec2mat,basis) )
+
+        W = V.span_of_basis( vector_basis, check=check_axioms)
+
+        # 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 )
+
+        # 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):
             for j in range(i+1):
-                # do another mat2vec because the multiplication
-                # table is in terms of vectors
-                elt = _mat2vec(jordan_product(basis[i],basis[j]))
+                # 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)
+
                 elt = W.coordinate_vector(elt)
                 mult_table[i][j] = elt
-                mult_table[j][i] = elt
-                ip = inner_product(basis[i],basis[j])
                 ip_table[i][j] = ip
-                ip_table[j][i] = ip
-
-        self._inner_product_matrix = matrix(field,ip_table)
+                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:
@@ -1095,32 +1322,20 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge
 
         super().__init__(field,
                          mult_table,
+                         ip_table,
                          prefix,
                          category,
                          basis, # matrix basis
                          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:
 
@@ -1138,29 +1353,31 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr
             Algebraic Real Field
 
         """
-        if self.base_ring() is QQ:
+        if 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
 
+        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.
 
@@ -1213,7 +1430,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.
 
@@ -1221,134 +1438,31 @@ 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)))
+        # 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)
 
 
-    @cached_method
-    def _charpoly_coefficients(self):
+class MatrixEuclideanJordanAlgebra:
+    @staticmethod
+    def dimension_over_reals():
         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 )
+        The dimension of this matrix's base ring over the reals.
 
+        The reals are dimension one over themselves, obviously; that's
+        just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
+        have dimension two. Finally, the quaternions have dimension
+        four over the reals.
 
-    @staticmethod
-    def real_embed(M):
+        This is used to determine the size of the matrix returned from
+        :meth:`real_embed`, among other things.
+        """
+        raise NotImplementedError
+
+    @classmethod
+    def real_embed(cls,M):
         """
         Embed the matrix ``M`` into a space of real matrices.
 
@@ -1361,18 +1475,83 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
           real_embed(M*N) = real_embed(M)*real_embed(N)
 
         """
-        raise NotImplementedError
+        if M.ncols() != M.nrows():
+            raise ValueError("the matrix 'M' must be square")
+        return M
 
 
-    @staticmethod
-    def real_unembed(M):
+    @classmethod
+    def real_unembed(cls,M):
         """
         The inverse of :meth:`real_embed`.
         """
-        raise NotImplementedError
+        if M.ncols() != M.nrows():
+            raise ValueError("the matrix 'M' must be square")
+        if not ZZ(M.nrows()).mod(cls.dimension_over_reals()).is_zero():
+            raise ValueError("the matrix 'M' must be a real embedding")
+        return M
+
+    @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):
+        r"""
+        Compute the trace inner-product of two real-embeddings.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
+            ....:                                  ComplexHermitianEJA,
+            ....:                                  QuaternionHermitianEJA)
+
+        EXAMPLES::
+
+        This gives the same answer as it would if we computed the trace
+        from the unembedded (original) matrices::
+
+            sage: set_random_seed()
+            sage: J = RealSymmetricEJA.random_instance()
+            sage: x,y = J.random_elements(2)
+            sage: Xe = x.to_matrix()
+            sage: Ye = y.to_matrix()
+            sage: X = J.real_unembed(Xe)
+            sage: Y = J.real_unembed(Ye)
+            sage: expected = (X*Y).trace()
+            sage: actual = J.trace_inner_product(Xe,Ye)
+            sage: actual == expected
+            True
+
+        ::
+
+            sage: set_random_seed()
+            sage: J = ComplexHermitianEJA.random_instance()
+            sage: x,y = J.random_elements(2)
+            sage: Xe = x.to_matrix()
+            sage: Ye = y.to_matrix()
+            sage: X = J.real_unembed(Xe)
+            sage: Y = J.real_unembed(Ye)
+            sage: expected = (X*Y).trace().real()
+            sage: actual = J.trace_inner_product(Xe,Ye)
+            sage: actual == expected
+            True
+
+        ::
+
+            sage: set_random_seed()
+            sage: J = QuaternionHermitianEJA.random_instance()
+            sage: x,y = J.random_elements(2)
+            sage: Xe = x.to_matrix()
+            sage: Ye = y.to_matrix()
+            sage: X = J.real_unembed(Xe)
+            sage: Y = J.real_unembed(Ye)
+            sage: expected = (X*Y).trace().coefficient_tuple()[0]
+            sage: actual = J.trace_inner_product(Xe,Ye)
+            sage: actual == expected
+            True
+
+        """
         Xu = cls.real_unembed(X)
         Yu = cls.real_unembed(Y)
         tr = (Xu*Yu).trace()
@@ -1389,24 +1568,12 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
 
 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
     @staticmethod
-    def real_embed(M):
-        """
-        The identity function, for embedding real matrices into real
-        matrices.
-        """
-        return M
-
-    @staticmethod
-    def real_unembed(M):
-        """
-        The identity function, for unembedding real matrices from real
-        matrices.
-        """
-        return M
+    def dimension_over_reals():
+        return 1
 
 
-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
@@ -1429,9 +1596,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
 
@@ -1472,7 +1639,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.
 
@@ -1484,7 +1651,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
 
@@ -1494,13 +1661,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
@@ -1508,25 +1675,39 @@ 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, field=AA, **kwargs):
-        basis = self._denormalized_basis(n, field)
-        super(RealSymmetricEJA, self).__init__(field,
-                                               basis,
-                                               check_axioms=False,
+    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(RealSymmetricEJA, self).__init__(self._denormalized_basis(n),
+                                               self.jordan_product,
+                                               self.trace_inner_product,
                                                **kwargs)
+
+        # TODO: this could be factored out somehow, but is left here
+        # because the MatrixEuclideanJordanAlgebra is not presently
+        # a subclass of the FDEJA class that defines rank() and one().
         self.rank.set_cache(n)
+        idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
+        self.one.set_cache(self(idV))
+
 
 
 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
     @staticmethod
-    def real_embed(M):
+    def dimension_over_reals():
+        return 2
+
+    @classmethod
+    def real_embed(cls,M):
         """
         Embed the n-by-n complex matrix ``M`` into the space of real
         matrices of size 2n-by-2n via the map the sends each entry `z = a +
@@ -1568,9 +1749,8 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             True
 
         """
+        super(ComplexMatrixEuclideanJordanAlgebra,cls).real_embed(M)
         n = M.nrows()
-        if M.ncols() != n:
-            raise ValueError("the matrix 'M' must be square")
 
         # We don't need any adjoined elements...
         field = M.base_ring().base_ring()
@@ -1584,8 +1764,8 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         return matrix.block(field, n, blocks)
 
 
-    @staticmethod
-    def real_unembed(M):
+    @classmethod
+    def real_unembed(cls,M):
         """
         The inverse of _embed_complex_matrix().
 
@@ -1616,31 +1796,37 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             True
 
         """
+        super(ComplexMatrixEuclideanJordanAlgebra,cls).real_unembed(M)
         n = ZZ(M.nrows())
-        if M.ncols() != n:
-            raise ValueError("the matrix 'M' must be square")
-        if not n.mod(2).is_zero():
-            raise ValueError("the matrix 'M' must be a complex embedding")
+        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())
         i = F.gen()
 
         # Go top-left to bottom-right (reading order), converting every
         # 2-by-2 block we see to a single complex element.
         elements = []
-        for k in range(n/2):
-            for j in range(n/2):
-                submat = M[2*k:2*k+2,2*j:2*j+2]
+        for k in range(n/d):
+            for j in range(n/d):
+                submat = M[d*k:d*k+d,d*j:d*j+d]
                 if submat[0,0] != submat[1,1]:
                     raise ValueError('bad on-diagonal submatrix')
                 if submat[0,1] != -submat[1,0]:
@@ -1648,42 +1834,11 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
                 z = submat[0,0] + submat[0,1]*i
                 elements.append(z)
 
-        return matrix(F, n/2, elements)
-
-
-    @classmethod
-    def matrix_inner_product(cls,X,Y):
-        """
-        Compute a matrix inner product in this algebra directly from
-        its real embedding.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
-
-        TESTS:
+        return matrix(F, n/d, elements)
 
-        This gives the same answer as the slow, default method implemented
-        in :class:`MatrixEuclideanJordanAlgebra`::
 
-            sage: set_random_seed()
-            sage: J = ComplexHermitianEJA.random_instance()
-            sage: x,y = J.random_elements(2)
-            sage: Xe = x.to_matrix()
-            sage: Ye = y.to_matrix()
-            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 == expected
-            True
-
-        """
-        return RealMatrixEuclideanJordanAlgebra.matrix_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,
@@ -1698,9 +1853,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
 
@@ -1742,7 +1897,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.
 
@@ -1761,15 +1916,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:
         #
@@ -1792,32 +1948,44 @@ 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, **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(ComplexHermitianEJA,self).__init__(field,
-                                                 basis,
-                                                 check_axioms=False,
-                                                 **kwargs)
+        super(ComplexHermitianEJA, self).__init__(self._denormalized_basis(n),
+                                                  self.jordan_product,
+                                                  self.trace_inner_product,
+                                                  **kwargs)
+        # TODO: this could be factored out somehow, but is left here
+        # because the MatrixEuclideanJordanAlgebra is not presently
+        # a subclass of the FDEJA class that defines rank() and one().
         self.rank.set_cache(n)
+        idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
+        self.one.set_cache(self(idV))
 
     @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
-    def real_embed(M):
+    def dimension_over_reals():
+        return 4
+
+    @classmethod
+    def real_embed(cls,M):
         """
         Embed the n-by-n quaternion matrix ``M`` into the space of real
         matrices of size 4n-by-4n by first sending each quaternion entry `z
@@ -1856,10 +2024,9 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             True
 
         """
+        super(QuaternionMatrixEuclideanJordanAlgebra,cls).real_embed(M)
         quaternions = M.base_ring()
         n = M.nrows()
-        if M.ncols() != n:
-            raise ValueError("the matrix 'M' must be square")
 
         F = QuadraticField(-1, 'I')
         i = F.gen()
@@ -1882,8 +2049,8 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
 
 
 
-    @staticmethod
-    def real_unembed(M):
+    @classmethod
+    def real_unembed(cls,M):
         """
         The inverse of _embed_quaternion_matrix().
 
@@ -1913,11 +2080,9 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             True
 
         """
+        super(QuaternionMatrixEuclideanJordanAlgebra,cls).real_unembed(M)
         n = ZZ(M.nrows())
-        if M.ncols() != n:
-            raise ValueError("the matrix 'M' must be square")
-        if not n.mod(4).is_zero():
-            raise ValueError("the matrix 'M' must be a quaternion embedding")
+        d = cls.dimension_over_reals()
 
         # Use the base ring of the matrix to ensure that its entries can be
         # multiplied by elements of the quaternion algebra.
@@ -1929,10 +2094,10 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
         # quaternion block.
         elements = []
-        for l in range(n/4):
-            for m in range(n/4):
+        for l in range(n/d):
+            for m in range(n/d):
                 submat = ComplexMatrixEuclideanJordanAlgebra.real_unembed(
-                    M[4*l:4*l+4,4*m:4*m+4] )
+                    M[d*l:d*l+d,d*m:d*m+d] )
                 if submat[0,0] != submat[1,1].conjugate():
                     raise ValueError('bad on-diagonal submatrix')
                 if submat[0,1] != -submat[1,0].conjugate():
@@ -1943,42 +2108,11 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
                 z += submat[0,1].imag()*k
                 elements.append(z)
 
-        return matrix(Q, n/4, elements)
-
-
-    @classmethod
-    def matrix_inner_product(cls,X,Y):
-        """
-        Compute a matrix inner product in this algebra directly from
-        its real embedding.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
-
-        TESTS:
-
-        This gives the same answer as the slow, default method implemented
-        in :class:`MatrixEuclideanJordanAlgebra`::
-
-            sage: set_random_seed()
-            sage: J = QuaternionHermitianEJA.random_instance()
-            sage: x,y = J.random_elements(2)
-            sage: Xe = x.to_matrix()
-            sage: Ye = y.to_matrix()
-            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 == expected
-            True
-
-        """
-        return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/4
+        return matrix(Q, n/d, elements)
 
 
-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
@@ -1993,9 +2127,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
 
@@ -2036,7 +2170,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.
 
@@ -2054,11 +2188,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()
 
@@ -2088,16 +2223,25 @@ 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)
+    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(QuaternionHermitianEJA, self).__init__(self._denormalized_basis(n),
+                                                     self.jordan_product,
+                                                     self.trace_inner_product,
+                                                     **kwargs)
+        # TODO: this could be factored out somehow, but is left here
+        # because the MatrixEuclideanJordanAlgebra is not presently
+        # a subclass of the FDEJA class that defines rank() and one().
         self.rank.set_cache(n)
+        idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
+        self.one.set_cache(self(idV))
+
 
     @staticmethod
     def _max_random_instance_size():
@@ -2107,16 +2251,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.
@@ -2156,17 +2299,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)
@@ -2185,16 +2336,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(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 =
@@ -2250,7 +2400,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)
@@ -2259,10 +2411,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)
@@ -2273,40 +2425,32 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
         sage: actual == expected
         True
     """
-    def __init__(self, B, field=AA, **kwargs):
-        n = B.nrows()
-
+    def __init__(self, B, **kwargs):
         if not B.is_positive_definite():
             raise ValueError("bilinear form is not positive-definite")
 
-        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
-
-        super(BilinearFormEJA, self).__init__(field,
-                                              mult_table,
-                                              check_axioms=False,
+        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 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
+
+        n = B.nrows()
+        standard_basis = FreeModule(ZZ, n).basis()
+        super(BilinearFormEJA, self).__init__(standard_basis,
+                                              jordan_product,
+                                              inner_product,
                                               **kwargs)
 
         # The rank of this algebra is two, unless we're in a
@@ -2327,28 +2471,28 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
         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):
@@ -2401,11 +2545,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():
@@ -2415,18 +2566,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.
 
@@ -2455,12 +2605,18 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra,
         0
 
     """
-    def __init__(self, field=AA, **kwargs):
-        mult_table = []
-        self._inner_product_matrix = matrix(field,0)
-        super(TrivialEJA, self).__init__(field,
-                                         mult_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.
@@ -2468,10 +2624,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"""
@@ -2504,8 +2660,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):
         ...
@@ -2522,20 +2678,25 @@ 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 = [ [ 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())
@@ -2553,8 +2714,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,
@@ -2683,8 +2844,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