]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: add "of" to the algebra characteristic_polynomial() method name.
[sage.d.git] / mjo / eja / eja_algebra.py
index 34fc0c3bc7b1dfe250fd5f7ea5de4fbbc52c66e0..9aa40eeacdacbb54c792158b1462d2f36d3dc4d6 100644 (file)
@@ -54,7 +54,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
     def __init__(self,
                  field,
                  mult_table,
-                 rank,
                  prefix='e',
                  category=None,
                  natural_basis=None,
@@ -91,7 +90,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
                 # a real embedding.
                 raise ValueError('field is not real')
 
-        self._rank = rank
         self._natural_basis = natural_basis
 
         if category is None:
@@ -194,6 +192,24 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         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):
         """
@@ -219,168 +235,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
     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():
-            return R.one()
-        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())
-
-
     @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
@@ -400,7 +261,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         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)
@@ -413,27 +274,28 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         any argument::
 
             sage: J = TrivialEJA()
-            sage: J.characteristic_polynomial()
+            sage: J.characteristic_polynomial_of()
             1
 
         """
         r = self.rank()
         n = self.dimension()
 
-        # The list of coefficient polynomials a_0, a_1, a_2, ..., a_n.
-        a = [ self._charpoly_coeff(i) for i in range(r+1) ]
+        # 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.
-        R = a[0].parent()
         S = PolynomialRing(self.base_ring(),'t')
         t = S.gen(0)
-        S = PolynomialRing(S, R.variable_names())
-        t = S(t)
+        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 range(len(a)) )
+        return (t**r + sum( a[k]*(t**k) for k in range(r) ))
 
 
     def inner_product(self, x, y):
@@ -578,8 +440,15 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         """
         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()
@@ -791,125 +660,35 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             True
 
         """
-        return  tuple( self.random_element() for idx in range(count) )
-
-
-    def _rank_computation1(self):
-        r"""
-        Compute the rank of this algebra using highly suspicious voodoo.
-
-        ALGORITHM:
-
-        We first compute the basis representation of the operator L_x
-        using polynomial indeterminates are placeholders for the
-        coordinates of "x", which is arbitrary. We then use that
-        matrix to compute the (polynomial) entries of x^0, x^1, ...,
-        x^d,... for increasing values of "d", starting at zero. The
-        idea is that. If we also add "coefficient variables" a_0,
-        a_1,...  to the ring, we can form the linear combination
-        a_0*x^0 + ... + a_d*x^d = 0, and ask what dimension the
-        solution space has as an affine variety. When "d" is smaller
-        than the rank, we expect that dimension to be the number of
-        coordinates of "x", since we can set *those* to whatever we
-        want, but linear independence forces the coefficients a_i to
-        be zero. Eventually, when "d" passes the rank, the dimension
-        of the solution space begins to grow, because we can *still*
-        set the coordinates of "x" arbitrarily, but now there are some
-        coefficients that make the sum zero as well. So, when the
-        dimension of the variety jumps, we return the corresponding
-        "d" as the rank of the algebra. This appears to work.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import (HadamardEJA,
-            ....:                                  JordanSpinEJA,
-            ....:                                  RealSymmetricEJA,
-            ....:                                  ComplexHermitianEJA,
-            ....:                                  QuaternionHermitianEJA)
-
-        EXAMPLES::
-
-            sage: J = HadamardEJA(5)
-            sage: J._rank_computation() == J.rank()
-            True
-            sage: J = JordanSpinEJA(5)
-            sage: J._rank_computation() == J.rank()
-            True
-            sage: J = RealSymmetricEJA(4)
-            sage: J._rank_computation() == J.rank()
-            True
-            sage: J = ComplexHermitianEJA(3)
-            sage: J._rank_computation() == J.rank()
-            True
-            sage: J = QuaternionHermitianEJA(2)
-            sage: J._rank_computation() == J.rank()
-            True
+        return tuple( self.random_element() for idx in range(count) )
 
+    @classmethod
+    def random_instance(cls, field=AA, **kwargs):
         """
-        n = self.dimension()
-        var_names = [ "X" + str(z) for z in range(1,n+1) ]
-        d = 0
-        ideal_dim = len(var_names)
-        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[d+k]*self.monomial(k).operator().matrix()[i,j]
-                        for k in range(n) )
+        Return a random instance of this type of algebra.
 
-        while ideal_dim == len(var_names):
-            coeff_names = [ "a" + str(z) for z in range(d) ]
-            R = PolynomialRing(self.base_ring(), coeff_names + var_names)
-            vars = R.gens()
-            L_x = matrix(R, n, n, L_x_i_j)
-            x_powers = [ vars[k]*(L_x**k)*self.one().to_vector()
-                         for k in range(d) ]
-            eqs = [ sum(x_powers[k][j] for k in range(d)) for j in range(n) ]
-            ideal_dim = R.ideal(eqs).dimension()
-            d += 1
-
-        # Subtract one because we increment one too many times, and
-        # subtract another one because "d" is one greater than the
-        # answer anyway; when d=3, we go up to x^2.
-        return d-2
-
-    def _rank_computation2(self):
-        r"""
-        Instead of using the dimension of an ideal, find the rank of a
-        matrix containing polynomials.
+        Beware, this will crash for "most instances" because the
+        constructor below looks wrong.
         """
-        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()
-
-        def L_x_i_j(i,j):
-            # From a result in my book, these are the entries of the
-            # basis representation of L_x.
-            return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
-                        for k in range(n) )
-
-        L_x = matrix(R, n, n, L_x_i_j)
-        x_powers = [ (vars[k]*(L_x**k)*self.one().to_vector()).row()
-                     for k in range(n) ]
+        if cls is TrivialEJA:
+            # The TrivialEJA class doesn't take an "n" argument because
+            # there's only one.
+            return cls(field)
 
-        from sage.matrix.constructor import block_matrix
-        M = block_matrix(n,1,x_powers)
-        return M.rank()
+        n = ZZ.random_element(cls._max_test_case_size() + 1)
+        return cls(n, field, **kwargs)
 
-    def _rank_computation3(self):
+    @cached_method
+    def _charpoly_coefficients(self):
         r"""
-        Similar to :meth:`_rank_computation2`, but it stops echelonizing
-        as soon as it hits the first zero row.
+        The `r` polynomial coefficients of the "characteristic polynomial
+        of" function.
         """
         n = self.dimension()
-        if n == 0:
-            return 0
-        elif n == 1:
-            return 1
-
         var_names = [ "X" + str(z) for z in range(1,n+1) ]
         R = PolynomialRing(self.base_ring(), var_names)
         vars = R.gens()
+        F = R.fraction_field()
 
         def L_x_i_j(i,j):
             # From a result in my book, these are the entries of the
@@ -917,39 +696,48 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
                         for k in range(n) )
 
-        L_x = matrix(R, n, n, L_x_i_j)
-        x_powers = [ vars[k]*(L_x**k)*self.one().to_vector()
-                     for k in range(n) ]
-
-        # Can assume n >= 2
-        rows = [x_powers[0]]
-        M = matrix(rows)
-        old_rank = 1
-
-        for d in range(1,n):
-            rows = M.rows() + [x_powers[d]]
-            M = matrix(rows)
-            M.echelonize()
-            new_rank = M.rank()
-            if new_rank == old_rank:
-                return new_rank
-            else:
-                old_rank = new_rank
+        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]
 
+    @cached_method
     def rank(self):
-        """
+        r"""
         Return the rank of this EJA.
 
-        ALGORITHM:
-
-        The author knows of no algorithm to compute the rank of an EJA
-        where only the multiplication table is known. In lieu of one, we
-        require the rank to be specified when the algebra is created,
-        and simply pass along that number here.
+        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::
 
-            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            ....:                                  JordanSpinEJA,
             ....:                                  RealSymmetricEJA,
             ....:                                  ComplexHermitianEJA,
             ....:                                  QuaternionHermitianEJA,
@@ -990,8 +778,43 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: r > 0 or (r == 0 and J.is_trivial())
             True
 
+        Ensure that computing the rank actually works, since the ranks
+        of all simple algebras are known and will be cached by default::
+
+            sage: J = HadamardEJA(4)
+            sage: J.rank.clear_cache()
+            sage: J.rank()
+            4
+
+        ::
+
+            sage: J = JordanSpinEJA(4)
+            sage: J.rank.clear_cache()
+            sage: J.rank()
+            2
+
+        ::
+
+            sage: J = RealSymmetricEJA(3)
+            sage: J.rank.clear_cache()
+            sage: J.rank()
+            3
+
+        ::
+
+            sage: J = ComplexHermitianEJA(2)
+            sage: J.rank.clear_cache()
+            sage: J.rank()
+            2
+
+        ::
+
+            sage: J = QuaternionHermitianEJA(2)
+            sage: J.rank.clear_cache()
+            sage: J.rank()
+            2
         """
-        return self._rank
+        return len(self._charpoly_coefficients())
 
 
     def vector_space(self):
@@ -1015,61 +838,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
     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=AA, **kwargs):
-        """
-        Return a random instance of this type of algebra.
-
-        Beware, this will crash for "most instances" because the
-        constructor below looks wrong.
-        """
-        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)
-
-
-class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
+class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra):
     """
     Return the Euclidean Jordan Algebra corresponding to the set
     `R^n` under the Hadamard product.
@@ -1115,7 +884,8 @@ class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
                        for i in range(n) ]
 
         fdeja = super(HadamardEJA, self)
-        return fdeja.__init__(field, mult_table, rank=n, **kwargs)
+        fdeja.__init__(field, mult_table, **kwargs)
+        self.rank.set_cache(n)
 
     def inner_product(self, x, y):
         """
@@ -1142,7 +912,7 @@ class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
         return x.to_vector().inner_product(y.to_vector())
 
 
-def random_eja(field=AA, nontrivial=False):
+def random_eja(field=AA):
     """
     Return a "random" finite-dimensional Euclidean Jordan Algebra.
 
@@ -1156,17 +926,17 @@ def random_eja(field=AA, nontrivial=False):
         Euclidean Jordan algebra of dimension...
 
     """
-    eja_classes = KnownRankEJA.__subclasses__()
-    if nontrivial:
-        eja_classes.remove(TrivialEJA)
-    classname = choice(eja_classes)
+    classname = choice([TrivialEJA,
+                        HadamardEJA,
+                        JordanSpinEJA,
+                        RealSymmetricEJA,
+                        ComplexHermitianEJA,
+                        QuaternionHermitianEJA])
     return classname.random_instance(field=field)
 
 
 
 
-
-
 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
     @staticmethod
     def _max_test_case_size():
@@ -1174,20 +944,20 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         # field can have dimension 4 (quaternions) too.
         return 2
 
-    def __init__(self, field, basis, rank, normalize_basis=True, **kwargs):
+    def __init__(self, field, basis, normalize_basis=True, **kwargs):
         """
         Compared to the superclass constructor, we take a basis instead of
         a multiplication table because the latter can be computed in terms
         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)
 
-        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.
@@ -1204,45 +974,50 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         Qs = self.multiplication_table_from_matrix_basis(basis)
 
         fdeja = super(MatrixEuclideanJordanAlgebra, self)
-        return fdeja.__init__(field,
-                              Qs,
-                              rank=rank,
-                              natural_basis=basis,
-                              **kwargs)
+        fdeja.__init__(field, Qs, natural_basis=basis, **kwargs)
+        return
 
 
     @cached_method
-    def _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.
-            return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coeff(i)
+            return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coefficients()
         else:
             basis = ( (b/n) for (b,n) in zip(self.natural_basis(),
                                              self._basis_normalizers) )
 
             # Do this over the rationals and convert back at the end.
+            # Only works because we know the entries of the basis are
+            # integers.
             J = MatrixEuclideanJordanAlgebra(QQ,
                                              basis,
-                                             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 = zip(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())
+            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
@@ -1260,6 +1035,9 @@ 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.
+        if len(basis) == 0:
+            return []
+
         field = basis[0].base_ring()
         dimension = basis[0].nrows()
 
@@ -1339,7 +1117,7 @@ class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         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
@@ -1417,6 +1195,11 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
         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):
@@ -1457,7 +1240,8 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
 
     def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n, field)
-        super(RealSymmetricEJA, self).__init__(field, basis, n, **kwargs)
+        super(RealSymmetricEJA, self).__init__(field, basis, **kwargs)
+        self.rank.set_cache(n)
 
 
 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
@@ -1619,7 +1403,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         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,
@@ -1689,6 +1473,11 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
         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
@@ -1747,7 +1536,8 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
 
     def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n,field)
-        super(ComplexHermitianEJA,self).__init__(field, basis, n, **kwargs)
+        super(ComplexHermitianEJA,self).__init__(field, basis, **kwargs)
+        self.rank.set_cache(n)
 
 
 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
@@ -1913,8 +1703,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         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
@@ -1984,6 +1773,11 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
         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):
@@ -2043,10 +1837,11 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
 
     def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n,field)
-        super(QuaternionHermitianEJA,self).__init__(field, basis, n, **kwargs)
+        super(QuaternionHermitianEJA,self).__init__(field, basis, **kwargs)
+        self.rank.set_cache(n)
 
 
-class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
+class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra):
     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 =
@@ -2126,7 +1921,8 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
         # one-dimensional ambient space (because the rank is bounded
         # by the ambient dimension).
         fdeja = super(BilinearFormEJA, self)
-        return fdeja.__init__(field, mult_table, rank=min(n,2), **kwargs)
+        fdeja.__init__(field, mult_table, **kwargs)
+        self.rank.set_cache(min(n,2))
 
     def inner_product(self, x, y):
         r"""
@@ -2220,7 +2016,7 @@ class JordanSpinEJA(BilinearFormEJA):
         return super(JordanSpinEJA, self).__init__(n, field, **kwargs)
 
 
-class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
+class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra):
     """
     The trivial Euclidean Jordan algebra consisting of only a zero element.
 
@@ -2254,4 +2050,5 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
         fdeja = super(TrivialEJA, self)
         # The rank is zero using my definition, namely the dimension of the
         # largest subalgebra generated by any element.
-        return fdeja.__init__(field, mult_table, rank=0, **kwargs)
+        fdeja.__init__(field, mult_table, **kwargs)
+        self.rank.set_cache(0)