]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: fix tests and pre-cache ranks.
[sage.d.git] / mjo / eja / eja_algebra.py
index 0f2b655f21795cb1feb0e8b41cfc878261a67962..19db8b0ccefef297623d81721d4a774f99a9d3d1 100644 (file)
@@ -17,7 +17,7 @@ from sage.misc.lazy_import import lazy_import
 from sage.misc.prandom import choice
 from sage.misc.table import table
 from sage.modules.free_module import FreeModule, VectorSpace
-from sage.rings.all import (ZZ, QQ, RR, RLF, CLF,
+from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
                             PolynomialRing,
                             QuadraticField)
 from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement
@@ -26,18 +26,34 @@ lazy_import('mjo.eja.eja_subalgebra',
 from mjo.eja.eja_utils import _mat2vec
 
 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
-    # This is an ugly hack needed to prevent the category framework
-    # from implementing a coercion from our base ring (e.g. the
-    # rationals) into the algebra. First of all -- such a coercion is
-    # nonsense to begin with. But more importantly, it tries to do so
-    # in the category of rings, and since our algebras aren't
-    # associative they generally won't be rings.
-    _no_generic_basering_coercion = True
+
+    def _coerce_map_from_base_ring(self):
+        """
+        Disable the map from the base ring into the algebra.
+
+        Performing a nonsense conversion like this automatically
+        is counterpedagogical. The fallback is to try the usual
+        element constructor, which should also fail.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import random_eja
+
+        TESTS::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: J(1)
+            Traceback (most recent call last):
+            ...
+            ValueError: not a naturally-represented algebra element
+
+        """
+        return None
 
     def __init__(self,
                  field,
                  mult_table,
-                 rank,
                  prefix='e',
                  category=None,
                  natural_basis=None,
@@ -74,7 +90,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
                 # a real embedding.
                 raise ValueError('field is not real')
 
-        self._rank = rank
         self._natural_basis = natural_basis
 
         if category is None:
@@ -94,8 +109,10 @@ 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.
-        self._multiplication_table = [ map(lambda x: self.from_vector(x), ls)
-                                       for ls in mult_table ]
+        self._multiplication_table = [
+            list(map(lambda x: self.from_vector(x), ls))
+            for ls in mult_table
+        ]
 
 
     def _element_constructor_(self, elt):
@@ -109,7 +126,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
-            ....:                                  RealCartesianProductEJA,
+            ....:                                  HadamardEJA,
             ....:                                  RealSymmetricEJA)
 
         EXAMPLES:
@@ -137,7 +154,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         vector representations) back and forth faithfully::
 
             sage: set_random_seed()
-            sage: J = RealCartesianProductEJA.random_instance()
+            sage: J = HadamardEJA.random_instance()
             sage: x = J.random_element()
             sage: J(x.to_vector().column()) == x
             True
@@ -147,15 +164,22 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             True
 
         """
+        msg = "not a naturally-represented algebra element"
         if elt == 0:
             # The superclass implementation of random_element()
             # needs to be able to coerce "0" into the algebra.
             return self.zero()
+        elif elt in self.base_ring():
+            # Ensure that no base ring -> algebra coercion is performed
+            # by this method. There's some stupidity in sage that would
+            # otherwise propagate to this method; for example, sage thinks
+            # that the integer 3 belongs to the space of 2-by-2 matrices.
+            raise ValueError(msg)
 
         natural_basis = self.natural_basis()
         basis_space = natural_basis[0].matrix_space()
         if elt not in basis_space:
-            raise ValueError("not a naturally-represented algebra element")
+            raise ValueError(msg)
 
         # Thanks for nothing! Matrix spaces aren't vector spaces in
         # Sage, so we have to figure out its natural-basis coordinates
@@ -181,8 +205,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         Ensure that it says what we think it says::
 
-            sage: JordanSpinEJA(2, field=QQ)
-            Euclidean Jordan algebra of dimension 2 over Rational Field
+            sage: JordanSpinEJA(2, field=AA)
+            Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
             sage: JordanSpinEJA(3, field=RDF)
             Euclidean Jordan algebra of dimension 3 over Real Double Field
 
@@ -407,7 +431,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         S = PolynomialRing(S, R.variable_names())
         t = S(t)
 
-        return sum( a[k]*(t**k) for k in xrange(len(a)) )
+        return sum( a[k]*(t**k) for k in range(len(a)) )
 
 
     def inner_product(self, x, y):
@@ -493,7 +517,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         """
         M = list(self._multiplication_table) # copy
-        for i in xrange(len(M)):
+        for i in range(len(M)):
             # M had better be "square"
             M[i] = [self.monomial(i)] + M[i]
         M = [["*"] + list(self.gens())] + M
@@ -525,8 +549,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             Finite family {0: e0, 1: e1, 2: e2}
             sage: J.natural_basis()
             (
-            [1 0]  [        0 1/2*sqrt2]  [0 0]
-            [0 0], [1/2*sqrt2         0], [0 1]
+            [1 0]  [                  0 0.7071067811865475?]  [0 0]
+            [0 0], [0.7071067811865475?                   0], [0 1]
             )
 
         ::
@@ -580,12 +604,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
             ....:                                  random_eja)
 
         EXAMPLES::
 
-            sage: J = RealCartesianProductEJA(5)
+            sage: J = HadamardEJA(5)
             sage: J.one()
             e0 + e1 + e2 + e3 + e4
 
@@ -731,7 +755,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
         J1 = trivial                          # eigenvalue one
 
-        for (eigval, eigspace) in c.operator().matrix().left_eigenspaces():
+        for (eigval, eigspace) in c.operator().matrix().right_eigenspaces():
             if eigval == ~(self.base_ring()(2)):
                 J5 = eigspace
             else:
@@ -765,23 +789,26 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             True
 
         """
-        return  tuple( self.random_element() for idx in xrange(count) )
-
+        return tuple( self.random_element() for idx in range(count) )
 
+    @cached_method
     def rank(self):
         """
         Return the rank of this EJA.
 
         ALGORITHM:
 
-        The author knows of no algorithm to compute the rank of an EJA
-        where only the multiplication table is known. In lieu of one, we
-        require the rank to be specified when the algebra is created,
-        and simply pass along that number here.
+        We first compute the polynomial "column matrices" `p_{k}` that
+        evaluate to `x^k` on the coordinates of `x`. Then, we begin
+        adding them to a matrix one at a time, and trying to solve the
+        system that makes `p_{0}`,`p_{1}`,..., `p_{s-1}` add up to
+        `p_{s}`. This will succeed only when `s` is the rank of the
+        algebra, as proven in a recent draft paper of mine.
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            ....:                                  JordanSpinEJA,
             ....:                                  RealSymmetricEJA,
             ....:                                  ComplexHermitianEJA,
             ....:                                  QuaternionHermitianEJA,
@@ -822,8 +849,80 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: r > 0 or (r == 0 and J.is_trivial())
             True
 
+        Ensure that computing the rank actually works, since the ranks
+        of all simple algebras are known and will be cached by default::
+
+            sage: J = HadamardEJA(4)
+            sage: J.rank.clear_cache()
+            sage: J.rank()
+            4
+
+        ::
+
+            sage: J = JordanSpinEJA(4)
+            sage: J.rank.clear_cache()
+            sage: J.rank()
+            2
+
+        ::
+
+            sage: J = RealSymmetricEJA(3)
+            sage: J.rank.clear_cache()
+            sage: J.rank()
+            3
+
+        ::
+
+            sage: J = ComplexHermitianEJA(2)
+            sage: J.rank.clear_cache()
+            sage: J.rank()
+            2
+
+        ::
+
+            sage: J = QuaternionHermitianEJA(2)
+            sage: J.rank.clear_cache()
+            sage: J.rank()
+            2
+
         """
-        return self._rank
+        n = self.dimension()
+        if n == 0:
+            return 0
+        elif n == 1:
+            return 1
+
+        var_names = [ "X" + str(z) for z in range(1,n+1) ]
+        R = PolynomialRing(self.base_ring(), var_names)
+        vars = R.gens()
+
+        def L_x_i_j(i,j):
+            # From a result in my book, these are the entries of the
+            # basis representation of L_x.
+            return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
+                        for k in range(n) )
+
+        L_x = matrix(R, n, n, L_x_i_j)
+        x_powers = [ vars[k]*(L_x**k)*self.one().to_vector()
+                     for k in range(n) ]
+
+        # Can assume n >= 2
+        M = matrix([x_powers[0]])
+        old_rank = 1
+
+        for d in range(1,n):
+            M = matrix(M.rows() + [x_powers[d]])
+            M.echelonize()
+            # TODO: we've basically solved the system here.
+            # We should save the echelonized matrix somehow
+            # so that it can be reused in the charpoly method.
+            new_rank = M.rank()
+            if new_rank == old_rank:
+                return new_rank
+            else:
+                old_rank = new_rank
+
+        return n
 
 
     def vector_space(self):
@@ -885,7 +984,7 @@ class KnownRankEJA(object):
         return 5
 
     @classmethod
-    def random_instance(cls, field=QQ, **kwargs):
+    def random_instance(cls, field=AA, **kwargs):
         """
         Return a random instance of this type of algebra.
 
@@ -901,8 +1000,7 @@ class KnownRankEJA(object):
         return cls(n, field, **kwargs)
 
 
-class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra,
-                              KnownRankEJA):
+class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
     """
     Return the Euclidean Jordan Algebra corresponding to the set
     `R^n` under the Hadamard product.
@@ -913,13 +1011,13 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra,
 
     SETUP::
 
-        sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
+        sage: from mjo.eja.eja_algebra import HadamardEJA
 
     EXAMPLES:
 
     This multiplication table can be verified by hand::
 
-        sage: J = RealCartesianProductEJA(3)
+        sage: J = HadamardEJA(3)
         sage: e0,e1,e2 = J.gens()
         sage: e0*e0
         e0
@@ -938,17 +1036,18 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra,
 
     We can change the generator prefix::
 
-        sage: RealCartesianProductEJA(3, prefix='r').gens()
+        sage: HadamardEJA(3, prefix='r').gens()
         (r0, r1, r2)
 
     """
-    def __init__(self, n, field=QQ, **kwargs):
+    def __init__(self, n, field=AA, **kwargs):
         V = VectorSpace(field, n)
-        mult_table = [ [ V.gen(i)*(i == j) for j in xrange(n) ]
-                       for i in xrange(n) ]
+        mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
+                       for i in range(n) ]
 
-        fdeja = super(RealCartesianProductEJA, self)
-        return fdeja.__init__(field, mult_table, rank=n, **kwargs)
+        fdeja = super(HadamardEJA, self)
+        fdeja.__init__(field, mult_table, **kwargs)
+        self.rank.set_cache(n)
 
     def inner_product(self, x, y):
         """
@@ -956,7 +1055,7 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra,
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
+            sage: from mjo.eja.eja_algebra import HadamardEJA
 
         TESTS:
 
@@ -964,7 +1063,7 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra,
         over `R^n`::
 
             sage: set_random_seed()
-            sage: J = RealCartesianProductEJA.random_instance()
+            sage: J = HadamardEJA.random_instance()
             sage: x,y = J.random_elements(2)
             sage: X = x.natural_representation()
             sage: Y = y.natural_representation()
@@ -975,7 +1074,7 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra,
         return x.to_vector().inner_product(y.to_vector())
 
 
-def random_eja(field=QQ, nontrivial=False):
+def random_eja(field=AA, nontrivial=False):
     """
     Return a "random" finite-dimensional Euclidean Jordan Algebra.
 
@@ -1007,7 +1106,7 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         # field can have dimension 4 (quaternions) too.
         return 2
 
-    def __init__(self, field, basis, rank, normalize_basis=True, **kwargs):
+    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
@@ -1020,7 +1119,7 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         # time to ensure that it isn't a generator expression.
         basis = tuple(basis)
 
-        if rank > 1 and normalize_basis:
+        if len(basis) > 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.
@@ -1037,12 +1136,31 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         Qs = self.multiplication_table_from_matrix_basis(basis)
 
         fdeja = super(MatrixEuclideanJordanAlgebra, self)
-        return fdeja.__init__(field,
-                              Qs,
-                              rank=rank,
-                              natural_basis=basis,
-                              **kwargs)
+        fdeja.__init__(field, Qs, natural_basis=basis, **kwargs)
+        return
+
+
+    @cached_method
+    def rank(self):
+        r"""
+        Override the parent method with something that tries to compute
+        over a faster (non-extension) field.
+        """
+        if self._basis_normalizers is None:
+            # We didn't normalize, so assume that the basis we started
+            # with had entries in a nice field.
+            return super(MatrixEuclideanJordanAlgebra, self).rank()
+        else:
+            basis = ( (b/n) for (b,n) in zip(self.natural_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.
+            J = MatrixEuclideanJordanAlgebra(QQ,
+                                             basis,
+                                             normalize_basis=False)
+            return J.rank()
 
     @cached_method
     def _charpoly_coeff(self, i):
@@ -1061,7 +1179,6 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
             # Do this over the rationals and convert back at the end.
             J = MatrixEuclideanJordanAlgebra(QQ,
                                              basis,
-                                             self.rank(),
                                              normalize_basis=False)
             (_,x,_,_) = J._charpoly_matrix_system()
             p = J._charpoly_coeff(i)
@@ -1099,9 +1216,9 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         V = VectorSpace(field, dimension**2)
         W = V.span_of_basis( _mat2vec(s) for s in basis )
         n = len(basis)
-        mult_table = [[W.zero() for j in xrange(n)] for i in xrange(n)]
-        for i in xrange(n):
-            for j in xrange(n):
+        mult_table = [[W.zero() for j in range(n)] for i in range(n)]
+        for i in range(n):
+            for j in range(n):
                 mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
                 mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
 
@@ -1195,8 +1312,8 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
 
     In theory, our "field" can be any subfield of the reals::
 
-        sage: RealSymmetricEJA(2, AA)
-        Euclidean Jordan algebra of dimension 3 over Algebraic Real Field
+        sage: RealSymmetricEJA(2, RDF)
+        Euclidean Jordan algebra of dimension 3 over Real Double Field
         sage: RealSymmetricEJA(2, RR)
         Euclidean Jordan algebra of dimension 3 over Real Field with
         53 bits of precision
@@ -1272,8 +1389,8 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
         # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
         # coordinates.
         S = []
-        for i in xrange(n):
-            for j in xrange(i+1):
+        for i in range(n):
+            for j in range(i+1):
                 Eij = matrix(field, n, lambda k,l: k==i and l==j)
                 if i == j:
                     Sij = Eij
@@ -1288,9 +1405,10 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
         return 4 # Dimension 10
 
 
-    def __init__(self, n, field=QQ, **kwargs):
+    def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n, field)
-        super(RealSymmetricEJA, self).__init__(field, basis, n, **kwargs)
+        super(RealSymmetricEJA, self).__init__(field, basis, **kwargs)
+        self.rank.set_cache(n)
 
 
 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
@@ -1308,7 +1426,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
 
         EXAMPLES::
 
-            sage: F = QuadraticField(-1, 'i')
+            sage: F = QuadraticField(-1, 'I')
             sage: x1 = F(4 - 2*i)
             sage: x2 = F(1 + 2*i)
             sage: x3 = F(-i)
@@ -1328,7 +1446,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             sage: set_random_seed()
             sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
             sage: n = ZZ.random_element(n_max)
-            sage: F = QuadraticField(-1, 'i')
+            sage: F = QuadraticField(-1, 'I')
             sage: X = random_matrix(F, n)
             sage: Y = random_matrix(F, n)
             sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
@@ -1371,15 +1489,15 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             ....:                 [ 9,  10, 11, 12],
             ....:                 [-10, 9, -12, 11] ])
             sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
-            [  2*i + 1   4*i + 3]
-            [ 10*i + 9 12*i + 11]
+            [  2*I + 1   4*I + 3]
+            [ 10*I + 9 12*I + 11]
 
         TESTS:
 
         Unembedding is the inverse of embedding::
 
             sage: set_random_seed()
-            sage: F = QuadraticField(-1, 'i')
+            sage: F = QuadraticField(-1, 'I')
             sage: M = random_matrix(F, 3)
             sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
             sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
@@ -1397,14 +1515,19 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         field = M.base_ring()
         R = PolynomialRing(field, 'z')
         z = R.gen()
-        F = field.extension(z**2 + 1, 'i', embedding=CLF(-1).sqrt())
+        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
+        else:
+            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 xrange(n/2):
-            for j in xrange(n/2):
+        for k in range(n/2):
+            for j in range(n/2):
                 submat = M[2*k:2*k+2,2*j:2*j+2]
                 if submat[0,0] != submat[1,1]:
                     raise ValueError('bad on-diagonal submatrix')
@@ -1438,7 +1561,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             sage: Ye = y.natural_representation()
             sage: X = ComplexHermitianEJA.real_unembed(Xe)
             sage: Y = ComplexHermitianEJA.real_unembed(Ye)
-            sage: expected = (X*Y).trace().vector()[0]
+            sage: expected = (X*Y).trace().real()
             sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
             sage: actual == expected
             True
@@ -1462,8 +1585,8 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
 
     In theory, our "field" can be any subfield of the reals::
 
-        sage: ComplexHermitianEJA(2, AA)
-        Euclidean Jordan algebra of dimension 4 over Algebraic Real Field
+        sage: ComplexHermitianEJA(2, RDF)
+        Euclidean Jordan algebra of dimension 4 over Real Double Field
         sage: ComplexHermitianEJA(2, RR)
         Euclidean Jordan algebra of dimension 4 over Real Field with
         53 bits of precision
@@ -1555,8 +1678,8 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
         #   * The diagonal will (as a result) be real.
         #
         S = []
-        for i in xrange(n):
-            for j in xrange(i+1):
+        for i in range(n):
+            for j in range(i+1):
                 Eij = matrix(F, n, lambda k,l: k==i and l==j)
                 if i == j:
                     Sij = cls.real_embed(Eij)
@@ -1573,9 +1696,10 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
         return ( s.change_ring(field) for s in S )
 
 
-    def __init__(self, n, field=QQ, **kwargs):
+    def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n,field)
-        super(ComplexHermitianEJA,self).__init__(field, basis, n, **kwargs)
+        super(ComplexHermitianEJA,self).__init__(field, basis, **kwargs)
+        self.rank.set_cache(n)
 
 
 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
@@ -1625,7 +1749,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         if M.ncols() != n:
             raise ValueError("the matrix 'M' must be square")
 
-        F = QuadraticField(-1, 'i')
+        F = QuadraticField(-1, 'I')
         i = F.gen()
 
         blocks = []
@@ -1693,18 +1817,18 @@ 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 xrange(n/4):
-            for m in xrange(n/4):
+        for l in range(n/4):
+            for m in range(n/4):
                 submat = ComplexMatrixEuclideanJordanAlgebra.real_unembed(
                     M[4*l:4*l+4,4*m:4*m+4] )
                 if submat[0,0] != submat[1,1].conjugate():
                     raise ValueError('bad on-diagonal submatrix')
                 if submat[0,1] != -submat[1,0].conjugate():
                     raise ValueError('bad off-diagonal submatrix')
-                z  = submat[0,0].vector()[0]   # real part
-                z += submat[0,0].vector()[1]*i # imag part
-                z += submat[0,1].vector()[0]*j # real part
-                z += submat[0,1].vector()[1]*k # imag part
+                z  = submat[0,0].real()
+                z += submat[0,0].imag()*i
+                z += submat[0,1].real()*j
+                z += submat[0,1].imag()*k
                 elements.append(z)
 
         return matrix(Q, n/4, elements)
@@ -1757,8 +1881,8 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
 
     In theory, our "field" can be any subfield of the reals::
 
-        sage: QuaternionHermitianEJA(2, AA)
-        Euclidean Jordan algebra of dimension 6 over Algebraic Real Field
+        sage: QuaternionHermitianEJA(2, RDF)
+        Euclidean Jordan algebra of dimension 6 over Real Double Field
         sage: QuaternionHermitianEJA(2, RR)
         Euclidean Jordan algebra of dimension 6 over Real Field with
         53 bits of precision
@@ -1846,8 +1970,8 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
         #   * The diagonal will (as a result) be real.
         #
         S = []
-        for i in xrange(n):
-            for j in xrange(i+1):
+        for i in range(n):
+            for j in range(i+1):
                 Eij = matrix(Q, n, lambda k,l: k==i and l==j)
                 if i == j:
                     Sij = cls.real_embed(Eij)
@@ -1869,16 +1993,136 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
         return ( s.change_ring(field) for s in S )
 
 
-    def __init__(self, n, field=QQ, **kwargs):
+    def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n,field)
-        super(QuaternionHermitianEJA,self).__init__(field, basis, n, **kwargs)
+        super(QuaternionHermitianEJA,self).__init__(field, basis, **kwargs)
+        self.rank.set_cache(n)
+
+
+class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
+    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 =
+    (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
+    symmetric positive-definite "bilinear form" matrix. It has
+    dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
+    when ``B`` is the identity matrix of order ``n-1``.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
+        ....:                                  JordanSpinEJA)
+
+    EXAMPLES:
 
+    When no bilinear form is specified, the identity matrix is used,
+    and the resulting algebra is the Jordan spin algebra::
 
-class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
+        sage: J0 = BilinearFormEJA(3)
+        sage: J1 = JordanSpinEJA(3)
+        sage: J0.multiplication_table() == J0.multiplication_table()
+        True
+
+    TESTS:
+
+    We can create a zero-dimensional algebra::
+
+        sage: J = BilinearFormEJA(0)
+        sage: J.basis()
+        Finite family {}
+
+    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::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(5)
+        sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
+        sage: B = M.transpose()*M
+        sage: J = BilinearFormEJA(n, B=B)
+        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()))
+        ....:         for ei in eis ]
+        sage: actual = [ sis[i]*sis[j]
+        ....:            for i in range(n-1)
+        ....:            for j in range(n-1) ]
+        sage: expected = [ J.one() if i == j else J.zero()
+        ....:              for i in range(n-1)
+        ....:              for j in range(n-1) ]
+        sage: actual == expected
+        True
+    """
+    def __init__(self, n, field=AA, B=None, **kwargs):
+        if B is None:
+            self._B = matrix.identity(field, max(0,n-1))
+        else:
+            self._B = B
+
+        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 = x0*y0 + (self._B*xbar).inner_product(ybar)
+                zbar = y0*xbar + x0*ybar
+                z = V([z0] + zbar.list())
+                mult_table[i][j] = z
+
+        # The rank of this algebra is two, unless we're in a
+        # one-dimensional ambient space (because the rank is bounded
+        # by the ambient dimension).
+        fdeja = super(BilinearFormEJA, self)
+        fdeja.__init__(field, mult_table, **kwargs)
+        self.rank.set_cache(min(n,2))
+
+    def inner_product(self, x, y):
+        r"""
+        Half of the trace inner product.
+
+        This is defined so that the special case of the Jordan spin
+        algebra gets the usual inner product.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import BilinearFormEJA
+
+        TESTS:
+
+        Ensure that this is one-half of the trace inner-product when
+        the algebra isn't just the reals (when ``n`` isn't one). This
+        is in Faraut and Koranyi, and also my "On the symmetry..."
+        paper::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(2,5)
+            sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
+            sage: B = M.transpose()*M
+            sage: J = BilinearFormEJA(n, B=B)
+            sage: x = J.random_element()
+            sage: y = J.random_element()
+            sage: x.inner_product(y) == (x*y).trace()/2
+            True
+
+        """
+        xvec = x.to_vector()
+        xbar = xvec[1:]
+        yvec = y.to_vector()
+        ybar = yvec[1:]
+        return x[0]*y[0] + (self._B*xbar).inner_product(ybar)
+
+
+class JordanSpinEJA(BilinearFormEJA):
     """
     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
     with the usual inner product and jordan product ``x*y =
-    (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
+    (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
     the reals.
 
     SETUP::
@@ -1911,42 +2155,9 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
         sage: JordanSpinEJA(2, prefix='B').gens()
         (B0, B1)
 
-    """
-    def __init__(self, n, field=QQ, **kwargs):
-        V = VectorSpace(field, n)
-        mult_table = [[V.zero() for j in xrange(n)] for i in xrange(n)]
-        for i in xrange(n):
-            for j in xrange(n):
-                x = V.gen(i)
-                y = V.gen(j)
-                x0 = x[0]
-                xbar = x[1:]
-                y0 = y[0]
-                ybar = y[1:]
-                # z = x*y
-                z0 = x.inner_product(y)
-                zbar = y0*xbar + x0*ybar
-                z = V([z0] + zbar.list())
-                mult_table[i][j] = z
-
-        # The rank of the spin algebra is two, unless we're in a
-        # one-dimensional ambient space (because the rank is bounded by
-        # the ambient dimension).
-        fdeja = super(JordanSpinEJA, self)
-        return fdeja.__init__(field, mult_table, rank=min(n,2), **kwargs)
-
-    def inner_product(self, x, y):
-        """
-        Faster to reimplement than to use natural representations.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import JordanSpinEJA
-
-        TESTS:
+    TESTS:
 
-        Ensure that this is the usual inner product for the algebras
-        over `R^n`::
+        Ensure that we have the usual inner product on `R^n`::
 
             sage: set_random_seed()
             sage: J = JordanSpinEJA.random_instance()
@@ -1956,8 +2167,11 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
             sage: x.inner_product(y) == J.natural_inner_product(X,Y)
             True
 
-        """
-        return x.to_vector().inner_product(y.to_vector())
+    """
+    def __init__(self, n, field=AA, **kwargs):
+        # This is a special case of the BilinearFormEJA with the identity
+        # matrix as its bilinear form.
+        return super(JordanSpinEJA, self).__init__(n, field, **kwargs)
 
 
 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
@@ -1984,14 +2198,15 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
         sage: J.one().norm()
         0
         sage: J.one().subalgebra_generated_by()
-        Euclidean Jordan algebra of dimension 0 over Rational Field
+        Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
         sage: J.rank()
         0
 
     """
-    def __init__(self, field=QQ, **kwargs):
+    def __init__(self, field=AA, **kwargs):
         mult_table = []
         fdeja = super(TrivialEJA, self)
         # The rank is zero using my definition, namely the dimension of the
         # largest subalgebra generated by any element.
-        return fdeja.__init__(field, mult_table, rank=0, **kwargs)
+        fdeja.__init__(field, mult_table, **kwargs)
+        self.rank.set_cache(0)