Ticket #8614: trac_8614.patch

File trac_8614.patch, 15.9 KB (added by was, 11 years ago)
  • module_list.py

    # HG changeset patch
    # User William Stein <wstein@gmail.com>
    # Date 1269661144 25200
    # Node ID 9e66b6e3bde8ff31dbc22188424080d6003cf14f
    # Parent  cf4cd2737516b7cc0462a7906cfb31744932c9c2
    trac 8614 -- Optimize creation of modular symbols spaces by speeding up quotienting out by 2-term relations
    
    diff -r cf4cd2737516 -r 9e66b6e3bde8 module_list.py
    a b  
    843843              extra_compile_args=["-std=c99",  "-D_XPG6"],
    844844              depends = [SAGE_ROOT + "/local/include/FLINT/flint.h"]),
    845845
     846    Extension('sage.modular.modsym.relation_matrix_pyx',
     847              sources = ['sage/modular/modsym/relation_matrix_pyx.pyx']),
     848
    846849    Extension('sage.modular.modsym.heilbronn',
    847850              sources = ['sage/modular/modsym/heilbronn.pyx'],
    848851              libraries = ["csage", "flint", "gmp", "gmpxx", "m", "stdc++"],
  • sage/modular/modform/ambient.py

    diff -r cf4cd2737516 -r 9e66b6e3bde8 sage/modular/modform/ambient.py
    a b  
    394394        self.__module = free_module.VectorSpace(self.base_ring(), d)
    395395        return self.__module
    396396
    397     def free_module(self): return self.module()
    398     # stupid thing: there are functions in classes ModularFormsSpace and
    399     # HeckeModule that both do much the same thing, and one has to override
    400     # both of them!
     397    # free_module -- stupid thing: there are functions in classes
     398    # ModularFormsSpace and HeckeModule that both do much the same
     399    # thing, and one has to override both of them!
     400    def free_module(self):
     401        """
     402        Return free module underlying this modular forms space.
     403
     404        EXAMPLES::
     405
     406            sage: ModularForms(37).free_module()
     407            Vector space of dimension 3 over Rational Field
     408        """
     409        return self.module()
    401410
    402411    def prec(self, new_prec=None):
    403412        """
  • sage/modular/modsym/ambient.py

    diff -r cf4cd2737516 -r 9e66b6e3bde8 sage/modular/modsym/ambient.py
    a b  
    9191import sage.structure.formal_sum as formal_sum
    9292import sage.categories.all as cat
    9393from sage.modular.cusps import Cusp
    94 from sage.rings.arith import binomial
    9594import sage.structure.all
    9695
    9796import boundary
  • sage/modular/modsym/apply.pyx

    diff -r cf4cd2737516 -r 9e66b6e3bde8 sage/modular/modsym/apply.pyx
    a b  
    2323
    2424cdef class Apply:
    2525    def __cinit__(self):
     26        """
     27        EXAMPLES:::
     28
     29            sage: import sage.modular.modsym.apply
     30            sage: sage.modular.modsym.apply.Apply()
     31            <sage.modular.modsym.apply.Apply object at ...>
     32        """
    2633        fmpz_poly_init(self.f)
    2734        fmpz_poly_init(self.g)
    2835        fmpz_poly_init(self.ff)
     
    5764cdef Apply A = Apply()       
    5865
    5966def apply_to_monomial(int i, int j, int a, int b, int c, int d):
    60     """
     67    r"""
    6168    Returns a list of the coefficients of
    6269    $$
    6370             (aX + bY)^i (cX + dY)^{j-i},
     
    7582        of $Y^j$, $Y^{j-1}*X$, \ldots, $X^j$, respectively.
    7683       
    7784    EXAMPLE:
    78     We compute that $(X+Y)^2(X-Y) = X^3 + X^2Y - XY^2 - Y^3$.
     85   
     86    We compute that $(X+Y)^2(X-Y) = X^3 + X^2Y - XY^2 - Y^3$::
     87   
    7988        sage: from sage.modular.modsym.manin_symbols import apply_to_monomial
    8089        sage: apply_to_monomial(2, 3, 1,1,1,-1)
    8190        [-1, -1, 1, 1]
  • sage/modular/modsym/manin_symbols.py

    diff -r cf4cd2737516 -r 9e66b6e3bde8 sage/modular/modsym/manin_symbols.py
    a b  
    181181
    182182        """
    183183        raise NotImplementedError, "Only implemented in derived classes"
     184
     185    def _apply_S_only_0pm1(self):
     186        """
     187        Return True if the coefficient when applying the S relation is
     188        always 0, 1, or -1.  This is useful for optimizing code in
     189        relation_matrix.py.
     190
     191        EXAMPLES::
     192
     193            sage: eps = DirichletGroup(4).gen(0)
     194            sage: from sage.modular.modsym.manin_symbols import ManinSymbolList_character
     195            sage: ManinSymbolList_character(eps,2)._apply_S_only_0pm1()
     196            True
     197            sage: eps = DirichletGroup(7).gen(0)
     198            sage: from sage.modular.modsym.manin_symbols import ManinSymbolList_character
     199            sage: ManinSymbolList_character(eps,2)._apply_S_only_0pm1()
     200            False
     201        """
     202        return False # derived classes could overload and put True
    184203       
    185204    def apply_S(self, j):
    186205        """
     
    500519        else:
    501520            return k, -1
    502521       
     522    def _apply_S_only_0pm1(self):
     523        """
     524        Return True if the coefficient when applying the S relation is
     525        always 0, 1, or -1.  This is useful for optimizing code in
     526        relation_matrix.py.
     527
     528        EXAMPLES::
     529
     530            sage: from sage.modular.modsym.manin_symbols import ManinSymbolList_gamma0
     531            sage: ManinSymbolList_gamma0(5,8)._apply_S_only_0pm1()
     532            True
     533        """
     534        return True
     535
    503536    def apply_I(self, j):
    504537        """
    505538        Apply the matrix `I=[-1,0,0,1]` to the `j`-th Manin symbol.
     
    575608        else:
    576609            s = -1
    577610        z = []
     611        a = rings.ZZ(k-2-i)
    578612        for j in range(k-2-i +1):
    579613            m = self.index((j, u, v))
    580             z.append((m,s*arith.binomial(k-2-i,j)))
     614            z.append((m, s * a.binomial(j)))
    581615            s *= -1
    582616        return z
    583617
     
    617651        else:
    618652            s = -1
    619653        z = []
     654        a = rings.ZZ(i)
    620655        for j in range(i+1):
    621656            m = self.index((k-2-i+j, u, v))
    622             z.append((m,s*arith.binomial(i,j)))
     657            z.append((m, s*a.binomial(j)))
    623658            s *= -1
    624659        return z
    625660       
     
    10231058        OUTPUT:
    10241059
    10251060        ``(k, s)`` where `k` is the index of the symbol obtained by acting
    1026         on the `j`'th symbol with `S`, and `s` is the parity of of the
     1061        on the `j`'th symbol with `S`, and `s` is the parity of the
    10271062        `j`'th symbol.
    10281063       
    10291064        EXAMPLE::
     
    10441079        else:
    10451080            return k, -s
    10461081       
     1082    def _apply_S_only_0pm1(self):
     1083        """
     1084        Return True if the coefficient when applying the S relation is
     1085        always 0, 1, or -1.  This is useful for optimizing code in
     1086        relation_matrix.py.
     1087
     1088        EXAMPLES::
     1089
     1090            sage: eps = DirichletGroup(4).gen(0)
     1091            sage: from sage.modular.modsym.manin_symbols import ManinSymbolList_character
     1092            sage: ManinSymbolList_character(eps,2)._apply_S_only_0pm1()
     1093            True
     1094            sage: ManinSymbolList_character(DirichletGroup(13).0,2)._apply_S_only_0pm1()
     1095            False
     1096        """
     1097        return self.__character.order() <= 2
     1098
    10471099    def apply_I(self, j):
    10481100        """
    10491101        Apply the matrix `I=[-1,0,0,1]` to the `j`-th Manin symbol.
     
    11111163        else:
    11121164            s = -r
    11131165        z = []
     1166        a = rings.ZZ(k-2-i)
    11141167        for j in range(k-2-i +1):
    11151168            m, r = self.index((j, u, v))
    1116             z.append((m,s*r*arith.binomial(k-2-i,j)))
     1169            z.append((m,s*r*a.binomial(j)))
    11171170            s *= -1
    11181171        return z
    11191172
     
    11521205        else:
    11531206            s = -r
    11541207        z = []
     1208        a = rings.ZZ(i)
    11551209        for j in range(i+1):
    11561210            m, r = self.index((k-2-i+j, u, v))
    1157             z.append((m,s*r*arith.binomial(i,j)))
     1211            z.append((m,s*r*a.binomial(j)))
    11581212            s *= -1
    11591213        return z
    11601214                                 
  • sage/modular/modsym/p1list.pyx

    diff -r cf4cd2737516 -r 9e66b6e3bde8 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 cf4cd2737516 -r 9e66b6e3bde8 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   
     
    510516        sage: sparse_2term_quotient(rels, n, QQ)
    511517        [(3, -1/3), (3, -1), (3, -1), (3, 1), (5, 1), (5, 1)]
    512518    """
    513    
    514519    if not isinstance(rels, set):
    515520        raise TypeError, "rels must be a set"
    516521    n = int(n)
    517     #if not isinstance(n, int):
    518     #    raise TypeError, "n must be an int"
    519522    if not isinstance(F, rings.Ring):
    520523        raise TypeError, "F must be a ring."
    521524   
    522     tm = misc.verbose()
     525    tm = misc.verbose("Starting sparse 2-term quotient...")
    523526    free = range(n)
    524527    ONE = F(1)
    525528    ZERO = F(0)
     
    528531    for v0, v1 in rels:
    529532        c0 = coef[v0[0]] * F(v0[1])
    530533        c1 = coef[v1[0]] * F(v1[1])
    531         # Mod out by the relation
    532         #    c1*x_free[t[0]] + c2*x_free[t[1]] = 0.
     534
     535        # Mod out by the following relation:
     536        #
     537        #    c0*free[v0[0]] + c1*free[v1[0]] = 0.
     538        #
    533539        die = None
    534540        if c0 == ZERO and c1 == ZERO:
    535541            pass
    536         elif c0 == ZERO and c1 != ZERO:  # free[t[1]] --> 0
     542        elif c0 == ZERO and c1 != ZERO:  # free[v1[0]] --> 0
    537543            die = free[v1[0]]
    538544        elif c1 == ZERO and c0 != ZERO:
    539545            die = free[v0[0]]
    540546        elif free[v0[0]] == free[v1[0]]:
    541             if c0+c1 != 0:
    542                 # all xi equal to free[t[0]] must now equal to zero.
     547            if c0 + c1 != 0:
     548                # all xi equal to free[v0[0]] must now equal to zero.
    543549                die = free[v0[0]]
    544550        else:  # x1 = -c1/c0 * x2.
    545551            x = free[v0[0]]
     
    550556                coef[i] *= coef[x]
    551557                related_to_me[free[v1[0]]].append(i)
    552558            related_to_me[free[v1[0]]].append(x)
    553         if die != None:
     559        if die is not None:
    554560            for i in related_to_me[die]:
    555561                free[i] = 0
    556562                coef[i] = ZERO
  • new file sage/modular/modsym/relation_matrix_pyx.pyx

    diff -r cf4cd2737516 -r 9e66b6e3bde8 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
  • sage/rings/integer.pyx

    diff -r cf4cd2737516 -r 9e66b6e3bde8 sage/rings/integer.pyx
    a b  
    50875087        """
    50885088        return self
    50895089
     5090    def binomial(self, unsigned long int y, algorithm='mpir'):
     5091        """
     5092        Return the binomial coefficient "self choose y".
     5093
     5094        INPUT:
     5095
     5096            - y -- an unsigned long integer
     5097
     5098            - algorithm -- 'mpir' (default) or 'pari'; 'mpir' is
     5099              faster for small y, and 'pari' tends to be faster for
     5100              large y
     5101
     5102        OUTPUT:
     5103
     5104            - integer
     5105
     5106        EXAMPLES::
     5107
     5108            sage: 10.binomial(2)
     5109            45
     5110            sage: 10.binomial(2, algorithm='pari')
     5111            45
     5112        """
     5113        cdef Integer x
     5114        if algorithm == 'mpir':
     5115            x = PY_NEW(Integer)
     5116            mpz_bin_ui(x.value, self.value, y)
     5117            return x
     5118        elif algorithm == 'pari':
     5119            return the_integer_ring(self._pari_().binomial(y))
     5120           
     5121       
     5122
    50905123ONE = Integer(1)
    50915124Py_INCREF(ONE)
    50925125