]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: make all tests work in trivial algebras.
[sage.d.git] / mjo / eja / eja_algebra.py
index 0a260653a7c5480acea6ffc00a9158169e4e7a18..055bbbda83d7576fda8f5e1794d759c2eee467bc 100644 (file)
@@ -54,7 +54,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
     def __init__(self,
                  field,
                  mult_table,
     def __init__(self,
                  field,
                  mult_table,
-                 rank,
                  prefix='e',
                  category=None,
                  natural_basis=None,
                  prefix='e',
                  category=None,
                  natural_basis=None,
@@ -91,7 +90,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
                 # a real embedding.
                 raise ValueError('field is not real')
 
                 # a real embedding.
                 raise ValueError('field is not real')
 
-        self._rank = rank
         self._natural_basis = natural_basis
 
         if category is None:
         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)
 
         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):
         """
@@ -219,163 +235,6 @@ 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():
-            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):
         """
     @cached_method
     def characteristic_polynomial(self):
         """
@@ -420,20 +279,21 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         r = self.rank()
         n = self.dimension()
 
         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.
 
         # 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)
+        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):
 
 
     def inner_product(self, x, y):
@@ -773,108 +633,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return (J0, J5, J1)
 
 
         return (J0, J5, J1)
 
 
-    def a_jordan_frame(self):
-        r"""
-        Generate a Jordan frame for this algebra.
-
-        This implementation is based on the so-called "central
-        orthogonal idempotents" implemented for (semisimple) centers
-        of SageMath ``FiniteDimensionalAlgebrasWithBasis``. Since all
-        Euclidean Jordan algebas are commutative (and thus equal to
-        their own centers) and semisimple, the method should work more
-        or less as implemented, if it ever worked in the first place.
-        (I don't know the justification for the original implementation.
-        yet).
-
-        How it works: we loop through the algebras generators, looking
-        for their eigenspaces. If there's more than one eigenspace,
-        and if they result in more than one subalgebra, then we split
-        those subalgebras recursively until we get to subalgebras of
-        dimension one (whose idempotent is the unit element). Why does
-        some generator have to produce at least two subalgebras? I
-        dunno. But it seems to work.
-
-        Beware that Koecher defines the "center" of a Jordan algebra to
-        be something else, because the usual definition is stupid in a
-        (necessarily commutative) Jordan algebra.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import (random_eja,
-            ....:                                  JordanSpinEJA,
-            ....:                                  TrivialEJA)
-
-        EXAMPLES:
-
-        A Jordan frame for the trivial algebra has to be empty
-        (zero-length) since its rank is zero. More to the point, there
-        are no non-zero idempotents in the trivial EJA. This does not
-        cause any problems so long as we adopt the convention that the
-        empty sum is zero, since then the sole element of the trivial
-        EJA has an (empty) spectral decomposition::
-
-            sage: J = TrivialEJA()
-            sage: J.a_jordan_frame()
-            ()
-
-        A one-dimensional algebra has rank one (equal to its dimension),
-        and only one primitive idempotent, namely the algebra's unit
-        element::
-
-            sage: J = JordanSpinEJA(1)
-            sage: J.a_jordan_frame()
-            (e0,)
-
-        TESTS::
-
-            sage: J = random_eja()
-            sage: c = J.a_jordan_frame()
-            sage: all( x^2 == x for x in c )
-            True
-            sage: r = len(c)
-            sage: all( c[i]*c[j] == c[i]*(i==j) for i in range(r)
-            ....:                               for j in range(r) )
-            True
-
-        """
-        if self.dimension() == 0:
-            return ()
-        if self.dimension() == 1:
-            return (self.one(),)
-
-        for g in self.gens():
-            eigenpairs = g.operator().matrix().right_eigenspaces()
-            if len(eigenpairs) >= 2:
-                subalgebras = []
-                for eigval, eigspace in eigenpairs:
-                    # Make sub-EJAs from the matrix eigenspaces...
-                    sb = tuple( self.from_vector(b) for b in eigspace.basis() )
-                    try:
-                        # This will fail if e.g. the eigenspace basis
-                        # contains two elements and their product
-                        # isn't a linear combination of the two of
-                        # them (i.e. the generated EJA isn't actually
-                        # two dimensional).
-                        s = FiniteDimensionalEuclideanJordanSubalgebra(self, sb)
-                        subalgebras.append(s)
-                    except ArithmeticError as e:
-                        if str(e) == "vector is not in free module":
-                            # Ignore only the "not a sub-EJA" error
-                            pass
-
-                if len(subalgebras) >= 2:
-                    # apply this method recursively.
-                    return tuple( c.superalgebra_element()
-                                  for subalgebra in subalgebras
-                                  for c in subalgebra.a_jordan_frame() )
-
-        # If we got here, the algebra didn't decompose, at least not when we looked at
-        # the eigenspaces corresponding only to basis elements of the algebra. The
-        # implementation I stole says that this should work because of Schur's Lemma,
-        # so I personally blame Schur's Lemma if it does not.
-        raise Exception("Schur's Lemma didn't work!")
-
-
     def random_elements(self, count):
         """
         Return ``count`` random elements as a tuple.
     def random_elements(self, count):
         """
         Return ``count`` random elements as a tuple.
@@ -893,23 +651,84 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             True
 
         """
             True
 
         """
-        return  tuple( self.random_element() for idx in range(count) )
+        return tuple( self.random_element() 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)
 
 
-        ALGORITHM:
+        n = ZZ.random_element(cls._max_test_case_size()) + 1
+        return cls(n, field, **kwargs)
 
 
-        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.
+    @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]
+
+    @cached_method
+    def rank(self):
+        r"""
+        Return the rank of this EJA.
+
+        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,
@@ -950,8 +769,43 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: r > 0 or (r == 0 and J.is_trivial())
             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):
@@ -975,61 +829,7 @@ 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=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.
     """
     Return the Euclidean Jordan Algebra corresponding to the set
     `R^n` under the Hadamard product.
@@ -1075,7 +875,8 @@ class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
                        for i in range(n) ]
 
         fdeja = super(HadamardEJA, self)
                        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):
         """
 
     def inner_product(self, x, y):
         """
@@ -1102,7 +903,7 @@ class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
         return x.to_vector().inner_product(y.to_vector())
 
 
         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.
 
     """
     Return a "random" finite-dimensional Euclidean Jordan Algebra.
 
@@ -1116,17 +917,17 @@ def random_eja(field=AA, nontrivial=False):
         Euclidean Jordan algebra of dimension...
 
     """
         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)
 
 
 
 
     return classname.random_instance(field=field)
 
 
 
 
-
-
 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
     @staticmethod
     def _max_test_case_size():
 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
     @staticmethod
     def _max_test_case_size():
@@ -1134,20 +935,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.
@@ -1164,45 +965,50 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         Qs = self.multiplication_table_from_matrix_basis(basis)
 
         fdeja = super(MatrixEuclideanJordanAlgebra, self)
         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
 
 
     @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:
             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.
         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,
             J = MatrixEuclideanJordanAlgebra(QQ,
                                              basis,
-                                             self.rank(),
                                              normalize_basis=False)
                                              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
 
 
     @staticmethod
@@ -1220,6 +1026,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.
         # 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()
 
         field = basis[0].base_ring()
         dimension = basis[0].nrows()
 
@@ -1299,7 +1108,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
@@ -1377,6 +1186,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):
@@ -1417,7 +1231,8 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
 
     def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n, field)
 
     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):
 
 
 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
@@ -1579,7 +1394,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,
@@ -1649,6 +1464,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
@@ -1707,7 +1527,8 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
 
     def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n,field)
 
     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):
 
 
 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
@@ -1873,8 +1694,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
@@ -1944,6 +1764,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):
@@ -2003,10 +1828,11 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
 
     def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n,field)
 
     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 =
     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 =
@@ -2086,7 +1912,8 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
         # one-dimensional ambient space (because the rank is bounded
         # by the ambient dimension).
         fdeja = super(BilinearFormEJA, self)
         # 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"""
 
     def inner_product(self, x, y):
         r"""
@@ -2180,7 +2007,7 @@ class JordanSpinEJA(BilinearFormEJA):
         return super(JordanSpinEJA, self).__init__(n, field, **kwargs)
 
 
         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.
 
     """
     The trivial Euclidean Jordan algebra consisting of only a zero element.
 
@@ -2214,4 +2041,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.
         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)