eja: begin generalizing the charpoly-over-QQ optimizations.

[sage.d.git] / mjo / eja / eja_algebra.py

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index c131c5f6b612a6ffcae8e1d3357264e5f016176e..1fc3618005eb0ac8be12e0e3e09849ab6f983819 100644 (file)
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -5,7 +5,7 @@ are used in optimization, and have some additional nice methods beyond
 what can be supported in a general Jordan Algebra.
 """
 
 what can be supported in a general Jordan Algebra.
 """
 

-from itertools import izip, repeat
+from itertools import repeat

 
 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
 from sage.categories.magmatic_algebras import MagmaticAlgebras
 
 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
 from sage.categories.magmatic_algebras import MagmaticAlgebras


@@ -13,38 +13,59 @@ from sage.combinat.free_module import CombinatorialFreeModule
 from sage.matrix.constructor import matrix
 from sage.matrix.matrix_space import MatrixSpace
 from sage.misc.cachefunc import cached_method
 from sage.matrix.constructor import matrix
 from sage.matrix.matrix_space import MatrixSpace
 from sage.misc.cachefunc import cached_method

+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.misc.prandom import choice
 from sage.misc.table import table
 from sage.modules.free_module import FreeModule, VectorSpace

-from sage.rings.integer_ring import ZZ
-from sage.rings.number_field.number_field import NumberField, QuadraticField
-from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
-from sage.rings.rational_field import QQ
-from sage.rings.real_lazy import CLF, RLF
-
+from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
+                            PolynomialRing,
+                            QuadraticField)

 from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement
 from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement

+lazy_import('mjo.eja.eja_subalgebra',
+            'FiniteDimensionalEuclideanJordanSubalgebra')

 from mjo.eja.eja_utils import _mat2vec
 
 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 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,
 
     def __init__(self,
                  field,
                  mult_table,

-                 rank,

                  prefix='e',
                  category=None,
                  prefix='e',
                  category=None,

-                 natural_basis=None):
+                 natural_basis=None,
+                 check_field=True,
+                 check_axioms=True):

         """
         SETUP::
 
         """
         SETUP::
 

-            sage: from mjo.eja.eja_algebra import random_eja
+            sage: from mjo.eja.eja_algebra import (
+            ....:   FiniteDimensionalEuclideanJordanAlgebra,
+            ....:   JordanSpinEJA,
+            ....:   random_eja)

 
         EXAMPLES:
 
 
         EXAMPLES:
 


@@ -56,8 +77,36 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: x*y == y*x
             True
 
             sage: x*y == y*x
             True
 

+        TESTS:
+
+        The ``field`` we're given must be real with ``check_field=True``::
+
+            sage: JordanSpinEJA(2,QQbar)
+            Traceback (most recent call last):
+            ...
+            ValueError: scalar field is not real
+
+        The multiplication table must be square with ``check_axioms=True``::
+
+            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()))
+            Traceback (most recent call last):
+            ...
+            ValueError: multiplication table is not square
+

         """
         """

-        self._rank = rank
+        if check_field:
+            if not field.is_subring(RR):
+                # Note: this does return true for the real algebraic
+                # field, the rationals, and any quadratic field where
+                # we've specified a real embedding.
+                raise ValueError("scalar field is not real")
+
+        # The multiplication table had better be square
+        n = len(mult_table)
+        if check_axioms:
+            if not all( len(l) == n for l in mult_table ):
+                raise ValueError("multiplication table is not square")
+

         self._natural_basis = natural_basis
 
         if category is None:
         self._natural_basis = natural_basis
 
         if category is None:


@@ -66,7 +115,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
         fda.__init__(field,
 
         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
         fda.__init__(field,

-                     range(len(mult_table)),
+                     range(n),

                      prefix=prefix,
                      category=category)
         self.print_options(bracket='')
                      prefix=prefix,
                      category=category)
         self.print_options(bracket='')


@@ -77,9 +126,18 @@ 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.
         # 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
+        ]
+
+        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():
+                raise ValueError("inner product is not associative")

 
     def _element_constructor_(self, elt):
         """
 
     def _element_constructor_(self, elt):
         """


@@ -92,7 +150,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,

-            ....:                                  RealCartesianProductEJA,
+            ....:                                  HadamardEJA,

             ....:                                  RealSymmetricEJA)
 
         EXAMPLES:
             ....:                                  RealSymmetricEJA)
 
         EXAMPLES:


@@ -120,7 +178,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         vector representations) back and forth faithfully::
 
             sage: set_random_seed()
         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
             sage: x = J.random_element()
             sage: J(x.to_vector().column()) == x
             True


@@ -130,15 +188,22 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             True
 
         """
             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()
         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:
 
         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
 
         # Thanks for nothing! Matrix spaces aren't vector spaces in
         # Sage, so we have to figure out its natural-basis coordinates


@@ -151,6 +216,24 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         coords =  W.coordinate_vector(_mat2vec(elt))
         return self.from_vector(coords)
 
         coords =  W.coordinate_vector(_mat2vec(elt))
         return self.from_vector(coords)
 

+    @staticmethod
+    def _max_test_case_size():
+        """
+        Return an integer "size" that is an upper bound on the size of
+        this algebra when it is used in a random test
+        case. Unfortunately, the term "size" is quite vague -- when
+        dealing with `R^n` under either the Hadamard or Jordan spin
+        product, the "size" refers to the dimension `n`. When dealing
+        with a matrix algebra (real symmetric or complex/quaternion
+        Hermitian), it refers to the size of the matrix, which is
+        far less than the dimension of the underlying vector space.
+
+        We default to five in this class, which is safe in `R^n`. The
+        matrix algebra subclasses (or any class where the "size" is
+        interpreted to be far less than the dimension) should override
+        with a smaller number.
+        """
+        return 5

 
     def _repr_(self):
         """
 
     def _repr_(self):
         """


@@ -164,8 +247,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         Ensure that it says what we think it says::
 
 
         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
 
             sage: JordanSpinEJA(3, field=RDF)
             Euclidean Jordan algebra of dimension 3 over Real Double Field
 


@@ -176,165 +259,74 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
     def product_on_basis(self, i, j):
         return self._multiplication_table[i][j]
 
     def product_on_basis(self, i, j):
         return self._multiplication_table[i][j]
 

-    def _a_regular_element(self):
-        """
-        Guess a regular element. Needed to compute the basis for our
-        characteristic polynomial coefficients.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import random_eja
-
-        TESTS:
-
-        Ensure that this hacky method succeeds for every algebra that we
-        know how to construct::
-
-            sage: set_random_seed()
-            sage: J = random_eja()
-            sage: J._a_regular_element().is_regular()
-            True
-
-        """
-        gs = self.gens()
-        z = self.sum( (i+1)*gs[i] for i in range(len(gs)) )
-        if not z.is_regular():
-            raise ValueError("don't know a regular element")
-        return z
-
-
-    @cached_method
-    def _charpoly_basis_space(self):
-        """
-        Return the vector space spanned by the basis used in our
-        characteristic polynomial coefficients. This is used not only to
-        compute those coefficients, but also any time we need to
-        evaluate the coefficients (like when we compute the trace or
-        determinant).
-        """
-        z = self._a_regular_element()
-        # Don't use the parent vector space directly here in case this
-        # happens to be a subalgebra. In that case, we would be e.g.
-        # two-dimensional but span_of_basis() would expect three
-        # coordinates.
-        V = VectorSpace(self.base_ring(), self.vector_space().dimension())
-        basis = [ (z**k).to_vector() for k in range(self.rank()) ]
-        V1 = V.span_of_basis( basis )
-        b =  (V1.basis() + V1.complement().basis())
-        return V.span_of_basis(b)
-
-
-
-    @cached_method
-    def _charpoly_coeff(self, i):
-        """
-        Return the coefficient polynomial "a_{i}" of this algebra's
-        general characteristic polynomial.
-
-        Having this be a separate cached method lets us compute and
-        store the trace/determinant (a_{r-1} and a_{0} respectively)
-        separate from the entire characteristic polynomial.
-        """
-        (A_of_x, x, xr, detA) = self._charpoly_matrix_system()
-        R = A_of_x.base_ring()
-        if i >= self.rank():
-            # Guaranteed by theory
-            return R.zero()
-
-        # Danger: the in-place modification is done for performance
-        # reasons (reconstructing a matrix with huge polynomial
-        # entries is slow), but I don't know how cached_method works,
-        # so it's highly possible that we're modifying some global
-        # list variable by reference, here. In other words, you
-        # probably shouldn't call this method twice on the same
-        # algebra, at the same time, in two threads
-        Ai_orig = A_of_x.column(i)
-        A_of_x.set_column(i,xr)
-        numerator = A_of_x.det()
-        A_of_x.set_column(i,Ai_orig)
-
-        # We're relying on the theory here to ensure that each a_i is
-        # indeed back in R, and the added negative signs are to make
-        # the whole charpoly expression sum to zero.
-        return R(-numerator/detA)
-
-
-    @cached_method
-    def _charpoly_matrix_system(self):
-        """
-        Compute the matrix whose entries A_ij are polynomials in
-        X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
-        corresponding to `x^r` and the determinent of the matrix A =
-        [A_ij]. In other words, all of the fixed (cachable) data needed
-        to compute the coefficients of the characteristic polynomial.
-        """
-        r = self.rank()
-        n = self.dimension()
-
-        # Turn my vector space into a module so that "vectors" can
-        # have multivatiate polynomial entries.
-        names = tuple('X' + str(i) for i in range(1,n+1))
-        R = PolynomialRing(self.base_ring(), names)
-
-        # Using change_ring() on the parent's vector space doesn't work
-        # here because, in a subalgebra, that vector space has a basis
-        # and change_ring() tries to bring the basis along with it. And
-        # that doesn't work unless the new ring is a PID, which it usually
-        # won't be.
-        V = FreeModule(R,n)
-
-        # Now let x = (X1,X2,...,Xn) be the vector whose entries are
-        # indeterminates...
-        x = V(names)
-
-        # And figure out the "left multiplication by x" matrix in
-        # that setting.
-        lmbx_cols = []
-        monomial_matrices = [ self.monomial(i).operator().matrix()
-                              for i in range(n) ] # don't recompute these!
-        for k in range(n):
-            ek = self.monomial(k).to_vector()
-            lmbx_cols.append(
-              sum( x[i]*(monomial_matrices[i]*ek)
-                   for i in range(n) ) )
-        Lx = matrix.column(R, lmbx_cols)
-
-        # Now we can compute powers of x "symbolically"
-        x_powers = [self.one().to_vector(), x]
-        for d in range(2, r+1):
-            x_powers.append( Lx*(x_powers[-1]) )
-
-        idmat = matrix.identity(R, n)
-
-        W = self._charpoly_basis_space()
-        W = W.change_ring(R.fraction_field())
-
-        # Starting with the standard coordinates x = (X1,X2,...,Xn)
-        # and then converting the entries to W-coordinates allows us
-        # to pass in the standard coordinates to the charpoly and get
-        # back the right answer. Specifically, with x = (X1,X2,...,Xn),
-        # we have
-        #
-        #   W.coordinates(x^2) eval'd at (standard z-coords)
-        #     =
-        #   W-coords of (z^2)
-        #     =
-        #   W-coords of (standard coords of x^2 eval'd at std-coords of z)
-        #
-        # We want the middle equivalent thing in our matrix, but use
-        # the first equivalent thing instead so that we can pass in
-        # standard coordinates.
-        x_powers = [ W.coordinate_vector(xp) for xp in x_powers ]
-        l2 = [idmat.column(k-1) for k in range(r+1, n+1)]
-        A_of_x = matrix.column(R, n, (x_powers[:r] + l2))
-        return (A_of_x, x, x_powers[r], A_of_x.det())
-
+    def _is_commutative(self):
+        r"""
+        Whether or not this algebra's multiplication table is commutative.
+
+        This method should of course always return ``True``, unless
+        this algebra was constructed with ``check_axioms=False`` and
+        passed an invalid multiplication table.
+        """
+        return all( self.product_on_basis(i,j) == self.product_on_basis(i,j)
+                    for i in range(self.dimension())
+                    for j in range(self.dimension()) )
+
+    def _is_jordanian(self):
+        r"""
+        Whether or not this algebra's multiplication table respects the
+        Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
+
+        We only check one arrangement of `x` and `y`, so for a
+        ``True`` result to be truly true, you should also check
+        :meth:`_is_commutative`. This method should of course always
+        return ``True``, unless this algebra was constructed with
+        ``check_axioms=False`` and passed an invalid multiplication table.
+        """
+        return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
+                    ==
+                    (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j))
+                    for i in range(self.dimension())
+                    for j in range(self.dimension()) )
+
+    def _inner_product_is_associative(self):
+        r"""
+        Return whether or not this algebra's inner product `B` is
+        associative; that is, whether or not `B(xy,z) = B(x,yz)`.
+
+        This method should of course always return ``True``, unless
+        this algebra was constructed with ``check_axioms=False`` and
+        passed an invalid multiplication table.
+        """
+
+        # 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.
+        epsilon = 1e-16
+
+        for i in range(self.dimension()):
+            for j in range(self.dimension()):
+                for k in range(self.dimension()):
+                    x = self.monomial(i)
+                    y = self.monomial(j)
+                    z = self.monomial(k)
+                    diff = (x*y).inner_product(z) - x.inner_product(y*z)
+
+                    if self.base_ring().is_exact():
+                        if diff != 0:
+                            return False
+                    else:
+                        if diff.abs() > epsilon:
+                            return False
+
+        return True

 
     @cached_method
 
     @cached_method

-    def characteristic_polynomial(self):
+    def characteristic_polynomial_of(self):

         """
         """

-        Return a characteristic polynomial that works for all elements
-        of this algebra.
+        Return the algebra's "characteristic polynomial of" function,
+        which is itself a multivariate polynomial that, when evaluated
+        at the coordinates of some algebra element, returns that
+        element's characteristic polynomial.

 
         The resulting polynomial has `n+1` variables, where `n` is the
         dimension of this algebra. The first `n` variables correspond to
 
         The resulting polynomial has `n+1` variables, where `n` is the
         dimension of this algebra. The first `n` variables correspond to


@@ -346,7 +338,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         SETUP::
 
 
         SETUP::
 

-            sage: from mjo.eja.eja_algebra import JordanSpinEJA
+            sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA

 
         EXAMPLES:
 
 
         EXAMPLES:
 


@@ -354,41 +346,41 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         Alizadeh, Example 11.11::
 
             sage: J = JordanSpinEJA(3)
         Alizadeh, Example 11.11::
 
             sage: J = JordanSpinEJA(3)

-            sage: p = J.characteristic_polynomial(); p
+            sage: p = J.characteristic_polynomial_of(); p

             X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
             sage: xvec = J.one().to_vector()
             sage: p(*xvec)
             t^2 - 2*t + 1
 
             X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
             sage: xvec = J.one().to_vector()
             sage: p(*xvec)
             t^2 - 2*t + 1
 

+        By definition, the characteristic polynomial is a monic
+        degree-zero polynomial in a rank-zero algebra. Note that
+        Cayley-Hamilton is indeed satisfied since the polynomial
+        ``1`` evaluates to the identity element of the algebra on
+        any argument::
+
+            sage: J = TrivialEJA()
+            sage: J.characteristic_polynomial_of()
+            1
+

         """
         r = self.rank()
         n = self.dimension()
 
         """
         r = self.rank()
         n = self.dimension()
 

-        # The list of coefficient polynomials a_1, a_2, ..., a_n.
-        a = [ self._charpoly_coeff(i) for i in range(n) ]
+        # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
+        a = self._charpoly_coefficients()

 
         # We go to a bit of trouble here to reorder the
         # indeterminates, so that it's easier to evaluate the
         # characteristic polynomial at x's coordinates and get back
         # something in terms of t, which is what we want.
 
         # We go to a bit of trouble here to reorder the
         # indeterminates, so that it's easier to evaluate the
         # characteristic polynomial at x's coordinates and get back
         # something in terms of t, which is what we want.

-        R = a[0].parent()

         S = PolynomialRing(self.base_ring(),'t')
         t = S.gen(0)
         S = PolynomialRing(self.base_ring(),'t')
         t = S.gen(0)

-        S = PolynomialRing(S, R.variable_names())
-        t = S(t)
-
-        # Note: all entries past the rth should be zero. The
-        # coefficient of the highest power (x^r) is 1, but it doesn't
-        # appear in the solution vector which contains coefficients
-        # for the other powers (to make them sum to x^r).
-        if (r < n):
-            a[r] = 1 # corresponds to x^r
-        else:
-            # When the rank is equal to the dimension, trying to
-            # assign a[r] goes out-of-bounds.
-            a.append(1) # corresponds to x^r
+        if r > 0:
+            R = a[0].parent()
+            S = PolynomialRing(S, R.variable_names())
+            t = S(t)

 
 

-        return sum( a[k]*(t**k) for k in xrange(len(a)) )
+        return (t**r + sum( a[k]*(t**k) for k in range(r) ))

 
 
     def inner_product(self, x, y):
 
 
     def inner_product(self, x, y):


@@ -428,15 +420,19 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         SETUP::
 
 
         SETUP::
 

-            sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
+            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
+            ....:                                  TrivialEJA)

 
         EXAMPLES::
 
             sage: J = ComplexHermitianEJA(3)
             sage: J.is_trivial()
             False
 
         EXAMPLES::
 
             sage: J = ComplexHermitianEJA(3)
             sage: J.is_trivial()
             False

-            sage: A = J.zero().subalgebra_generated_by()
-            sage: A.is_trivial()
+
+        ::
+
+            sage: J = TrivialEJA()
+            sage: J.is_trivial()

             True
 
         """
             True
 
         """


@@ -470,7 +466,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         """
         M = list(self._multiplication_table) # copy
 
         """
         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
             # M had better be "square"
             M[i] = [self.monomial(i)] + M[i]
         M = [["*"] + list(self.gens())] + M


@@ -502,8 +498,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             Finite family {0: e0, 1: e1, 2: e2}
             sage: J.natural_basis()
             (
             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]

             )
 
         ::
             )
 
         ::


@@ -529,8 +525,15 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         """
         Return the matrix space in which this algebra's natural basis
         elements live.
         """
         Return the matrix space in which this algebra's natural basis
         elements live.

+
+        Generally this will be an `n`-by-`1` column-vector space,
+        except when the algebra is trivial. There it's `n`-by-`n`
+        (where `n` is zero), to ensure that two elements of the
+        natural basis space (empty matrices) can be multiplied.

         """
         """

-        if self._natural_basis is None or len(self._natural_basis) == 0:
+        if self.is_trivial():
+            return MatrixSpace(self.base_ring(), 0)
+        elif self._natural_basis is None or len(self._natural_basis) == 0:

             return MatrixSpace(self.base_ring(), self.dimension(), 1)
         else:
             return self._natural_basis[0].matrix_space()
             return MatrixSpace(self.base_ring(), self.dimension(), 1)
         else:
             return self._natural_basis[0].matrix_space()


@@ -557,12 +560,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         SETUP::
 
 
         SETUP::
 

-            sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,

             ....:                                  random_eja)
 
         EXAMPLES::
 
             ....:                                  random_eja)
 
         EXAMPLES::
 

-            sage: J = RealCartesianProductEJA(5)
+            sage: J = HadamardEJA(5)

             sage: J.one()
             e0 + e1 + e2 + e3 + e4
 
             sage: J.one()
             e0 + e1 + e2 + e3 + e4
 


@@ -611,18 +614,206 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return self.linear_combination(zip(self.gens(), coeffs))
 
 
         return self.linear_combination(zip(self.gens(), coeffs))
 
 

-    def random_element(self):
-        # Temporary workaround for https://trac.sagemath.org/ticket/28327
-        if self.is_trivial():
-            return self.zero()
-        else:
-            s = super(FiniteDimensionalEuclideanJordanAlgebra, self)
-            return s.random_element()
+    def peirce_decomposition(self, c):
+        """
+        The Peirce decomposition of this algebra relative to the
+        idempotent ``c``.
+
+        In the future, this can be extended to a complete system of
+        orthogonal idempotents.
+
+        INPUT:
+
+          - ``c`` -- an idempotent of this algebra.
+
+        OUTPUT:
+
+        A triple (J0, J5, J1) containing two subalgebras and one subspace
+        of this algebra,
+
+          - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
+            corresponding to the eigenvalue zero.
+
+          - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
+            corresponding to the eigenvalue one-half.
+
+          - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
+            corresponding to the eigenvalue one.
+
+        These are the only possible eigenspaces for that operator, and this
+        algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
+        orthogonal, and are subalgebras of this algebra with the appropriate
+        restrictions.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
+
+        EXAMPLES:
+
+        The canonical example comes from the symmetric matrices, which
+        decompose into diagonal and off-diagonal parts::
+
+            sage: J = RealSymmetricEJA(3)
+            sage: C = matrix(QQ, [ [1,0,0],
+            ....:                  [0,1,0],
+            ....:                  [0,0,0] ])
+            sage: c = J(C)
+            sage: J0,J5,J1 = J.peirce_decomposition(c)
+            sage: J0
+            Euclidean Jordan algebra of dimension 1...
+            sage: J5
+            Vector space of degree 6 and dimension 2...
+            sage: J1
+            Euclidean Jordan algebra of dimension 3...
+            sage: J0.one().natural_representation()
+            [0 0 0]
+            [0 0 0]
+            [0 0 1]
+            sage: orig_df = AA.options.display_format
+            sage: AA.options.display_format = 'radical'
+            sage: J.from_vector(J5.basis()[0]).natural_representation()
+            [          0           0 1/2*sqrt(2)]
+            [          0           0           0]
+            [1/2*sqrt(2)           0           0]
+            sage: J.from_vector(J5.basis()[1]).natural_representation()
+            [          0           0           0]
+            [          0           0 1/2*sqrt(2)]
+            [          0 1/2*sqrt(2)           0]
+            sage: AA.options.display_format = orig_df
+            sage: J1.one().natural_representation()
+            [1 0 0]
+            [0 1 0]
+            [0 0 0]
+
+        TESTS:
+
+        Every algebra decomposes trivially with respect to its identity
+        element::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: J0,J5,J1 = J.peirce_decomposition(J.one())
+            sage: J0.dimension() == 0 and J5.dimension() == 0
+            True
+            sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
+            True
+
+        The decomposition is into eigenspaces, and its components are
+        therefore necessarily orthogonal. Moreover, the identity
+        elements in the two subalgebras are the projections onto their
+        respective subspaces of the superalgebra's identity element::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: x = J.random_element()
+            sage: if not J.is_trivial():
+            ....:     while x.is_nilpotent():
+            ....:         x = J.random_element()
+            sage: c = x.subalgebra_idempotent()
+            sage: J0,J5,J1 = J.peirce_decomposition(c)
+            sage: ipsum = 0
+            sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
+            ....:     w = w.superalgebra_element()
+            ....:     y = J.from_vector(y)
+            ....:     z = z.superalgebra_element()
+            ....:     ipsum += w.inner_product(y).abs()
+            ....:     ipsum += w.inner_product(z).abs()
+            ....:     ipsum += y.inner_product(z).abs()
+            sage: ipsum
+            0
+            sage: J1(c) == J1.one()
+            True
+            sage: J0(J.one() - c) == J0.one()
+            True
+
+        """
+        if not c.is_idempotent():
+            raise ValueError("element is not idempotent: %s" % c)
+
+        # Default these to what they should be if they turn out to be
+        # trivial, because eigenspaces_left() won't return eigenvalues
+        # corresponding to trivial spaces (e.g. it returns only the
+        # eigenspace corresponding to lambda=1 if you take the
+        # decomposition relative to the identity element).
+        trivial = FiniteDimensionalEuclideanJordanSubalgebra(self, ())
+        J0 = trivial                          # eigenvalue zero
+        J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
+        J1 = trivial                          # eigenvalue one
+
+        for (eigval, eigspace) in c.operator().matrix().right_eigenspaces():
+            if eigval == ~(self.base_ring()(2)):
+                J5 = eigspace
+            else:
+                gens = tuple( self.from_vector(b) for b in eigspace.basis() )
+                subalg = FiniteDimensionalEuclideanJordanSubalgebra(self,
+                                                                    gens,
+                                                                    check_axioms=False)
+                if eigval == 0:
+                    J0 = subalg
+                elif eigval == 1:
+                    J1 = subalg
+                else:
+                    raise ValueError("unexpected eigenvalue: %s" % eigval)
+
+        return (J0, J5, J1)
+
+
+    def random_element(self, thorough=False):
+        r"""
+        Return a random element of this algebra.
+
+        Our algebra superclass method only returns a linear
+        combination of at most two basis elements. We instead
+        want the vector space "random element" method that
+        returns a more diverse selection.
+
+        INPUT:
+
+        - ``thorough`` -- (boolean; default False) whether or not we
+          should generate irrational coefficients for the random
+          element when our base ring is irrational; this slows the
+          algebra operations to a crawl, but any truly random method
+          should include them
+
+        """
+        # For a general base ring... maybe we can trust this to do the
+        # right thing? Unlikely, but.
+        V = self.vector_space()
+        v = V.random_element()
+
+        if self.base_ring() is AA:
+            # The "random element" method of the algebraic reals is
+            # stupid at the moment, and only returns integers between
+            # -2 and 2, inclusive:
+            #
+            #   https://trac.sagemath.org/ticket/30875
+            #
+            # Instead, we implement our own "random vector" method,
+            # and then coerce that into the algebra. We use the vector
+            # space degree here instead of the dimension because a
+            # subalgebra could (for example) be spanned by only two
+            # vectors, each with five coordinates.  We need to
+            # generate all five coordinates.
+            if thorough:
+                v *= QQbar.random_element().real()
+            else:
+                v *= QQ.random_element()

 
 

-    def random_elements(self, count):
+        return self.from_vector(V.coordinate_vector(v))
+
+    def random_elements(self, count, thorough=False):

         """
         Return ``count`` random elements as a tuple.
 
         """
         Return ``count`` random elements as a tuple.
 

+        INPUT:
+
+        - ``thorough`` -- (boolean; default False) whether or not we
+          should generate irrational coefficients for the random
+          elements when our base ring is irrational; this slows the
+          algebra operations to a crawl, but any truly random method
+          should include them
+

         SETUP::
 
             sage: from mjo.eja.eja_algebra import JordanSpinEJA
         SETUP::
 
             sage: from mjo.eja.eja_algebra import JordanSpinEJA


@@ -637,23 +828,85 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             True
 
         """
             True
 
         """

-        return  tuple( self.random_element() for idx in xrange(count) )
+        return tuple( self.random_element(thorough)
+                      for idx in range(count) )

 
 

+    @classmethod
+    def random_instance(cls, field=AA, **kwargs):
+        """
+        Return a random instance of this type of algebra.

 
 

-    def rank(self):
+        Beware, this will crash for "most instances" because the
+        constructor below looks wrong.

         """
         """

-        Return the rank of this EJA.
+        if cls is TrivialEJA:
+            # The TrivialEJA class doesn't take an "n" argument because
+            # there's only one.
+            return cls(field)
+
+        n = ZZ.random_element(cls._max_test_case_size() + 1)
+        return cls(n, field, **kwargs)
+
+    @cached_method
+    def _charpoly_coefficients(self):
+        r"""
+        The `r` polynomial coefficients of the "characteristic polynomial
+        of" function.
+        """
+        n = self.dimension()
+        var_names = [ "X" + str(z) for z in range(1,n+1) ]
+        R = PolynomialRing(self.base_ring(), var_names)
+        vars = R.gens()
+        F = R.fraction_field()
+
+        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(F, n, n, L_x_i_j)
+
+        r = None
+        if self.rank.is_in_cache():
+            r = self.rank()
+            # There's no need to pad the system with redundant
+            # columns if we *know* they'll be redundant.
+            n = r
+
+        # Compute an extra power in case the rank is equal to
+        # the dimension (otherwise, we would stop at x^(r-1)).
+        x_powers = [ (L_x**k)*self.one().to_vector()
+                     for k in range(n+1) ]
+        A = matrix.column(F, x_powers[:n])
+        AE = A.extended_echelon_form()
+        E = AE[:,n:]
+        A_rref = AE[:,:n]
+        if r is None:
+            r = A_rref.rank()
+        b = x_powers[r]
+
+        # The theory says that only the first "r" coefficients are
+        # nonzero, and they actually live in the original polynomial
+        # ring and not the fraction field. We negate them because
+        # in the actual characteristic polynomial, they get moved
+        # to the other side where x^r lives.
+        return -A_rref.solve_right(E*b).change_ring(R)[:r]

 
 

-        ALGORITHM:
+    @cached_method
+    def rank(self):
+        r"""
+        Return the rank of this EJA.

 
 

-        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.
+        This is a cached method because we know the rank a priori for
+        all of the algebras we can construct. Thus we can avoid the
+        expensive ``_charpoly_coefficients()`` call unless we truly
+        need to compute the whole characteristic polynomial.

 
         SETUP::
 
 
         SETUP::
 

-            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            ....:                                  JordanSpinEJA,

             ....:                                  RealSymmetricEJA,
             ....:                                  ComplexHermitianEJA,
             ....:                                  QuaternionHermitianEJA,
             ....:                                  RealSymmetricEJA,
             ....:                                  ComplexHermitianEJA,
             ....:                                  QuaternionHermitianEJA,


@@ -683,15 +936,54 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         TESTS:
 
         Ensure that every EJA that we know how to construct has a
         TESTS:
 
         Ensure that every EJA that we know how to construct has a

-        positive integer rank::
+        positive integer rank, unless the algebra is trivial in
+        which case its rank will be zero::

 
             sage: set_random_seed()
 
             sage: set_random_seed()

-            sage: r = random_eja().rank()
-            sage: r in ZZ and r > 0
+            sage: J = random_eja()
+            sage: r = J.rank()
+            sage: r in ZZ

             True
             True

+            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
+        return len(self._charpoly_coefficients())

 
 
     def vector_space(self):
 
 
     def vector_space(self):


@@ -715,155 +1007,11 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
     Element = FiniteDimensionalEuclideanJordanAlgebraElement
 
 
     Element = FiniteDimensionalEuclideanJordanAlgebraElement
 
 

-class KnownRankEJA(object):
-    """
-    A class for algebras that we actually know we can construct.  The
-    main issue is that, for most of our methods to make sense, we need
-    to know the rank of our algebra. Thus we can't simply generate a
-    "random" algebra, or even check that a given basis and product
-    satisfy the axioms; because even if everything looks OK, we wouldn't
-    know the rank we need to actuallty build the thing.
-
-    Not really a subclass of FDEJA because doing that causes method
-    resolution errors, e.g.
-
-      TypeError: Error when calling the metaclass bases
-      Cannot create a consistent method resolution
-      order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra,
-      KnownRankEJA

 
 

-    """
-    @staticmethod
-    def _max_test_case_size():
-        """
-        Return an integer "size" that is an upper bound on the size of
-        this algebra when it is used in a random test
-        case. Unfortunately, the term "size" is quite vague -- when
-        dealing with `R^n` under either the Hadamard or Jordan spin
-        product, the "size" refers to the dimension `n`. When dealing
-        with a matrix algebra (real symmetric or complex/quaternion
-        Hermitian), it refers to the size of the matrix, which is
-        far less than the dimension of the underlying vector space.
-
-        We default to five in this class, which is safe in `R^n`. The
-        matrix algebra subclasses (or any class where the "size" is
-        interpreted to be far less than the dimension) should override
-        with a smaller number.
-        """
-        return 5
-
-    @classmethod
-    def random_instance(cls, field=QQ, **kwargs):
-        """
-        Return a random instance of this type of algebra.
-
-        Beware, this will crash for "most instances" because the
-        constructor below looks wrong.
-        """
-        n = ZZ.random_element(cls._max_test_case_size()) + 1
-        return cls(n, field, **kwargs)
-
-
-class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra,
-                              KnownRankEJA):
-    """
-    Return the Euclidean Jordan Algebra corresponding to the set
-    `R^n` under the Hadamard product.
-
-    Note: this is nothing more than the Cartesian product of ``n``
-    copies of the spin algebra. Once Cartesian product algebras
-    are implemented, this can go.
-
-    SETUP::
-
-        sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
-
-    EXAMPLES:
-
-    This multiplication table can be verified by hand::
-
-        sage: J = RealCartesianProductEJA(3)
-        sage: e0,e1,e2 = J.gens()
-        sage: e0*e0
-        e0
-        sage: e0*e1
-        0
-        sage: e0*e2
-        0
-        sage: e1*e1
-        e1
-        sage: e1*e2
-        0
-        sage: e2*e2
-        e2
-
-    TESTS:
-
-    We can change the generator prefix::
-
-        sage: RealCartesianProductEJA(3, prefix='r').gens()
-        (r0, r1, r2)
-
-    """
-    def __init__(self, n, field=QQ, **kwargs):
-        V = VectorSpace(field, n)
-        mult_table = [ [ V.gen(i)*(i == j) for j in xrange(n) ]
-                       for i in xrange(n) ]
-
-        fdeja = super(RealCartesianProductEJA, self)
-        return fdeja.__init__(field, mult_table, rank=n, **kwargs)
-
-    def inner_product(self, x, y):
-        """
-        Faster to reimplement than to use natural representations.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
-
-        TESTS:
-
-        Ensure that this is the usual inner product for the algebras
-        over `R^n`::
-
-            sage: set_random_seed()
-            sage: J = RealCartesianProductEJA.random_instance()
-            sage: x,y = J.random_elements(2)
-            sage: X = x.natural_representation()
-            sage: Y = y.natural_representation()
-            sage: x.inner_product(y) == J.natural_inner_product(X,Y)
-            True
-
-        """
-        return x.to_vector().inner_product(y.to_vector())
-
-
-def random_eja():
+def random_eja(field=AA):

     """
     Return a "random" finite-dimensional Euclidean Jordan Algebra.
 
     """
     Return a "random" finite-dimensional Euclidean Jordan Algebra.
 

-    ALGORITHM:
-
-    For now, we choose a random natural number ``n`` (greater than zero)
-    and then give you back one of the following:
-
-      * The cartesian product of the rational numbers ``n`` times; this is
-        ``QQ^n`` with the Hadamard product.
-
-      * The Jordan spin algebra on ``QQ^n``.
-
-      * The ``n``-by-``n`` rational symmetric matrices with the symmetric
-        product.
-
-      * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
-        in the space of ``2n``-by-``2n`` real symmetric matrices.
-
-      * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
-        in the space of ``4n``-by-``4n`` real symmetric matrices.
-
-    Later this might be extended to return Cartesian products of the
-    EJAs above.
-

     SETUP::
 
         sage: from mjo.eja.eja_algebra import random_eja
     SETUP::
 
         sage: from mjo.eja.eja_algebra import random_eja


@@ -874,12 +1022,53 @@ def random_eja():
         Euclidean Jordan algebra of dimension...
 
     """
         Euclidean Jordan algebra of dimension...
 
     """

-    classname = choice(KnownRankEJA.__subclasses__())
-    return classname.random_instance()
+    classname = choice([TrivialEJA,
+                        HadamardEJA,
+                        JordanSpinEJA,
+                        RealSymmetricEJA,
+                        ComplexHermitianEJA,
+                        QuaternionHermitianEJA])
+    return classname.random_instance(field=field)
+
+

 
 
 
 

+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.
+        """
+        if self.base_ring() is QQ:
+            # There's no need to construct *another* algebra over the
+            # rationals if this one is already over the rationals.
+            superclass = super(RationalBasisEuclideanJordanAlgebra, self)
+            return superclass._charpoly_coefficients()
+
+        mult_table = 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)
+        return J._charpoly_coefficients()

 
 
 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
 
 
 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):


@@ -889,20 +1078,20 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         # field can have dimension 4 (quaternions) too.
         return 2
 
         # 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
         of the former when the product is known (like it is here).
         """
         """
         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_coeff() override.
+        # 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)
 
         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)
 

-        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.
             # 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.


@@ -910,54 +1099,64 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
             z = R.gen()
             p = z**2 - 2
             if p.is_irreducible():
             z = R.gen()
             p = z**2 - 2
             if p.is_irreducible():

-                field = NumberField(p, 'sqrt2', embedding=RLF(2).sqrt())
+                field = field.extension(p, 'sqrt2', embedding=RLF(2).sqrt())

                 basis = tuple( s.change_ring(field) for s in basis )
             self._basis_normalizers = tuple(
                 ~(self.natural_inner_product(s,s).sqrt()) for s in basis )
                 basis = tuple( s.change_ring(field) for s in basis )
             self._basis_normalizers = tuple(
                 ~(self.natural_inner_product(s,s).sqrt()) for s in basis )

-            basis = tuple(s*c for (s,c) in izip(basis,self._basis_normalizers))
+            basis = tuple(s*c for (s,c) in zip(basis,self._basis_normalizers))

 
         Qs = self.multiplication_table_from_matrix_basis(basis)
 
 
         Qs = self.multiplication_table_from_matrix_basis(basis)
 

-        fdeja = super(MatrixEuclideanJordanAlgebra, self)
-        return fdeja.__init__(field,
-                              Qs,
-                              rank=rank,
-                              natural_basis=basis,
-                              **kwargs)
+        super(MatrixEuclideanJordanAlgebra, self).__init__(field,
+                                                           Qs,
+                                                           natural_basis=basis,
+                                                           **kwargs)

 
 
     @cached_method
 
 
     @cached_method

-    def _charpoly_coeff(self, i):
-        """
+    def _charpoly_coefficients(self):
+        r"""

         Override the parent method with something that tries to compute
         over a faster (non-extension) field.
         """
         if self._basis_normalizers is None:
             # We didn't normalize, so assume that the basis we started
             # with had entries in a nice field.
         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)._charpoly_coeff(i)
+            return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coefficients()

         else:
         else:

-            basis = ( (b/n) for (b,n) in izip(self.natural_basis(),
-                                              self._basis_normalizers) )
+            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.
 
             # 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,
             J = MatrixEuclideanJordanAlgebra(QQ,
                                              basis,

-                                             self.rank(),
-                                             normalize_basis=False)
-            (_,x,_,_) = J._charpoly_matrix_system()
-            p = J._charpoly_coeff(i)
-            # p might be missing some vars, have to substitute "optionally"
-            pairs = izip(x.base_ring().gens(), self._basis_normalizers)
-            substitutions = { v: v*c for (v,c) in pairs }
-            result = p.subs(substitutions)
-
-            # The result of "subs" can be either a coefficient-ring
-            # element or a polynomial. Gotta handle both cases.
-            if result in QQ:
-                return self.base_ring()(result)
-            else:
-                return result.change_ring(self.base_ring())
+                                             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 )

 
 
     @staticmethod
 
 
     @staticmethod


@@ -975,15 +1174,18 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         # is supposed to hold the entire long vector, and the subspace W
         # of V will be spanned by the vectors that arise from symmetric
         # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
         # is supposed to hold the entire long vector, and the subspace W
         # of V will be spanned by the vectors that arise from symmetric
         # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.

+        if len(basis) == 0:
+            return []
+

         field = basis[0].base_ring()
         dimension = basis[0].nrows()
 
         V = VectorSpace(field, dimension**2)
         W = V.span_of_basis( _mat2vec(s) for s in basis )
         n = len(basis)
         field = basis[0].base_ring()
         dimension = basis[0].nrows()
 
         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))
 
                 mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
                 mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
 


@@ -1020,16 +1222,12 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         Xu = cls.real_unembed(X)
         Yu = cls.real_unembed(Y)
         tr = (Xu*Yu).trace()
         Xu = cls.real_unembed(X)
         Yu = cls.real_unembed(Y)
         tr = (Xu*Yu).trace()

-        if tr in RLF:
-            # It's real already.
-            return tr

 
 

-        # Otherwise, try the thing that works for complex numbers; and
-        # if that doesn't work, the thing that works for quaternions.

         try:
         try:

-            return tr.vector()[0] # real part, imag part is index 1
+            # Works in QQ, AA, RDF, et cetera.
+            return tr.real()

         except AttributeError:
         except AttributeError:

-            # A quaternions doesn't have a vector() method, but does
+            # A quaternion doesn't have a real() method, but does

             # have coefficient_tuple() method that returns the
             # coefficients of 1, i, j, and k -- in that order.
             return tr.coefficient_tuple()[0]
             # have coefficient_tuple() method that returns the
             # coefficients of 1, i, j, and k -- in that order.
             return tr.coefficient_tuple()[0]


@@ -1053,7 +1251,7 @@ class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         return M
 
 
         return M
 
 

-class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
+class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra):

     """
     The rank-n simple EJA consisting of real symmetric n-by-n
     matrices, the usual symmetric Jordan product, and the trace inner
     """
     The rank-n simple EJA consisting of real symmetric n-by-n
     matrices, the usual symmetric Jordan product, and the trace inner


@@ -1074,6 +1272,14 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
         sage: e2*e2
         e2
 
         sage: e2*e2
         e2
 

+    In theory, our "field" can be any subfield of the reals::
+
+        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
+

     TESTS:
 
     The dimension of this algebra is `(n^2 + n) / 2`::
     TESTS:
 
     The dimension of this algebra is `(n^2 + n) / 2`::


@@ -1123,6 +1329,11 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
         sage: x.operator().matrix().is_symmetric()
         True
 
         sage: x.operator().matrix().is_symmetric()
         True
 

+    We can construct the (trivial) algebra of rank zero::
+
+        sage: RealSymmetricEJA(0)
+        Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
+

     """
     @classmethod
     def _denormalized_basis(cls, n, field):
     """
     @classmethod
     def _denormalized_basis(cls, n, field):


@@ -1145,8 +1356,8 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
         # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
         # coordinates.
         S = []
         # 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
                 Eij = matrix(field, n, lambda k,l: k==i and l==j)
                 if i == j:
                     Sij = Eij


@@ -1161,9 +1372,13 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
         return 4 # Dimension 10
 
 
         return 4 # Dimension 10
 
 

-    def __init__(self, n, field=QQ, **kwargs):
+    def __init__(self, n, field=AA, **kwargs):

         basis = self._denormalized_basis(n, field)
         basis = self._denormalized_basis(n, field)

-        super(RealSymmetricEJA, self).__init__(field, basis, n, **kwargs)
+        super(RealSymmetricEJA, self).__init__(field,
+                                               basis,
+                                               check_axioms=False,
+                                               **kwargs)
+        self.rank.set_cache(n)

 
 
 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
 
 
 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):


@@ -1181,7 +1396,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
 
         EXAMPLES::
 
 
         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)
             sage: x1 = F(4 - 2*i)
             sage: x2 = F(1 + 2*i)
             sage: x3 = F(-i)


@@ -1201,7 +1416,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             sage: set_random_seed()
             sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
             sage: n = ZZ.random_element(n_max)
             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)
             sage: X = random_matrix(F, n)
             sage: Y = random_matrix(F, n)
             sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)


@@ -1214,15 +1429,17 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         n = M.nrows()
         if M.ncols() != n:
             raise ValueError("the matrix 'M' must be square")
         n = M.nrows()
         if M.ncols() != n:
             raise ValueError("the matrix 'M' must be square")

-        field = M.base_ring()
+
+        # We don't need any adjoined elements...
+        field = M.base_ring().base_ring()
+

         blocks = []
         for z in M.list():
         blocks = []
         for z in M.list():

-            a = z.vector()[0] # real part, I guess
-            b = z.vector()[1] # imag part, I guess
+            a = z.list()[0] # real part, I guess
+            b = z.list()[1] # imag part, I guess

             blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
 
             blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
 

-        # We can drop the imaginaries here.
-        return matrix.block(field.base_ring(), n, blocks)
+        return matrix.block(field, n, blocks)

 
 
     @staticmethod
 
 
     @staticmethod


@@ -1242,15 +1459,15 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             ....:                 [ 9,  10, 11, 12],
             ....:                 [-10, 9, -12, 11] ])
             sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
             ....:                 [ 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()
 
         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
             sage: M = random_matrix(F, 3)
             sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
             sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M


@@ -1263,17 +1480,24 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         if not n.mod(2).is_zero():
             raise ValueError("the matrix 'M' must be a complex embedding")
 
         if not n.mod(2).is_zero():
             raise ValueError("the matrix 'M' must be a complex embedding")
 

-        field = QQ
+        # 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()
         R = PolynomialRing(field, 'z')
         z = R.gen()

-        F = NumberField(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 = []
         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')
                 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')


@@ -1307,7 +1531,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             sage: Ye = y.natural_representation()
             sage: X = ComplexHermitianEJA.real_unembed(Xe)
             sage: Y = ComplexHermitianEJA.real_unembed(Ye)
             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
             sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
             sage: actual == expected
             True


@@ -1316,7 +1540,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/2
 
 
         return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/2
 
 

-class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
+class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra):

     """
     The rank-n simple EJA consisting of complex Hermitian n-by-n
     matrices over the real numbers, the usual symmetric Jordan product,
     """
     The rank-n simple EJA consisting of complex Hermitian n-by-n
     matrices over the real numbers, the usual symmetric Jordan product,


@@ -1327,6 +1551,16 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
 
         sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
 
 
         sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
 

+    EXAMPLES:
+
+    In theory, our "field" can be any subfield of the reals::
+
+        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
+

     TESTS:
 
     The dimension of this algebra is `n^2`::
     TESTS:
 
     The dimension of this algebra is `n^2`::


@@ -1376,6 +1610,11 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
         sage: x.operator().matrix().is_symmetric()
         True
 
         sage: x.operator().matrix().is_symmetric()
         True
 

+    We can construct the (trivial) algebra of rank zero::
+
+        sage: ComplexHermitianEJA(0)
+        Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
+

     """
 
     @classmethod
     """
 
     @classmethod


@@ -1403,9 +1642,9 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
             True
 
         """
             True
 
         """

-        R = PolynomialRing(QQ, 'z')
+        R = PolynomialRing(field, 'z')

         z = R.gen()
         z = R.gen()

-        F = NumberField(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
+        F = field.extension(z**2 + 1, 'I')

         I = F.gen()
 
         # This is like the symmetric case, but we need to be careful:
         I = F.gen()
 
         # This is like the symmetric case, but we need to be careful:


@@ -1414,8 +1653,8 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
         #   * The diagonal will (as a result) be real.
         #
         S = []
         #   * 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)
                 Eij = matrix(F, n, lambda k,l: k==i and l==j)
                 if i == j:
                     Sij = cls.real_embed(Eij)


@@ -1432,9 +1671,13 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
         return ( s.change_ring(field) for s in S )
 
 
         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)
         basis = self._denormalized_basis(n,field)

-        super(ComplexHermitianEJA,self).__init__(field, basis, n, **kwargs)
+        super(ComplexHermitianEJA,self).__init__(field,
+                                                 basis,
+                                                 check_axioms=False,
+                                                 **kwargs)
+        self.rank.set_cache(n)

 
 
 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
 
 
 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):


@@ -1484,7 +1727,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         if M.ncols() != n:
             raise ValueError("the matrix 'M' must be square")
 
         if M.ncols() != n:
             raise ValueError("the matrix 'M' must be square")
 

-        F = QuadraticField(-1, 'i')
+        F = QuadraticField(-1, 'I')

         i = F.gen()
 
         blocks = []
         i = F.gen()
 
         blocks = []


@@ -1540,7 +1783,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         if M.ncols() != n:
             raise ValueError("the matrix 'M' must be square")
         if not n.mod(4).is_zero():
         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 complex embedding")
+            raise ValueError("the matrix 'M' must be a quaternion embedding")

 
         # Use the base ring of the matrix to ensure that its entries can be
         # multiplied by elements of the quaternion algebra.
 
         # Use the base ring of the matrix to ensure that its entries can be
         # multiplied by elements of the quaternion algebra.


@@ -1552,18 +1795,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 = []
         # 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')
                 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)
                 elements.append(z)
 
         return matrix(Q, n/4, elements)


@@ -1600,8 +1843,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/4
 
 
         return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/4
 
 

-class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
-                             KnownRankEJA):
+class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra):

     """
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the
     """
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the


@@ -1612,6 +1854,16 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
 
         sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
 
 
         sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
 

+    EXAMPLES:
+
+    In theory, our "field" can be any subfield of the reals::
+
+        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
+

     TESTS:
 
     The dimension of this algebra is `2*n^2 - n`::
     TESTS:
 
     The dimension of this algebra is `2*n^2 - n`::


@@ -1661,6 +1913,11 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
         sage: x.operator().matrix().is_symmetric()
         True
 
         sage: x.operator().matrix().is_symmetric()
         True
 

+    We can construct the (trivial) algebra of rank zero::
+
+        sage: QuaternionHermitianEJA(0)
+        Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
+

     """
     @classmethod
     def _denormalized_basis(cls, n, field):
     """
     @classmethod
     def _denormalized_basis(cls, n, field):


@@ -1695,8 +1952,8 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
         #   * The diagonal will (as a result) be real.
         #
         S = []
         #   * 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)
                 Eij = matrix(Q, n, lambda k,l: k==i and l==j)
                 if i == j:
                     Sij = cls.real_embed(Eij)


@@ -1718,84 +1975,252 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
         return ( s.change_ring(field) for s in S )
 
 
         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)
         basis = self._denormalized_basis(n,field)

-        super(QuaternionHermitianEJA,self).__init__(field, basis, n, **kwargs)
+        super(QuaternionHermitianEJA,self).__init__(field,
+                                                    basis,
+                                                    check_axioms=False,
+                                                    **kwargs)
+        self.rank.set_cache(n)

 
 
 
 

-class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
+class HadamardEJA(RationalBasisEuclideanJordanAlgebra):

     """
     """

-    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
-    the reals.
+    Return the Euclidean Jordan Algebra corresponding to the set
+    `R^n` under the Hadamard product.
+
+    Note: this is nothing more than the Cartesian product of ``n``
+    copies of the spin algebra. Once Cartesian product algebras
+    are implemented, this can go.

 
     SETUP::
 
 
     SETUP::
 

-        sage: from mjo.eja.eja_algebra import JordanSpinEJA
+        sage: from mjo.eja.eja_algebra import HadamardEJA

 
     EXAMPLES:
 
     This multiplication table can be verified by hand::
 
 
     EXAMPLES:
 
     This multiplication table can be verified by hand::
 

-        sage: J = JordanSpinEJA(4)
-        sage: e0,e1,e2,e3 = J.gens()
+        sage: J = HadamardEJA(3)
+        sage: e0,e1,e2 = J.gens()

         sage: e0*e0
         e0
         sage: e0*e1
         sage: e0*e0
         e0
         sage: e0*e1

-        e1
-        sage: e0*e2
-        e2
-        sage: e0*e3
-        e3
-        sage: e1*e2

         0
         0

-        sage: e1*e3
+        sage: e0*e2

         0
         0

-        sage: e2*e3
+        sage: e1*e1
+        e1
+        sage: e1*e2

         0
         0

+        sage: e2*e2
+        e2
+
+    TESTS:

 
     We can change the generator prefix::
 
 
     We can change the generator prefix::
 

-        sage: JordanSpinEJA(2, prefix='B').gens()
-        (B0, B1)
+        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)
         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):
+        mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
+                       for i in range(n) ]
+
+        super(HadamardEJA, self).__init__(field,
+                                          mult_table,
+                                          check_axioms=False,
+                                          **kwargs)
+        self.rank.set_cache(n)
+
+    def inner_product(self, x, y):
+        """
+        Faster to reimplement than to use natural representations.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import HadamardEJA
+
+        TESTS:
+
+        Ensure that this is the usual inner product for the algebras
+        over `R^n`::
+
+            sage: set_random_seed()
+            sage: J = HadamardEJA.random_instance()
+            sage: x,y = J.random_elements(2)
+            sage: X = x.natural_representation()
+            sage: Y = y.natural_representation()
+            sage: x.inner_product(y) == J.natural_inner_product(X,Y)
+            True
+
+        """
+        return x.to_vector().inner_product(y.to_vector())
+
+
+class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
+    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::
+
+        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:]
                 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)
+                z0 = x0*y0 + (self._B*xbar).inner_product(ybar)

                 zbar = y0*xbar + x0*ybar
                 z = V([z0] + zbar.list())
                 mult_table[i][j] = z
 
                 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)
+        # 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).
+        super(BilinearFormEJA, self).__init__(field,
+                                              mult_table,
+                                              check_axioms=False,
+                                              **kwargs)
+        self.rank.set_cache(min(n,2))

 
     def inner_product(self, x, y):
 
     def inner_product(self, x, y):

-        """
-        Faster to reimplement than to use natural representations.
+        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::
 
 
         SETUP::
 

-            sage: from mjo.eja.eja_algebra import JordanSpinEJA
+            sage: from mjo.eja.eja_algebra import BilinearFormEJA

 
         TESTS:
 
 
         TESTS:
 

-        Ensure that this is the usual inner product for the algebras
-        over `R^n`::
+        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,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
+    the reals.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import JordanSpinEJA
+
+    EXAMPLES:
+
+    This multiplication table can be verified by hand::
+
+        sage: J = JordanSpinEJA(4)
+        sage: e0,e1,e2,e3 = J.gens()
+        sage: e0*e0
+        e0
+        sage: e0*e1
+        e1
+        sage: e0*e2
+        e2
+        sage: e0*e3
+        e3
+        sage: e1*e2
+        0
+        sage: e1*e3
+        0
+        sage: e2*e3
+        0
+
+    We can change the generator prefix::
+
+        sage: JordanSpinEJA(2, prefix='B').gens()
+        (B0, B1)
+
+    TESTS:
+
+        Ensure that we have the usual inner product on `R^n`::

 
             sage: set_random_seed()
             sage: J = JordanSpinEJA.random_instance()
 
             sage: set_random_seed()
             sage: J = JordanSpinEJA.random_instance()


@@ -1805,5 +2230,97 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
             sage: x.inner_product(y) == J.natural_inner_product(X,Y)
             True
 
             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.
+        super(JordanSpinEJA, self).__init__(n, field, **kwargs)
+
+
+class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra):
+    """
+    The trivial Euclidean Jordan algebra consisting of only a zero element.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import TrivialEJA
+
+    EXAMPLES::
+
+        sage: J = TrivialEJA()
+        sage: J.dimension()
+        0
+        sage: J.zero()
+        0
+        sage: J.one()
+        0
+        sage: 7*J.one()*12*J.one()
+        0
+        sage: J.one().inner_product(J.one())
+        0
+        sage: J.one().norm()
+        0
+        sage: J.one().subalgebra_generated_by()
+        Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
+        sage: J.rank()
+        0
+
+    """
+    def __init__(self, field=AA, **kwargs):
+        mult_table = []
+        super(TrivialEJA, self).__init__(field,
+                                         mult_table,
+                                         check_axioms=False,
+                                         **kwargs)
+        # The rank is zero using my definition, namely the dimension of the
+        # largest subalgebra generated by any element.
+        self.rank.set_cache(0)
+
+
+class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
+    r"""
+    The external (orthogonal) direct sum of two other Euclidean Jordan
+    algebras. Essentially the Cartesian product of its two factors.
+    Every Euclidean Jordan algebra decomposes into an orthogonal
+    direct sum of simple Euclidean Jordan algebras, so no generality
+    is lost by providing only this construction.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import (HadamardEJA,
+        ....:                                  RealSymmetricEJA,
+        ....:                                  DirectSumEJA)
+
+    EXAMPLES::
+
+        sage: J1 = HadamardEJA(2)
+        sage: J2 = RealSymmetricEJA(3)
+        sage: J = DirectSumEJA(J1,J2)
+        sage: J.dimension()
+        8
+        sage: J.rank()
+        5
+
+    """
+    def __init__(self, J1, J2, field=AA, **kwargs):
+        n1 = J1.dimension()
+        n2 = J2.dimension()
+        n = n1+n2
+        V = VectorSpace(field, n)
+        mult_table = [ [ V.zero() for j in range(n) ]
+                       for i in range(n) ]
+        for i in range(n1):
+            for j in range(n1):
+                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):
+                p = (J2.monomial(i)*J2.monomial(j)).to_vector()
+                mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
+
+        super(DirectSumEJA, self).__init__(field,
+                                           mult_table,
+                                           check_axioms=False,
+                                           **kwargs)
+        self.rank.set_cache(J1.rank() + J2.rank())