]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: replace rank computation with length of charpoly coefficients.
[sage.d.git] / mjo / eja / eja_algebra.py
index 19db8b0ccefef297623d81721d4a774f99a9d3d1..d1c0fa2eafd5d3d5e17ccab175e16ee44c3e9d82 100644 (file)
@@ -192,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):
         """
@@ -791,9 +809,62 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         """
         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.
+
+        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)
+
     @cached_method
-    def rank(self):
+    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)
+        # 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]
+        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:
@@ -886,43 +957,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             2
 
         """
-        n = self.dimension()
-        if n == 0:
-            return 0
-        elif n == 1:
-            return 1
-
-        var_names = [ "X" + str(z) for z in range(1,n+1) ]
-        R = PolynomialRing(self.base_ring(), var_names)
-        vars = R.gens()
-
-        def L_x_i_j(i,j):
-            # From a result in my book, these are the entries of the
-            # basis representation of L_x.
-            return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
-                        for k in range(n) )
-
-        L_x = matrix(R, n, n, L_x_i_j)
-        x_powers = [ vars[k]*(L_x**k)*self.one().to_vector()
-                     for k in range(n) ]
-
-        # Can assume n >= 2
-        M = matrix([x_powers[0]])
-        old_rank = 1
-
-        for d in range(1,n):
-            M = matrix(M.rows() + [x_powers[d]])
-            M.echelonize()
-            # TODO: we've basically solved the system here.
-            # We should save the echelonized matrix somehow
-            # so that it can be reused in the charpoly method.
-            new_rank = M.rank()
-            if new_rank == old_rank:
-                return new_rank
-            else:
-                old_rank = new_rank
-
-        return n
+        return len(self._charpoly_coefficients())
 
 
     def vector_space(self):
@@ -946,61 +981,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.
@@ -1088,9 +1069,13 @@ def random_eja(field=AA, nontrivial=False):
         Euclidean Jordan algebra of dimension...
 
     """
-    eja_classes = KnownRankEJA.__subclasses__()
-    if nontrivial:
-        eja_classes.remove(TrivialEJA)
+    eja_classes = [HadamardEJA,
+                   JordanSpinEJA,
+                   RealSymmetricEJA,
+                   ComplexHermitianEJA,
+                   QuaternionHermitianEJA]
+    if not nontrivial:
+        eja_classes.append(TrivialEJA)
     classname = choice(eja_classes)
     return classname.random_instance(field=field)
 
@@ -1289,7 +1274,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
@@ -1570,7 +1555,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,
@@ -1865,8 +1850,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
@@ -1999,7 +1983,7 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
         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 =
@@ -2174,7 +2158,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.