]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: fix the Hadamard/JordanSpinEJA fast path for charpoly_coefficients.
[sage.d.git] / mjo / eja / eja_algebra.py
index ec3dd1141a9f497acca26f6d09c830420dfb5d5d..3d49005bda1ad3c987c654b11be1ea901da5af7b 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``::
 
@@ -291,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:
@@ -313,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))
@@ -393,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()):
@@ -490,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,orthonormalize=False)
+            sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
             sage: J.coordinate_polynomial_ring()
             Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
 
@@ -1155,16 +1162,26 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr
 
     """
     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
 
@@ -1183,39 +1200,33 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr
 
         V = VectorSpace(field, degree)
 
-        # If we were asked to orthonormalize, and if the orthonormal
-        # basis is different from the given one, then we also want to
-        # compute multiplication and inner-product tables for the
-        # deorthonormalized basis. These can be used later to
-        # construct a deorthonormalized copy of this algebra over QQ
-        # in which several operations are much faster.
+        # Save a copy of an algebra with the original, rational basis
+        # and over QQ where computations are fast.
         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(
-                                          QQ,
-                                          basis,
-                                          jordan_product,
-                                          inner_product,
-                                          orthonormalize=False,
-                                          prefix=prefix,
-                                          category=category,
-                                          check_field=False,
-                                          check_axioms=False)
+        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.
-            W = V.span_of_basis( vector_basis )
+            # 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
@@ -1251,22 +1262,20 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr
             else:
                 vector_basis = gram_schmidt(vector_basis, inner_product)
 
-            W = V.span_of_basis( vector_basis )
-
             # Normalize the "matrix" basis, too!
             basis = vector_basis
 
             if basis_is_matrices:
                 basis = tuple( map(_vec2mat,basis) )
 
-        W = V.span_of_basis( vector_basis )
+        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 )
+        U = V.span_of_basis( deortho_vector_basis, check=check_axioms)
         self._deortho_matrix = matrix( U.coordinate_vector(q)
                                        for q in vector_basis )
 
@@ -1349,9 +1358,11 @@ 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()
 
@@ -1441,7 +1452,22 @@ class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra):
 
 class MatrixEuclideanJordanAlgebra:
     @staticmethod
-    def real_embed(M):
+    def dimension_over_reals():
+        r"""
+        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.
+
+        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.
 
@@ -1454,15 +1480,21 @@ class MatrixEuclideanJordanAlgebra:
           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):
@@ -1470,6 +1502,61 @@ class MatrixEuclideanJordanAlgebra:
 
     @classmethod
     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()
@@ -1486,20 +1573,8 @@ class MatrixEuclideanJordanAlgebra:
 
 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(ConcreteEuclideanJordanAlgebra,
@@ -1526,9 +1601,9 @@ class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra,
 
     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
 
@@ -1569,7 +1644,7 @@ class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra,
 
     """
     @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.
 
@@ -1581,7 +1656,7 @@ class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra,
 
             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
 
@@ -1591,7 +1666,7 @@ class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra,
         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:
@@ -1605,27 +1680,39 @@ class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra,
         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,
+    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)
-        self.one.set_cache(self(matrix.identity(field,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 +
@@ -1667,9 +1754,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()
@@ -1683,8 +1769,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().
 
@@ -1715,31 +1801,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]:
@@ -1747,38 +1839,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
                 z = submat[0,0] + submat[0,1]*i
                 elements.append(z)
 
-        return matrix(F, n/2, elements)
-
-
-    @classmethod
-    def trace_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:
-
-        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.trace_inner_product(Xe,Ye)
-            sage: actual == expected
-            True
-
-        """
-        return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/2
+        return matrix(F, n/d, elements)
 
 
 class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra,
@@ -1797,9 +1858,9 @@ class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra,
 
     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
 
@@ -1841,7 +1902,7 @@ class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra,
     """
 
     @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.
 
@@ -1860,15 +1921,16 @@ class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra,
             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:
         #
@@ -1894,31 +1956,41 @@ class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra,
         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,
-                                               self.jordan_product,
-                                               self.trace_inner_product,
-                                               **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)
+        # 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)
-        # TODO: pre-cache the identity!
+        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
@@ -1957,10 +2029,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()
@@ -1983,8 +2054,8 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
 
 
 
-    @staticmethod
-    def real_unembed(M):
+    @classmethod
+    def real_unembed(cls,M):
         """
         The inverse of _embed_quaternion_matrix().
 
@@ -2014,11 +2085,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.
@@ -2030,10 +2099,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():
@@ -2044,38 +2113,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
                 z += submat[0,1].imag()*k
                 elements.append(z)
 
-        return matrix(Q, n/4, elements)
-
-
-    @classmethod
-    def trace_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.trace_inner_product(Xe,Ye)
-            sage: actual == expected
-            True
-
-        """
-        return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/4
+        return matrix(Q, n/d, elements)
 
 
 class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra,
@@ -2094,9 +2132,9 @@ class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra,
 
     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
 
@@ -2137,7 +2175,7 @@ class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra,
 
     """
     @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.
 
@@ -2155,11 +2193,12 @@ class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra,
 
             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()
 
@@ -2192,15 +2231,22 @@ class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra,
         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,
+    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)
-        # TODO: cache one()!
+        idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
+        self.one.set_cache(self(idV))
+
 
     @staticmethod
     def _max_random_instance_size():
@@ -2210,12 +2256,12 @@ class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra,
         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(ConcreteEuclideanJordanAlgebra):
@@ -2258,17 +2304,25 @@ class HadamardEJA(ConcreteEuclideanJordanAlgebra):
         (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)
@@ -2287,12 +2341,12 @@ class HadamardEJA(ConcreteEuclideanJordanAlgebra):
         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(ConcreteEuclideanJordanAlgebra):
@@ -2376,27 +2430,30 @@ class BilinearFormEJA(ConcreteEuclideanJordanAlgebra):
         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)
@@ -2419,28 +2476,28 @@ class BilinearFormEJA(ConcreteEuclideanJordanAlgebra):
         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):
@@ -2493,11 +2550,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():
@@ -2507,14 +2571,14 @@ 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(ConcreteEuclideanJordanAlgebra):
@@ -2546,12 +2610,16 @@ class TrivialEJA(ConcreteEuclideanJordanAlgebra):
         0
 
     """
-    def __init__(self, field=AA, **kwargs):
+    def __init__(self, **kwargs):
         jordan_product = lambda x,y: x
-        inner_product = lambda x,y: field(0)
+        inner_product = lambda x,y: 0
         basis = ()
-        super(TrivialEJA, self).__init__(field,
-                                         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)
@@ -2561,10 +2629,10 @@ class TrivialEJA(ConcreteEuclideanJordanAlgebra):
         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"""
@@ -2597,8 +2665,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,orthonormalize=False)
+        sage: J1 = random_eja(field=AA)
+        sage: J2 = random_eja(field=QQ,orthonormalize=False)
         sage: J = DirectSumEJA(J1,J2)
         Traceback (most recent call last):
         ...
@@ -2651,8 +2719,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,
@@ -2781,8 +2849,8 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
 
         EXAMPLE::
 
-            sage: J1 = HadamardEJA(3,QQ)
-            sage: J2 = QuaternionHermitianEJA(2,QQ,orthonormalize=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