Ticket #8614: trac-8614-optimize-modular-symbol-relations-rebase.patch

module_list.py

@@ -1028 +1028 @@
               extra_compile_args=["-std=c99",  "-D_XPG6"],
               depends = flint_depends),
 
+    Extension('sage.modular.modsym.relation_matrix_pyx',
+              sources = ['sage/modular/modsym/relation_matrix_pyx.pyx']),
+
     Extension('sage.modular.modsym.heilbronn',
               sources = ['sage/modular/modsym/heilbronn.pyx'],
               libraries = ["csage", "flint", "gmp", "gmpxx", "m", "stdc++"],

sage/modular/modsym/manin_symbols.py

@@ -182 +182 @@
 
         """
         raise NotImplementedError, "Only implemented in derived classes"
+
+    def _apply_S_only_0pm1(self):
+        """
+        Return True if the coefficient when applying the S relation is
+        always 0, 1, or -1.  This is useful for optimizing code in
+        relation_matrix.py.
+
+        EXAMPLES::
+
+            sage: eps = DirichletGroup(4).gen(0)
+            sage: from sage.modular.modsym.manin_symbols import ManinSymbolList_character
+            sage: ManinSymbolList_character(eps,2)._apply_S_only_0pm1()
+            True
+            sage: eps = DirichletGroup(7).gen(0)
+            sage: from sage.modular.modsym.manin_symbols import ManinSymbolList_character
+            sage: ManinSymbolList_character(eps,2)._apply_S_only_0pm1()
+            False
+        """
+        return False # derived classes could overload and put True
         
     def apply_S(self, j):
         """
@@ -501 +520 @@
         else:
             return k, -1
         
+    def _apply_S_only_0pm1(self):
+        """
+        Return True if the coefficient when applying the S relation is
+        always 0, 1, or -1.  This is useful for optimizing code in
+        relation_matrix.py.
+
+        EXAMPLES::
+
+            sage: from sage.modular.modsym.manin_symbols import ManinSymbolList_gamma0
+            sage: ManinSymbolList_gamma0(5,8)._apply_S_only_0pm1()
+            True
+        """
+        return True
+
     def apply_I(self, j):
         """
         Apply the matrix `I=[-1,0,0,1]` to the `j`-th Manin symbol.
@@ -1047 +1080 @@
         else:
             return k, -s
         
+    def _apply_S_only_0pm1(self):
+        """
+        Return True if the coefficient when applying the S relation is
+        always 0, 1, or -1.  This is useful for optimizing code in
+        relation_matrix.py.
+
+        EXAMPLES::
+
+            sage: eps = DirichletGroup(4).gen(0)
+            sage: from sage.modular.modsym.manin_symbols import ManinSymbolList_character
+            sage: ManinSymbolList_character(eps,2)._apply_S_only_0pm1()
+            True
+            sage: ManinSymbolList_character(DirichletGroup(13).0,2)._apply_S_only_0pm1()
+            False
+        """
+        return self.__character.order() <= 2
+
     def apply_I(self, j):
         """
         Apply the matrix `I=[-1,0,0,1]` to the `j`-th Manin symbol.

sage/modular/modsym/p1list.pyx

@@ -678 +678 @@
             sage: L
             The projective line over the integers modulo 120
         """
+        cdef int i
+
         self.__N = N
         if N <= 46340:
             self.__list = p1list_int(N)
@@ -700 +702 @@
         if not self.t: raise MemoryError
 
         # Initialize xgcd table
-        cdef int i
         cdef llong ll_s, ll_t, ll_N = N
 
         if N <= 46340:

sage/modular/modsym/relation_matrix.py

@@ -459 +459 @@
     if sign != 0:
         # Let rels = rels union I relations.
         rels.update(modI_relations(syms,sign))
-    mod = sparse_2term_quotient(rels, len(syms), field)
+
+    if syms._apply_S_only_0pm1() and rings.is_RationalField(field):
+        import relation_matrix_pyx
+        mod = relation_matrix_pyx.sparse_2term_quotient_only_pm1(rels, len(syms))
+    else:
+        mod = sparse_2term_quotient(rels, len(syms), field)
+        
     R = T_relation_matrix_wtk_g0(syms, mod, field, sparse)
     return R, mod
     

new file sage/modular/modsym/relation_matrix_pyx.pyx

@@ +1 @@
+"""
+Optimized Cython code for computing relation matrices in certain cases.
+"""
+
+#############################################################################
+#       Copyright (C) 2010 William Stein <wstein@gmail.com>
+#  Distributed under the terms of the GNU General Public License (GPL) v2+.
+#  The full text of the GPL is available at:
+#                  http://www.gnu.org/licenses/
+#############################################################################
+
+import sage.misc.misc as misc
+from sage.rings.rational cimport Rational
+
+def sparse_2term_quotient_only_pm1(rels, n):
+    r"""
+    Performs Sparse Gauss elimination on a matrix all of whose columns
+    have at most 2 nonzero entries with relations all 1 or -1.
+
+    This algorithm is more subtle than just "identify symbols in pairs",
+    since complicated relations can cause generators to equal 0.
+
+    NOTE: Note the condition on the s,t coefficients in the relations
+    being 1 or -1 for this optimized function.  There is a more
+    general function in relation_matrix.py, which is much, much
+    slower.
+    
+    INPUT:
+    
+        - ``rels`` - set of pairs ((i,s), (j,t)). The pair represents
+          the relation s*x_i + t*x_j = 0, where the i, j must be
+          Python int's, and the s,t must all be 1 or -1.
+
+        - ``n`` - int, the x_i are x_0, ..., x_n-1.
+    
+    OUTPUT:
+    
+        -  ``mod`` - list such that mod[i] = (j,s), which means
+           that x_i is equivalent to s*x_j, where the x_j are a basis for
+           the quotient.
+
+    EXAMPLES::
+
+        sage: v = [((0,1), (1,-1)), ((1,1), (3,1)), ((2,1),(3,1)), ((4,1),(5,-1))]
+        sage: rels = set(v)
+        sage: n = 6
+        sage: from sage.modular.modsym.relation_matrix_pyx import sparse_2term_quotient_only_pm1
+        sage: sparse_2term_quotient_only_pm1(rels, n)
+        [(3, -1), (3, -1), (3, -1), (3, 1), (5, 1), (5, 1)]        
+    """
+    if not isinstance(rels, set):
+        raise TypeError, "rels must be a set"
+
+    n = int(n)
+    
+    tm = misc.verbose("Starting optimized integer sparse 2-term quotient...")
+
+    cdef int c0, c1, i, die
+    free = range(n)
+    coef = [1]*n
+    related_to_me = [[] for i in range(n)]
+    
+    for v0, v1 in rels:
+        c0 = coef[v0[0]] * v0[1]
+        c1 = coef[v1[0]] * v1[1]
+
+        # Mod out by the following relation:
+        #
+        #    c0*free[v0[0]] + c1*free[v1[0]] = 0.
+        #
+        die = -1
+        if c0 == 0 and c1 == 0:
+            pass
+        elif c0 == 0 and c1 != 0:  # free[v1[0]] --> 0
+            die = free[v1[0]]
+        elif c1 == 0 and c0 != 0:
+            die = free[v0[0]]
+        elif free[v0[0]] == free[v1[0]]:
+            if c0 + c1 != 0:
+                # all xi equal to free[v0[0]] must now equal to zero.
+                die = free[v0[0]]
+        else:  # x1 = -c1/c0 * x2.
+            x = free[v0[0]]
+            free[x] = free[v1[0]]
+            if c0 != 1 and c0 != -1:
+                raise ValueError, "coefficients must all be -1 or 1."
+            coef[x] = -c1/c0
+            for i in related_to_me[x]:
+                free[i] = free[x]
+                coef[i] *= coef[x]
+                related_to_me[free[v1[0]]].append(i)
+            related_to_me[free[v1[0]]].append(x)
+        if die != -1:
+            for i in related_to_me[die]:
+                free[i] = 0
+                coef[i] = 0
+            free[die] = 0
+            coef[die] = 0
+
+    # Special casing the rationals leads to a huge speedup,
+    # actually.  (All the code above is slower than just this line
+    # without this special case.)
+    mod = [(free[i], Rational(coef[i])) for i in range(len(free))]
+        
+    misc.verbose("finished",tm)
+    return mod
+