Python arith.gcd函数代码示例

    def iter_positive_forms_with_content(self) :
        if self.__disc is infinity :
            raise ValueError, "infinity is not a true filter index"
        
        
        if self.__reduced :        
            for a in xrange(1,isqrt(self.__disc // 3) + 1) :
                for b in xrange(a+1) :
                    g = gcd(a, b)
                    for c in xrange(a, (b**2 + (self.__disc - 1))//(4*a) + 1) :
                        yield (a,b,c), gcd(g,c)
        else :
            maxtrace = floor(5*self.__disc / 15 + sqrt(self.__disc)/2)
            for a in xrange(1, maxtrace + 1) :
                for c in xrange(1, maxtrace - a + 1) :
                    g = gcd(a,c)
                    
                    Bu = isqrt(4*a*c - 1)

                    di = 4*a*c - self.__disc
                    if di >= 0 :
                        Bl = isqrt(di) + 1 
                    else :
                        Bl = 0
                    
                    for b in xrange(-Bu, -Bl + 1) :
                        yield (a,b,c), gcd(g,b)
                    for b in xrange(Bl, Bu + 1) :
                        yield (a,b,c), gcd(g,b)
        #! if self.__reduced

        raise StopIteration

    def _find_cusps(self):
        r"""
        Return an ordered list of inequivalent cusps for self, i.e. a
        set of representatives for the orbits of self on
        `\mathbf{P}^1(\QQ)`.  These are returned in a reduced
        form; see self.reduce_cusp for the definition of reduced.

        ALGORITHM:
            Lemma 3.2 in Cremona's 1997 book shows that for the action
            of Gamma1(N) on "signed projective space"
            `\Q^2 / (\Q_{\geq 0}^+)`, we have `u_1/v_1 \sim u_2 / v_2`
            if and only if `v_1 = v_2 \bmod N` and `u_1 = u_2 \bmod
            gcd(v_1, N)`. It follows that every orbit has a
            representative `u/v` with `v \le N` and `0 \le u \le
            gcd(v, N)`.  We iterate through all pairs `(u,v)`
            satisfying this.

            Having found a set containing at least one of every
            equivalence class modulo Gamma1(N), we can be sure of
            picking up every class modulo GammaH(N) since this
            contains Gamma1(N); and the reduce_cusp call does the
            checking to make sure we don't get any duplicates.

        EXAMPLES::

            sage: Gamma1(5)._find_cusps()
            [0, 2/5, 1/2, Infinity]
            sage: Gamma1(35)._find_cusps()
            [0, 2/35, 1/17, 1/16, 1/15, 1/14, 1/13, 1/12, 3/35, 1/11, 1/10, 1/9, 4/35, 1/8, 2/15, 1/7, 1/6, 6/35, 1/5, 3/14, 8/35, 1/4, 9/35, 4/15, 2/7, 3/10, 11/35, 1/3, 12/35, 5/14, 13/35, 2/5, 3/7, 16/35, 17/35, 1/2, 8/15, 4/7, 3/5, 9/14, 7/10, 5/7, 11/14, 4/5, 6/7, 9/10, 13/14, Infinity]
            sage: Gamma1(24)._find_cusps() == Gamma1(24).cusps(algorithm='modsym')
            True
            sage: GammaH(24, [13,17])._find_cusps() == GammaH(24,[13,17]).cusps(algorithm='modsym')
            True
        """

        s = []
        hashes = []
        N = self.level()

        for d in xrange(1, 1+N):
            w = N.gcd(d)
            M = int(w) if w > 1 else 2
            for a in xrange(1,M):
                if gcd(a, w) != 1:
                    continue
                while gcd(a, d) != 1:
                    a += w
                c = self.reduce_cusp(Cusp(a,d))
                h = hash(c)
                if not h in hashes:
                    hashes.append(h)
                    s.append(c)
        return sorted(s)

def num_cusps_of_width(N, d):
    r"""
    Return the number of cusps on `X_0(N)` of width d.

    INPUT:


    -  ``N`` - (integer): the level

    -  ``d`` - (integer): an integer dividing N, the cusp
       width


    EXAMPLES::

        sage: [num_cusps_of_width(18,d) for d in divisors(18)]
        [1, 1, 2, 2, 1, 1]
    """
    try:
        N = ZZ(N)
        d = ZZ(d)
        assert N>0
        assert d>0
        assert ((N % d) == 0)
    except TypeError:
        raise TypeError, "N and d must be integers"
    except AssertionError:
        raise AssertionError, "N and d must be positive integers with d|N"

    return euler_phi(gcd(d, N//d))

def _normalize_H(H, level):
    """
    Normalize representatives for a given subgroup H of the units
    modulo level.

    NOTE: This function does *not* make any attempt to find a minimal
    set of generators for H. It simply normalizes the inputs for use
    in hashing.

    EXAMPLES::

        sage: sage.modular.arithgroup.congroup_gammaH._normalize_H([23], 10)
        [3]
        sage: sage.modular.arithgroup.congroup_gammaH._normalize_H([1,5], 7)
        [5]
        sage: sage.modular.arithgroup.congroup_gammaH._normalize_H([4,18], 14)
        Traceback (most recent call last):
        ...
        ArithmeticError: The generators [4, 18] must be units modulo 14
        sage: sage.modular.arithgroup.congroup_gammaH._normalize_H([3,17], 14)
        [3]
        sage: sage.modular.arithgroup.congroup_gammaH._normalize_H([-1,7,9], 10)
        [7, 9]
    """
    H = [ZZ(h) for h in H]
    for h in H:
        if gcd(h, level) > 1:
            raise ArithmeticError, 'The generators %s must be units modulo %s'%(H, level)
    H = list(set([h%level for h in H]))
    H.sort()
    if 1 in H:
        H.remove(1)
    return H

def _list_subgroup(N, gens):
    r"""
    Given an integer ``N`` and a list of integers ``gens``, return a list of
    the elements of the subgroup of `(\ZZ / N\ZZ)^\times` generated by the
    elements of ``gens``.

    EXAMPLE::

        sage: sage.modular.arithgroup.congroup_gammaH._list_subgroup(11, [3])
        [1, 3, 4, 5, 9]
    """
    H = set([1])
    N = int(N)
    for g in gens:
        if gcd(g, N) != 1:
            raise ValueError, "gen (=%s) is not in (Z/%sZ)^*"%(g,N)
        gk = int(g) % N
        sbgrp = [gk]
        while not (gk in H):
            gk = (gk * g)%N
            sbgrp.append(gk)
        H = set([(x*h)%N for x in sbgrp for h in H])
    H = list(H)
    H.sort()
    return H

def test_Hecke_relations(a,b,C):
    r"""Testing Hecke relations for the Fourier coefficients in C



    INPUT:
    -''C'' -- dictionary of complex (Fourier coefficients)
    -''a'' -- integer
    -''b'' -- integer

    OUTPUT:
    -''diff'' -- real : |C(a)C(b)-C(ab)| if (a,b)=1

    EXAMPLE::

    
    sage: S=MaassWaveForms(Gamma0(1))
    sage: R=mpmath.mpf(9.53369526135355755434423523592877032382125639510725198237579046413534)
    sage: Y=mpmath.mpf(0.85)
    sage: C=coefficients_for_Maass_waveforms(S,R,Y,10,20,12)
    sage: d=test_Hecke_relations(C,2,3); mppr(d)
    '9.29e-8'
    sage: C=coefficients_for_Maass_waveforms(S,R,Y,30,50,20)
    sage: d=test_Hecke_relations(C,2,3); mppr(d)
    '3.83e-43'
    
    
    """
    c=gcd(Integer(a),Integer(b))
    lhs=C[0][a]*C[0][b]
    rhs=mpmath.mpf(0)
    for d in divisors(c):
        rhs=rhs+C[0][Integer(a*b/d/d)]
    return abs(rhs-lhs)

    def order_at_cusp(self, cusp):
        r"""
        Return the order of vanishing of self at the given cusp.

        INPUT:


        -  ``cusp`` -  a CuspFamily object


        OUTPUT:


        - an integer


        EXAMPLES::

            sage: e = EtaProduct(2, {2:24, 1:-24})
            sage: e.order_at_cusp(CuspFamily(2, 1)) # cusp at infinity
            1
            sage: e.order_at_cusp(CuspFamily(2, 2)) # cusp 0
            -1
        """
        if not isinstance(cusp, CuspFamily):
            raise TypeError("Argument (=%s) should be a CuspFamily" % cusp)
        if cusp.level() != self.level():
            raise ValueError("Cusp not on right curve!")
        return 1/ZZ(24)/gcd(cusp.width(), self.level()//cusp.width()) * sum( [ell*self.r(ell)/cusp.width() * (gcd(cusp.width(), self.level()//ell))**2  for ell in self._keys] )

def num_cusps_of_width(N, d):
    r"""
    Return the number of cusps on `X_0(N)` of width d.

    INPUT:


    -  ``N`` - (integer): the level

    -  ``d`` - (integer): an integer dividing N, the cusp
       width


    EXAMPLES::

        sage: [num_cusps_of_width(18,d) for d in divisors(18)]
        [1, 1, 2, 2, 1, 1]
        sage: num_cusps_of_width(4,8)
        Traceback (most recent call last):
        ...
        ValueError: N and d must be positive integers with d|N
    """
    N = ZZ(N)
    d = ZZ(d)
    if N <= 0 or d <= 0 or (N % d) != 0:
        raise ValueError("N and d must be positive integers with d|N")

    return euler_phi(gcd(d, N//d))

    def _test__jacobi_predicted_taylor_coefficients(fs, q_precision) :
        r"""
        Given a list of power series, which are the corrected Taylor coefficients
        of a Jacobi form, return the renormalized uncorrected ones, assuming that
        all but one `f` vanish.
        
        INPUT:
        
        - ``fs`` -- A list of power series.
        
        - ``q_precision`` -- An integer.
        
        OUPUT:
        
        - A list of power series.
        
        TESTS:
        
        See jacobi_form_by_taylor_expansion.
        """
        from sage.rings.arith import gcd
        
        R = PowerSeriesRing(ZZ, 'q'); q = R.gen(0)
        
        diff = lambda f: f.derivative().shift(1)
        normalize = lambda f: f / gcd(f.list()) if f != 0 else f
        diffnorm = lambda f,l: normalize(reduce(lambda a, g: g(a), l*[diff], f))

        taylor_coefficients = list()
        allf = R(0)
        for f in fs :
            allf = f(q_precision) + diffnorm(allf, 1)            
            taylor_coefficients.append(allf)

        return taylor_coefficients

    def random_element(self, bound=100, *args, **kwds):
        r"""
        Return a random element of `{\rm SL}_2(\ZZ)` with entries whose
        absolute value is strictly less than bound (default 100).
        Additional arguments and keywords are passed to the random_element
        method of ZZ.

        (Algorithm: Generate a random pair of integers at most bound. If they
        are not coprime, throw them away and start again. If they are, find an
        element of `{\rm SL}_2(\ZZ)` whose bottom row is that, and
        left-multiply it by `\begin{pmatrix} 1 & w \\ 0 & 1\end{pmatrix}` for
        an integer `w` randomly chosen from a small enough range that the
        answer still has entries at most bound.)

        It is, unfortunately, not true that all elements of SL2Z with entries <
        bound appear with equal probability; those with larger bottom rows are
        favoured, because there are fewer valid possibilities for w.

        EXAMPLES::

            sage: SL2Z.random_element()
            [60 13]
            [83 18]
            sage: SL2Z.random_element(5)
            [-1  3]
            [ 1 -4]

        Passes extra positional or keyword arguments through::

            sage: SL2Z.random_element(5, distribution='1/n')
            [ 1 -4]
            [ 0  1]
        """
        if bound <= 1: raise ValueError("bound must be greater than 1")
        c = ZZ.random_element(1-bound, bound, *args, **kwds)
        d = ZZ.random_element(1-bound, bound, *args, **kwds)
        if gcd(c,d) != 1: # try again
            return self.random_element(bound, *args, **kwds)
        else:
            a,b,c,d = lift_to_sl2z(c,d,0)
            whi = bound
            wlo = bound
            if c > 0:
                whi = min(whi, ((bound - a)/ZZ(c)).ceil())
                wlo = min(wlo, ((bound + a)/ZZ(c)).ceil())
            elif c < 0:
                whi = min(whi, ((bound + a)/ZZ(-c)).ceil())
                wlo = min(wlo, ((bound - a)/ZZ(-c)).ceil())

            if d > 0:
                whi = min(whi, ((bound - b)/ZZ(d)).ceil())
                wlo = min(wlo, ((bound + b)/ZZ(d)).ceil())
            elif d < 0:
                whi = min(whi, ((bound + b)/ZZ(-d)).ceil())
                wlo = min(wlo, ((bound - b)/ZZ(-d)).ceil())

            w = ZZ.random_element(1-wlo, whi, *args, **kwds)
            a += c*w
            b += d*w
            return self([a,b,c,d])

 def eval(self, expansion, weight = None) :
     if weight is None :
         try :
             weight = expansion.weight()
         except AttributeError :
             raise ValueError, "weight must be defined for the Hecke action"
     
     precision = expansion.precision()
     if precision.is_infinite() :
         precision = expansion._bounding_precision()
     else :
         precision = precision._hecke_operator(self.__l)
     characters = expansion.non_zero_components()
     expansion_level = precision.level()
     
     if gcd(expansion_level, self.__l) != 1 :
         raise ValueError, "Level of expansion and of Hecke operator must be coprime."
     
     hecke_expansion = dict()
     for ch in characters :
         hecke_expansion[ch] = dict( (k, self.hecke_coeff(expansion, ch, k, weight, expansion_level))
                                     for k in precision )
     
     result = expansion.parent()._element_constructor_(hecke_expansion)
     result._set_precision(expansion.precision()._hecke_operator(self.__l))
     
     return result

    def get_representatives( self, t, N) :
        r"""
        A helper function used in hecke_coeff that computes the right
        coset representatives of $\Gamma^0(t) \cap \Gamma_0(N) \ Gamma_0(N)$ where
        $\Gamma^0(t)$ is the subgroup of $SL(2,Z)$ where the upper right hand
        corner is divisible by $t$.

        NOTE
            We use the bijection $\Gamma^0(t)\SL(2,Z) \rightarrow P^1(\Z/t\Z)$
            given by $A \mapsto [1:0]A$.
        """
        if t == 1 : return [(1,0,0,1)]
        
        rep_list = []
        
        for (x,y) in P1List(t):
            ## we know that (N, x, y) = 1
            N1 = gcd(N, x)
            if N1 != 1 :
                x = x + t * N // N1

            ## we calculate a pair c,d satisfying a minimality condition
            ## to make later multiplications cheaper
            (_, d, c) = Integer(x)._xgcd(Integer(N * y), minimal=True)
            
            #print (x, y, -N * c, d)
            rep_list.append((x, y, -N * c, d))
                
        return rep_list

def circle_drops(A,B):
    # Drops going around the unit circle for those A and B.
    # See http://user.math.uzh.ch/dehaye/thesis_students/Nicolas_Wider-Integrality_of_factorial_ratios.pdf
    # for longer description (not original, better references exist)    
    marks = lcm(lcm(A),lcm(B))    
    tmp = [0 for i in range(marks)]
    for a in A:
#        print tmp
        for i in range(a):
            if gcd(i, a) == 1:
                tmp[i*marks/a] -= 1
    for b in B:
#        print tmp
        for i in range(b):
            if gcd(i, b) == 1:
                tmp[i*marks/b] += 1
#    print tmp
    return [sum(tmp[:j]) for j in range(marks)]

def blift(LF,Li,p,S=None):
    r"""
    Search for a solution to the given list of inequalities. If found, lift the solution to
    an appropriate valuation. See Lemma 3.3.6 in [Molnar, M.Sc. thesis]

    INPUT:

    - ``LF`` -- a list of integer polynomials in one variable (the normalized coefficients)

    - ``Li`` -- an integer, the bound on coefficients

    - ``p`` -- a prime

    OUTPUT:

    - Boolean -- whether or not the lift is successful

    - integer -- the lift

    EXAMPLES::

        sage: R.<b> = PolynomialRing(QQ)
        sage: from sage.schemes.projective.endPN_minimal_model import blift
        sage: blift([8*b^3 + 12*b^2 + 6*b + 1, 48*b^2 + 483*b + 117, 72*b + 1341, -24*b^2 + 411*b + 99, -144*b + 1233, -216*b], 2, 3)
        (True, 4)
    """

    P = LF[0].parent()
    #Determine which inequalities are trivial, and scale the rest, so that we only lift
    #as many times as needed.
    keepScaledIneqs = [scale(P(coeff),Li,p) for coeff in LF if coeff != 0]
    keptVals = [i[2] for i in keepScaledIneqs if i[0]]
    if keptVals != []:
        #Determine the valuation to lift until.
        liftval = max(keptVals)
    else:
        #All inequalities are satisfied.
        return True,1
    if S is None:
        S = PolynomialRing(Zmod(p),'b')
    keptScaledIneqs = [S(i[1]) for i in keepScaledIneqs if i[0]]
    #We need a solution for each polynomial on the left hand side of the inequalities,
    #so we need only find a solution for their gcd.
    g = gcd(keptScaledIneqs)
    rts = g.roots(multiplicities=False)
    for r in rts:
        #Recursively try to lift each root
        r_initial=QQ(r)
        newInput = P([r_initial,p])
        LG = [F(newInput) for F in LF]
        lift,lifted = blift(LG,Li,p,S=S)
        if lift:
            #Lift successful.
            return True,r_initial+ p*lifted
    #Lift non successful.
    return False,0

 def __contains__(self, f) :
     if self.__disc is infinity :
         return True
     
     (a, b, c) = f
     disc = 4*a*c - b**2
     if disc == 0 :
         return gcd([a, b, c]) < self._indefinite_content_bound()
     else :
         return disc < self.__disc