def compute_G(p, F):
    r"""
    Given a power series `F \in R[[q]]^\times`, for some ring `R`, and an
    integer `p`, compute the quotient

    .. math::

        \frac{F(q)}{F(q^p)}.

    Used by :func:`level1_UpGj` and by :func:`higher_level_UpGj`, with `F` equal
    to the Eisenstein series `E_{p-1}`.

    INPUT:

    - ``p`` -- integer
    - ``F`` -- power series (with invertible constant term)

    OUTPUT:

    the power series `F(q) / F(q^p)`, to the same precision as `F`

    EXAMPLE::

        sage: E = sage.modular.overconvergent.hecke_series.eisenstein_series_qexp(2, 12, Zmod(9),normalization="constant")
        sage: sage.modular.overconvergent.hecke_series.compute_G(3, E)
        1 + 3*q + 3*q^4 + 6*q^7 + O(q^12)
    """
    Fp = (F.truncate_powerseries(ceil(F.prec() / ZZ(p)))).V(p)
    return F / Fp

def sector(ray1, ray2, **extra_options):
    r"""
    Plot a sector between ``ray1`` and ``ray2`` centered at the origin.

    .. NOTE::

        This function was intended for plotting strictly convex cones, so it
        plots the smaller sector between ``ray1`` and ``ray2`` and, therefore,
        they cannot be opposite. If you do want to use this function for bigger
        regions, split them into several parts.

    .. NOTE::

        As of version 4.6 Sage does not have a graphic primitive for sectors in
        3-dimensional space, so this function will actually approximate them
        using polygons (the number of vertices used depends on the angle
        between rays).

    INPUT:

    - ``ray1``, ``ray2`` -- rays in 2- or 3-dimensional space of the same
      length;

    - ``extra_options`` -- a dictionary of options that should be passed to
      lower level plotting functions.

    OUTPUT:

    - a plot.

    EXAMPLES::

        sage: from sage.geometry.toric_plotter import sector
        sage: sector((1,0), (0,1))
        Graphics object consisting of 1 graphics primitive
        sage: sector((3,2,1), (1,2,3))
        Graphics3d Object
    """
    ray1 = vector(RDF, ray1)
    ray2 = vector(RDF, ray2)
    r = ray1.norm()
    if len(ray1) == 2:
        # Plot an honest sector
        phi1 = arctan2(ray1[1], ray1[0])
        phi2 = arctan2(ray2[1], ray2[0])
        if phi1 > phi2:
            phi1, phi2 = phi2, phi1
        if phi2 - phi1 > pi:
            phi1, phi2 = phi2, phi1 + 2 * pi
        return disk((0,0), r, (phi1, phi2), **extra_options)
    else:
        # Plot a polygon, 30 vertices per radian.
        vertices_per_radian = 30
        n = ceil(arccos(ray1 * ray2 / r**2) * vertices_per_radian)
        dr = (ray2 - ray1) / n
        points = (ray1 + i * dr for i in range(n + 1))
        points = [r / pt.norm() * pt for pt in points]
        points.append(vector(RDF, 3))
        return polygon(points, **extra_options)

def native_two_isogeny_descent_work(E, two_tor_rk):
    """
    Prepares the output from two-descent by two-isogeny.

    INPUT:

        - ``E`` - an elliptic curve

        - ``two_tor_rk`` - its two-torsion rank

    OUTPUT:

        - a lower bound on the rank

        - an upper bound on the rank

        - a lower bound on the rank of Sha[2]

        - an upper bound on the rank of Sha[2]

        - a list of the generators found (currently None, since we don't store them)

    EXAMPLES::

        sage: from sage.schemes.elliptic_curves.BSD import native_two_isogeny_descent_work
        sage: E = EllipticCurve('14a')
        sage: native_two_isogeny_descent_work(E, E.two_torsion_rank())
        (0, 0, 0, 0, None)
        sage: E = EllipticCurve('65a')
        sage: native_two_isogeny_descent_work(E, E.two_torsion_rank())
        (1, 1, 0, 0, None)

    """
    from sage.schemes.elliptic_curves.descent_two_isogeny import two_descent_by_two_isogeny
    n1, n2, n1p, n2p = two_descent_by_two_isogeny(E)
    # bring n1 and n1p up to the nearest power of two
    two = ZZ(2) # otherwise "log" is symbolic >.<
    e1  = ceil(ZZ(n1).log(two))
    e1p = ceil(ZZ(n1p).log(two))
    e2  = ZZ(n2).log(two)
    e2p = ZZ(n2p).log(two)
    rank_lower_bd = e1 + e1p - 2
    rank_upper_bd = e2 + e2p - 2
    sha_upper_bd = e2 + e2p - e1 - e1p
    gens = None # right now, we are not keeping track of them
    return rank_lower_bd, rank_upper_bd, 0, sha_upper_bd, gens

    def set_params(lam, k):
        n = pow(2, ceil(log(lam**2 * k)/log(2))) # dim of poly ring, closest power of 2 to k(lam^2)
        q = next_prime(ZZ(2)**(8*k*lam) * n**k, proof=False) # prime modulus

        sigma = int(sqrt(lam * n))
        sigma_prime = lam * int(n**(1.5))

        return (n, q, sigma, sigma_prime, k)

    def initialize_fundom(self):
        F=self._F
        Pts=[self.embed_coords(a,b) for a in [_sage_const_0 ,_sage_const_1 ] for b in [_sage_const_0 ,_sage_const_1 ]]
        xmin=min([P.x for P in Pts])
        xmax=max([P.x for P in Pts])
        ymin=min([P.y for P in Pts])
        ymax=max([P.y for P in Pts])

        self._aplusbomega=dict([])
        rangea=self.rangea_gen(xmin-_sage_const_1 ,ymin-_sage_const_1 ,xmax+_sage_const_1 ,ymax+_sage_const_1 )
        for b in self._ranget:
            bw=b*self._w
            for a in rangea(b):
                self._aplusbomega[(a,b)]=a+bw

        dx=RealField()(_sage_const_0p49 )
        dy=dx
        Nx=ceil((xmax-xmin)/dx)
        Ny=ceil((ymax-ymin)/dy)
        self._holes=deque([])
        for ii in range(Nx):
            for jj in range(Ny):
                self._holes.append(Hole(self,[xmin+ii*dx,ymin+jj*dy,xmin+(ii+_sage_const_1 )*dx,ymin+(jj+_sage_const_1 )*dy]))

    def evaluate_hole(self,h,verifying=False,maxdepth=-_sage_const_1 ):
        if not h.overlaps_fundamental_domain():
            return _sage_const_0 
        regs=[[],[],[],[]]
        data=[]
        used=self._used_regions
        corners=h.corners()
        atws=[]
        w=self._w
        if(verifying):
            at_least_one=[False,False,False,False]

        for reg in used:
            passed=True
            for ii,P in enumerate(corners):
                if not reg[_sage_const_2 ].contains_point(P):
                    passed=False
                    if(not verifying):
                        break
                elif verifying:
                    at_least_one[ii]=True
            if(passed):
                return _sage_const_0 

        if(verifying):
            if(all(at_least_one)):
                return -_sage_const_1 
            else:
                if(h.depth()>maxdepth):
                    return -_sage_const_2 
                else:
                    return -_sage_const_1 

        rangea=self.rangea_gen(h.xmin,h.ymin,h.xmax,h.ymax)
        at_least_one=[False,False,False,False]

        l=_sage_const_100 
        B=ceil(_sage_const_4 /(l*(h.size**_sage_const_2 )))
        B=max([self._Nboundmin,B])
        #print 'B=',B,'Size ~ 10^',RealField(prec=10)((log(h.size)/log(10)))
        regions=self._master_regs.get_regions(min([B,self._Nbound]))
        nreg=self.test_regions(corners,at_least_one,embed_coords=self.embed_coords,regions=regions,Region=Region,aplusbomega=self._aplusbomega, ranget=self._ranget,rangea=rangea)
        if(nreg!=None):
            # used.append(nreg)
            used.extend(nreg)
            return _sage_const_0 
        return -_sage_const_1 

    def plot_lattice(self):
        r"""
        Plot the lattice (i.e. its points in the cut-off bounds of ``self``).

        OUTPUT:

        - a plot.

        EXAMPLES::

            sage: from sage.geometry.toric_plotter import ToricPlotter
            sage: tp = ToricPlotter(dict(), 2)
            sage: print tp.plot_lattice()
            Graphics object consisting of 1 graphics primitive
        """
        if not self.show_lattice:
            # Plot the origin anyway, otherwise rays/generators may look ugly.
            return self.plot_points([self.origin])
        d = self.dimension
        extra_options = self.extra_options
        if d == 1:
            points = ((x, 0) for x in range(ceil(self.xmin), floor(self.xmax) + 1))
        elif d == 2:
            points = (
                (x, y)
                for x in range(ceil(self.xmin), floor(self.xmax) + 1)
                for y in range(ceil(self.ymin), floor(self.ymax) + 1)
            )
        elif d == 3:
            points = (
                (x, y, z)
                for x in range(ceil(self.xmin), floor(self.xmax) + 1)
                for y in range(ceil(self.ymin), floor(self.ymax) + 1)
                for z in range(ceil(self.zmin), floor(self.zmax) + 1)
            )
        if self.mode == "round":
            r = 1.01 * self.radius  # To make sure integer values work OK.
            points = (pt for pt in points if vector(pt).norm() <= r)
        f = self.lattice_filter
        if f is not None:
            points = (pt for pt in points if f(pt))
        return self.plot_points(tuple(points))

        def coproduct_iterated(self, n=1):
            r"""
            Apply ``n`` coproducts to ``self``.

            .. TODO::

                Remove dependency on ``modules_with_basis`` methods.

            EXAMPLES::

                sage: Psi = NonCommutativeSymmetricFunctions(QQ).Psi()
                sage: Psi[2,2].coproduct_iterated(0)
                Psi[2, 2]
                sage: Psi[2,2].coproduct_iterated(2)
                Psi[] # Psi[] # Psi[2, 2] + 2*Psi[] # Psi[2] # Psi[2]
                 + Psi[] # Psi[2, 2] # Psi[] + 2*Psi[2] # Psi[] # Psi[2]
                 + 2*Psi[2] # Psi[2] # Psi[] + Psi[2, 2] # Psi[] # Psi[]

            TESTS::

                sage: p = SymmetricFunctions(QQ).p()
                sage: p[5,2,2].coproduct_iterated()
                p[] # p[5, 2, 2] + 2*p[2] # p[5, 2] + p[2, 2] # p[5]
                 + p[5] # p[2, 2] + 2*p[5, 2] # p[2] + p[5, 2, 2] # p[]
                sage: p([]).coproduct_iterated(3)
                p[] # p[] # p[] # p[]

            ::

                sage: Psi = NonCommutativeSymmetricFunctions(QQ).Psi()
                sage: Psi[2,2].coproduct_iterated(0)
                Psi[2, 2]
                sage: Psi[2,2].coproduct_iterated(3)
                Psi[] # Psi[] # Psi[] # Psi[2, 2] + 2*Psi[] # Psi[] # Psi[2] # Psi[2]
                 + Psi[] # Psi[] # Psi[2, 2] # Psi[] + 2*Psi[] # Psi[2] # Psi[] # Psi[2]
                 + 2*Psi[] # Psi[2] # Psi[2] # Psi[] + Psi[] # Psi[2, 2] # Psi[] # Psi[]
                 + 2*Psi[2] # Psi[] # Psi[] # Psi[2] + 2*Psi[2] # Psi[] # Psi[2] # Psi[]
                 + 2*Psi[2] # Psi[2] # Psi[] # Psi[] + Psi[2, 2] # Psi[] # Psi[] # Psi[]

            ::

                sage: m = SymmetricFunctionsNonCommutingVariables(QQ).m()
                sage: m[[1,3],[2]].coproduct_iterated(2)
                m{} # m{} # m{{1, 3}, {2}} + m{} # m{{1}} # m{{1, 2}}
                 + m{} # m{{1, 2}} # m{{1}} + m{} # m{{1, 3}, {2}} # m{}
                 + m{{1}} # m{} # m{{1, 2}} + m{{1}} # m{{1, 2}} # m{}
                 + m{{1, 2}} # m{} # m{{1}} + m{{1, 2}} # m{{1}} # m{}
                 + m{{1, 3}, {2}} # m{} # m{}
                sage: m[[]].coproduct_iterated(3), m[[1,3],[2]].coproduct_iterated(0)
                (m{} # m{} # m{} # m{}, m{{1, 3}, {2}})
            """
            if n < 0:
                raise ValueError("cannot take fewer than 0 coproduct iterations: %s < 0" % str(n))
            if n == 0:
                return self
            if n == 1:
                return self.coproduct()
            from sage.functions.all import floor, ceil
            from sage.rings.all import Integer

            # Use coassociativity of `\Delta` to perform many coproducts simultaneously.
            fn = floor(Integer(n-1)/2); cn = ceil(Integer(n-1)/2)
            split = lambda a,b: tensor([a.coproduct_iterated(fn), b.coproduct_iterated(cn)])
            return self.coproduct().apply_multilinear_morphism(split)

def level1_UpGj(p,klist,m):
    r"""
    Returns a list `[A_k]` of square matrices over ``IntegerRing(p^m)``
    parameterised by the weights k in ``klist``. The matrix `A_k` is the finite
    square matrix which occurs on input p,k and m in Step 6 of Algorithm 1 in
    [AGBL]_. Notational change from paper: In Step 1 following Wan we defined
    j by `k = k_0 + j(p-1)` with `0 \le k_0 < p-1`. Here we replace j by
    ``kdiv`` so that we may use j as a column index for matrices.

    INPUT:

    - ``p`` -- prime at least 5.
    - ``klist`` -- list of integers congruent modulo `(p-1)` (the weights).
    - ``m`` -- positive integer.

    OUTPUT:

    - list of square matrices.

    EXAMPLES::

        sage: from sage.modular.overconvergent.hecke_series import level1_UpGj
        sage: level1_UpGj(7,[100],5)
        [
        [    1   980  4802     0     0]
        [    0 13727 14406     0     0]
        [    0 13440  7203     0     0]
        [    0  1995  4802     0     0]
        [    0  9212 14406     0     0]
        ]
    """
    # Step 1
    t = cputime()

    k0 = klist[0] % (p-1)
    n = floor(((p+1)/(p-1)) * (m+1))
    ell = dimension_modular_forms(1, k0 + n*(p-1))
    ellp = ell*p
    mdash = m + ceil(n/(p+1))

    verbose("done step 1", t)
    t = cputime()
    # Steps 2 and 3

    e,Ep1 = katz_expansions(k0,p,ellp,mdash,n)

    verbose("done steps 2+3", t)
    t=cputime()
    # Step 4

    G = compute_G(p, Ep1)
    Alist = []

    verbose("done step 4a", t)
    t=cputime()
    for k in klist:
        k = ZZ(k) # convert to sage integer
        kdiv = k // (p-1)
        Gkdiv = G**kdiv
        u = []
        for i in xrange(0,ell):
            ei = e[i]
            ui = Gkdiv*ei
            u.append(ui)

        verbose("done step 4b", t)
        t = cputime()
        # Step 5 and computation of T in Step 6

        S = e[0][0].parent()
        T = matrix(S,ell,ell)

        for i in xrange(0,ell):
            for j in xrange(0,ell):
                T[i,j] = u[i][p*j]

        verbose("done step 5", t)
        t = cputime()
        # Step 6: solve T = AE using fact E is upper triangular.
        # Warning: assumes that T = AE (rather than pT = AE) has
        # a solution over Z/(p^mdash). This has always been the case in
        # examples computed by the author, see Note 3.1.

        A = matrix(S,ell,ell)
        verbose("solving a square matrix problem of dimension %s" % ell, t)

        for i in xrange(0,ell):
            Ti = T[i]
            for j in xrange(0,ell):
                ej = Ti.parent()([e[j][l] for l in xrange(0,ell)])
                lj = ZZ(ej[j])
                A[i,j] = S(ZZ(Ti[j])/lj)
                Ti = Ti - A[i,j]*ej

        Alist.append(MatrixSpace(Zmod(p**m),ell,ell)(A))
        verbose("done step 6", t)

    return Alist

def higher_level_UpGj(p,N,klist,m,modformsring,bound):
    r"""
    Returns a list ``[A_k]`` of square matrices over ``IntegerRing(p^m)``
    parameterised by the weights k in ``klist``. The matrix `A_k` is the finite
    square matrix which occurs on input p,k,N and m in Step 6 of Algorithm 2 in
    [AGBL]_. Notational change from paper: In Step 1 following Wan we defined
    j by `k = k_0 + j(p-1)` with `0 \le k_0 < p-1`. Here we replace j by
    ``kdiv`` so that we may use j as a column index for matrices.)

    INPUT:

    - ``p`` -- prime at least 5.
    - ``N`` -- integer at least 2 and not divisible by p (level).
    - ``klist`` -- list of integers all congruent modulo (p-1) (the weights).
    - ``m`` -- positive integer.
    - ``modformsring`` -- True or False.
    - ``bound`` -- (even) positive integer.

    OUTPUT:

    - list of square matrices.

    EXAMPLES::

        sage: from sage.modular.overconvergent.hecke_series import higher_level_UpGj
        sage: higher_level_UpGj(5,3,[4],2,true,6)
        [
        [ 1  0  0  0  0  0]
        [ 0  1  0  0  0  0]
        [ 0  7  0  0  0  0]
        [ 0  5 10 20  0  0]
        [ 0  7 20  0 20  0]
        [ 0  1 24  0 20  0]
        ]

    """
    t = cputime()
    # Step 1

    k0 = klist[0] % (p-1)
    n = floor(((p+1)/(p-1)) * (m+1))
    elldash = compute_elldash(p,N,k0,n)
    elldashp = elldash*p
    mdash = m + ceil(n/(p+1))

    verbose("done step 1",t)
    t = cputime()
    # Steps 2 and 3

    e,Ep1 = higher_level_katz_exp(p,N,k0,m,mdash,elldash,elldashp,modformsring,bound)
    ell = dimension(transpose(e)[0].parent())
    S = e[0,0].parent()

    verbose("done steps 2+3", t)
    t = cputime()
    # Step 4

    R = Ep1.parent()
    G = compute_G(p, Ep1)
    Alist = []

    verbose("done step 4a", t)
    t = cputime()
    for k in klist:
        k = ZZ(k) # convert to sage integer
        kdiv = k // (p-1)
        Gkdiv = G**kdiv

        T = matrix(S,ell,elldash)
        for i in xrange(ell):
            ei = R(e[i].list())
            Gkdivei = Gkdiv*ei; # act by G^kdiv
            for j in xrange(0, elldash):
                T[i,j] = Gkdivei[p*j]

        verbose("done steps 4b and 5", t)
        t = cputime()

        # Step 6: solve T = AE using fact E is upper triangular.
        # Warning: assumes that T = AE (rather than pT = AE) has
        # a solution over Z/(p^mdash). This has always been the case in
        # examples computed by the author, see Note 3.1.

        A = matrix(S,ell,ell)
        verbose("solving a square matrix problem of dimension %s" % ell)
        verbose("elldash is %s" % elldash)

        for i in xrange(0,ell):
            Ti = T[i]
            for j in xrange(0,ell):
                ej = Ti.parent()([e[j][l] for l in xrange(0,elldash)])
                ejleadpos = ej.nonzero_positions()[0]
                lj = ZZ(ej[ejleadpos])
                A[i,j] = S(ZZ(Ti[j])/lj)
                Ti = Ti - A[i,j]*ej

        Alist.append(MatrixSpace(Zmod(p**m),ell,ell)(A))
        verbose("done step 6", t)

    return Alist

def theta_by_cholesky(self, q_prec):
    r"""
    Uses the real Cholesky decomposition to compute (the `q`-expansion of) the
    theta function of the quadratic form as a power series in `q` with terms
    correct up to the power `q^{\text{q\_prec}}`. (So its error is `O(q^
    {\text{q\_prec} + 1})`.)

    REFERENCE:

        From Cohen's "A Course in Computational Algebraic Number Theory" book,
        p 102.

    EXAMPLES::

        ## Check the sum of 4 squares form against Jacobi's formula
        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1])
        sage: Theta = Q.theta_by_cholesky(10)
        sage: Theta
        1 + 8*q + 24*q^2 + 32*q^3 + 24*q^4 + 48*q^5 + 96*q^6 + 64*q^7 + 24*q^8 + 104*q^9 + 144*q^10
        sage: Expected =  [1] + [8*sum([d for d in divisors(n) if d%4 != 0])  for n in range(1,11)]
        sage: Expected
        [1, 8, 24, 32, 24, 48, 96, 64, 24, 104, 144]
        sage: Theta.list() == Expected
        True

    ::

        ## Check the form x^2 + 3y^2 + 5z^2 + 7w^2 represents everything except 2 and 22.
        sage: Q = DiagonalQuadraticForm(ZZ, [1,3,5,7])
        sage: Theta = Q.theta_by_cholesky(50)
        sage: Theta_list = Theta.list()
        sage: [m  for m in range(len(Theta_list))  if Theta_list[m] == 0]
        [2, 22]

    """
    ## RAISE AN ERROR -- This routine is deprecated!
    #raise NotImplementedError, "This routine is deprecated.  Try theta_series(), which uses theta_by_pari()."


    n = self.dim()
    theta = [0 for i in range(q_prec+1)]
    PS = PowerSeriesRing(ZZ, 'q')

    bit_prec = 53                                       ## TO DO: Set this precision to reflect the appropriate roundoff
    Cholesky = self.cholesky_decomposition(bit_prec)     ## error estimate, to be confident through our desired q-precision.
    Q = Cholesky      ##  <----  REDUNDANT!!!
    R = RealField(bit_prec)
    half = R(0.5)



    ## 1. Initialize
    i = n - 1
    T = [R(0)  for j in range(n)]
    U = [R(0)  for j in range(n)]
    T[i] = R(q_prec)
    U[i] = 0
    L = [0 for j in range (n)]
    x = [0 for j in range (n)]


    ## 2. Compute bounds
    #Z = sqrt(T[i] / Q[i,i])      ## IMPORTANT NOTE: sqrt preserves the precision of the real number it's given... which is not so good... =|
    #L[i] = floor(Z - U[i])       ## Note: This is a Sage Integer
    #x[i] = ceil(-Z - U[i]) - 1   ## Note: This is a Sage Integer too


    done_flag = False
    from_step4_flag = False
    from_step3_flag = True        ## We start by pretending this, since then we get to run through 2 and 3a once. =)

    #double Q_val_double;
    #unsigned long Q_val;                 // WARNING: Still need a good way of checking overflow for this value...



    ## Big loop which runs through all vectors
    while (done_flag == False):

        ## Loop through until we get to i=1 (so we defined a vector x)
        while from_step3_flag or from_step4_flag:              ## IMPORTANT WARNING:  This replaces a do...while loop, so it may have to be adjusted!

            ## Go to directly to step 3 if we're coming from step 4, otherwise perform step 2.
            if from_step4_flag:
                from_step4_flag = False
            else:
                ## 2. Compute bounds
                from_step3_flag = False
                Z = sqrt(T[i] / Q[i,i])
                L[i] = floor(Z - U[i])
                x[i] = ceil(-Z - U[i]) - 1



            ## 3a. Main loop

            ## DIAGNOSTIC
            #print
            #print "  L = ", L
            #print "  x = ", x
#.........这里部分代码省略.........

 def rangea(t):
     tp0=t*wprime0
     tp1=t*wprime1
     return sorted(list(union(range(floor(xmin-_sage_const_1 -tp0-wprime0),ceil(xmax-tp0)+_sage_const_1 ),range(floor(ymin-_sage_const_1 -tp1-wprime1),ceil(ymax-tp1)+_sage_const_1 ))),key=abs)