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

File trac-8614-optimize-modular-symbol-relations-rebase.patch, 8.7 KB (added by mraum, 10 years ago)
  • module_list.py

    # HG changeset patch
    # User Martin Raum <Martin.Raum@matha.rwth-aachen.de>
    # Date 1300754596 -3600
    # Node ID 14d5de3ffc211021f81b0a4416b2a3a426718a0c
    # Parent  cefe7ffc4322f0dc727cf001d09a6e7dd6659d83
    #Trac 8614: Optimize creation of modular symbols spaces by speeding up quotienting out by 2-term relations
    
    diff -r cefe7ffc4322 -r 14d5de3ffc21 module_list.py
    a b  
    10281028              extra_compile_args=["-std=c99",  "-D_XPG6"],
    10291029              depends = flint_depends),
    10301030
     1031    Extension('sage.modular.modsym.relation_matrix_pyx',
     1032              sources = ['sage/modular/modsym/relation_matrix_pyx.pyx']),
     1033
    10311034    Extension('sage.modular.modsym.heilbronn',
    10321035              sources = ['sage/modular/modsym/heilbronn.pyx'],
    10331036              libraries = ["csage", "flint", "gmp", "gmpxx", "m", "stdc++"],
  • sage/modular/modsym/manin_symbols.py

    diff -r cefe7ffc4322 -r 14d5de3ffc21 sage/modular/modsym/manin_symbols.py
    a b  
    182182
    183183        """
    184184        raise NotImplementedError, "Only implemented in derived classes"
     185
     186    def _apply_S_only_0pm1(self):
     187        """
     188        Return True if the coefficient when applying the S relation is
     189        always 0, 1, or -1.  This is useful for optimizing code in
     190        relation_matrix.py.
     191
     192        EXAMPLES::
     193
     194            sage: eps = DirichletGroup(4).gen(0)
     195            sage: from sage.modular.modsym.manin_symbols import ManinSymbolList_character
     196            sage: ManinSymbolList_character(eps,2)._apply_S_only_0pm1()
     197            True
     198            sage: eps = DirichletGroup(7).gen(0)
     199            sage: from sage.modular.modsym.manin_symbols import ManinSymbolList_character
     200            sage: ManinSymbolList_character(eps,2)._apply_S_only_0pm1()
     201            False
     202        """
     203        return False # derived classes could overload and put True
    185204       
    186205    def apply_S(self, j):
    187206        """
     
    501520        else:
    502521            return k, -1
    503522       
     523    def _apply_S_only_0pm1(self):
     524        """
     525        Return True if the coefficient when applying the S relation is
     526        always 0, 1, or -1.  This is useful for optimizing code in
     527        relation_matrix.py.
     528
     529        EXAMPLES::
     530
     531            sage: from sage.modular.modsym.manin_symbols import ManinSymbolList_gamma0
     532            sage: ManinSymbolList_gamma0(5,8)._apply_S_only_0pm1()
     533            True
     534        """
     535        return True
     536
    504537    def apply_I(self, j):
    505538        """
    506539        Apply the matrix `I=[-1,0,0,1]` to the `j`-th Manin symbol.
     
    10471080        else:
    10481081            return k, -s
    10491082       
     1083    def _apply_S_only_0pm1(self):
     1084        """
     1085        Return True if the coefficient when applying the S relation is
     1086        always 0, 1, or -1.  This is useful for optimizing code in
     1087        relation_matrix.py.
     1088
     1089        EXAMPLES::
     1090
     1091            sage: eps = DirichletGroup(4).gen(0)
     1092            sage: from sage.modular.modsym.manin_symbols import ManinSymbolList_character
     1093            sage: ManinSymbolList_character(eps,2)._apply_S_only_0pm1()
     1094            True
     1095            sage: ManinSymbolList_character(DirichletGroup(13).0,2)._apply_S_only_0pm1()
     1096            False
     1097        """
     1098        return self.__character.order() <= 2
     1099
    10501100    def apply_I(self, j):
    10511101        """
    10521102        Apply the matrix `I=[-1,0,0,1]` to the `j`-th Manin symbol.
  • sage/modular/modsym/p1list.pyx

    diff -r cefe7ffc4322 -r 14d5de3ffc21 sage/modular/modsym/p1list.pyx
    a b  
    678678            sage: L
    679679            The projective line over the integers modulo 120
    680680        """
     681        cdef int i
     682
    681683        self.__N = N
    682684        if N <= 46340:
    683685            self.__list = p1list_int(N)
     
    700702        if not self.t: raise MemoryError
    701703
    702704        # Initialize xgcd table
    703         cdef int i
    704705        cdef llong ll_s, ll_t, ll_N = N
    705706
    706707        if N <= 46340:
  • sage/modular/modsym/relation_matrix.py

    diff -r cefe7ffc4322 -r 14d5de3ffc21 sage/modular/modsym/relation_matrix.py
    a b  
    459459    if sign != 0:
    460460        # Let rels = rels union I relations.
    461461        rels.update(modI_relations(syms,sign))
    462     mod = sparse_2term_quotient(rels, len(syms), field)
     462
     463    if syms._apply_S_only_0pm1() and rings.is_RationalField(field):
     464        import relation_matrix_pyx
     465        mod = relation_matrix_pyx.sparse_2term_quotient_only_pm1(rels, len(syms))
     466    else:
     467        mod = sparse_2term_quotient(rels, len(syms), field)
     468       
    463469    R = T_relation_matrix_wtk_g0(syms, mod, field, sparse)
    464470    return R, mod
    465471   
  • new file sage/modular/modsym/relation_matrix_pyx.pyx

    diff -r cefe7ffc4322 -r 14d5de3ffc21 sage/modular/modsym/relation_matrix_pyx.pyx
    - +  
     1"""
     2Optimized Cython code for computing relation matrices in certain cases.
     3"""
     4
     5#############################################################################
     6#       Copyright (C) 2010 William Stein <wstein@gmail.com>
     7#  Distributed under the terms of the GNU General Public License (GPL) v2+.
     8#  The full text of the GPL is available at:
     9#                  http://www.gnu.org/licenses/
     10#############################################################################
     11
     12import sage.misc.misc as misc
     13from sage.rings.rational cimport Rational
     14
     15def sparse_2term_quotient_only_pm1(rels, n):
     16    r"""
     17    Performs Sparse Gauss elimination on a matrix all of whose columns
     18    have at most 2 nonzero entries with relations all 1 or -1.
     19
     20    This algorithm is more subtle than just "identify symbols in pairs",
     21    since complicated relations can cause generators to equal 0.
     22
     23    NOTE: Note the condition on the s,t coefficients in the relations
     24    being 1 or -1 for this optimized function.  There is a more
     25    general function in relation_matrix.py, which is much, much
     26    slower.
     27   
     28    INPUT:
     29   
     30        - ``rels`` - set of pairs ((i,s), (j,t)). The pair represents
     31          the relation s*x_i + t*x_j = 0, where the i, j must be
     32          Python int's, and the s,t must all be 1 or -1.
     33
     34        - ``n`` - int, the x_i are x_0, ..., x_n-1.
     35   
     36    OUTPUT:
     37   
     38        -  ``mod`` - list such that mod[i] = (j,s), which means
     39           that x_i is equivalent to s*x_j, where the x_j are a basis for
     40           the quotient.
     41
     42    EXAMPLES::
     43
     44        sage: v = [((0,1), (1,-1)), ((1,1), (3,1)), ((2,1),(3,1)), ((4,1),(5,-1))]
     45        sage: rels = set(v)
     46        sage: n = 6
     47        sage: from sage.modular.modsym.relation_matrix_pyx import sparse_2term_quotient_only_pm1
     48        sage: sparse_2term_quotient_only_pm1(rels, n)
     49        [(3, -1), (3, -1), (3, -1), (3, 1), (5, 1), (5, 1)]       
     50    """
     51    if not isinstance(rels, set):
     52        raise TypeError, "rels must be a set"
     53
     54    n = int(n)
     55   
     56    tm = misc.verbose("Starting optimized integer sparse 2-term quotient...")
     57
     58    cdef int c0, c1, i, die
     59    free = range(n)
     60    coef = [1]*n
     61    related_to_me = [[] for i in range(n)]
     62   
     63    for v0, v1 in rels:
     64        c0 = coef[v0[0]] * v0[1]
     65        c1 = coef[v1[0]] * v1[1]
     66
     67        # Mod out by the following relation:
     68        #
     69        #    c0*free[v0[0]] + c1*free[v1[0]] = 0.
     70        #
     71        die = -1
     72        if c0 == 0 and c1 == 0:
     73            pass
     74        elif c0 == 0 and c1 != 0:  # free[v1[0]] --> 0
     75            die = free[v1[0]]
     76        elif c1 == 0 and c0 != 0:
     77            die = free[v0[0]]
     78        elif free[v0[0]] == free[v1[0]]:
     79            if c0 + c1 != 0:
     80                # all xi equal to free[v0[0]] must now equal to zero.
     81                die = free[v0[0]]
     82        else:  # x1 = -c1/c0 * x2.
     83            x = free[v0[0]]
     84            free[x] = free[v1[0]]
     85            if c0 != 1 and c0 != -1:
     86                raise ValueError, "coefficients must all be -1 or 1."
     87            coef[x] = -c1/c0
     88            for i in related_to_me[x]:
     89                free[i] = free[x]
     90                coef[i] *= coef[x]
     91                related_to_me[free[v1[0]]].append(i)
     92            related_to_me[free[v1[0]]].append(x)
     93        if die != -1:
     94            for i in related_to_me[die]:
     95                free[i] = 0
     96                coef[i] = 0
     97            free[die] = 0
     98            coef[die] = 0
     99
     100    # Special casing the rationals leads to a huge speedup,
     101    # actually.  (All the code above is slower than just this line
     102    # without this special case.)
     103    mod = [(free[i], Rational(coef[i])) for i in range(len(free))]
     104       
     105    misc.verbose("finished",tm)
     106    return mod
     107