def _draw_funddom_d(coset_reps,format="MP",z0=I):
    r""" Draw a fundamental domain for self in the circle model
    INPUT:
    - ''format''  -- (default 'Disp') How to present the f.d.
    =  'S'  -- Display directly on the screen
    - z0          -- (default I) the upper-half plane is mapped to the disk by z-->(z-z0)/(z-z0.conjugate())
    EXAMPLES::
        

    sage: G=MySubgroup(Gamma0(3))
    sage: G._draw_funddom_d()
        
    """
    # The fundamental domain consists of copies of the standard fundamental domain
    pi=RR.pi()
    from sage.plot.plot import (Graphics,line)
    g=Graphics()
    bdcirc=_circ_arc(0 ,2 *pi,0 ,1 ,1000 )
    g=g+bdcirc
    # Corners
    x1=-RR(0.5) ; y1=RR(sqrt(3 )/2)
    x2=RR(0.5) ; y2=RR(sqrt(3 )/2)
    z_inf=1 
    l1 = _geodesic_between_two_points_d(x1,y1,x1,infinity)
    l2 = _geodesic_between_two_points_d(x2,y2,x2,infinity)
    c0 = _geodesic_between_two_points_d(x1,y1,x2,y2)
    tri=c0+l1+l2
    g=g+tri
    for A in coset_reps:
        [a,b,c,d]=A
        if(a==1  and b==0  and c==0  and d==1 ):
            continue
        if(a<0 ):
            a=-a; b=-b; c=-c; d=-1 
        if(c==0 ): # then this is easier
            l1 = _geodesic_between_two_points_d(x1+b,y1,x1+b,infinity)
            l2 = _geodesic_between_two_points_d(x2+b,y2,x2+b,infinity)
            c0 = _geodesic_between_two_points_d(x1+b,y1,x2+b,y2)
            # c0=line(L0); l1=line(L1); l2=line(L2); l3=line(L3)
            tri=c0+l1+l2
            g=g+tri
        else:
            den=(c*x1+d)**2 +c**2 *y1**2 
            x1_t=(a*c*(x1**2 +y1**2 )+(a*d+b*c)*x1+b*d)/den
            y1_t=y1/den
            den=(c*x2+d)**2 +c**2 *y2**2 
            x2_t=(a*c*(x2**2 +y2**2 )+(a*d+b*c)*x2+b*d)/den
            y2_t=y2/den
            inf_t=a/c
            c0=_geodesic_between_two_points_d(x1_t,y1_t,x2_t,y2_t)
            c1=_geodesic_between_two_points_d(x1_t,y1_t,inf_t,0.0 )
            c2=_geodesic_between_two_points_d(x2_t,y2_t,inf_t,0.0 )
            tri=c0+c1+c2
            g=g+tri
    g.xmax(1 )
    g.ymax(1 )
    g.xmin(-1 )
    g.ymin(-1 )
    g.set_aspect_ratio(1 )
    return g

def _draw_funddom(coset_reps,format="S"):
    r""" Draw a fundamental domain for G.
    
    INPUT:
    
    - ``format``  -- (default 'Disp') How to present the f.d.
    -   ``S`` -- Display directly on the screen
    
    EXAMPLES::        


    sage: G=MySubgroup(Gamma0(3))
    sage: G._draw_funddom()
        
    """
    pi=RR.pi()
    pi_3 = pi / RR(3.0)
    from sage.plot.plot import (Graphics,line)
    from sage.functions.trig import (cos,sin)
    g=Graphics()
    x1=RR(-0.5) ; y1=RR(sqrt(3 )/2 )
    x2=RR(0.5) ; y2=RR(sqrt(3 )/2 )
    xmax=RR(20.0) 
    l1 = line([[x1,y1],[x1,xmax]])
    l2 = line([[x2,y2],[x2,xmax]])
    l3 = line([[x2,xmax],[x1,xmax]]) # This is added to make a closed contour
    c0=_circ_arc(RR(pi/3.0) ,RR(2.0*pi)/RR(3.0) ,0 ,1 ,100 )
    tri=c0+l1+l3+l2
    g=g+tri
    for A in coset_reps:
        [a,b,c,d]=A
        if(a==1  and b==0  and c==0  and d==1 ):
            continue
        if(a<0 ):
            a=RR(-a); b=RR(-b); c=RR(-c); d=RR(-d) 
        else:
            a=RR(a); b=RR(b); c=RR(c); d=RR(d) 
        if(c==0 ): # then this is easier
            L0 = [[cos(pi_3*RR(i/100.0))+b,sin(pi_3*RR(i/100.0))] for i in range(100 ,201 )]
            L1 = [[x1+b,y1],[x1+b,xmax]]
            L2 = [[x2+b,y2],[x2+b,xmax]]
            L3 = [[x2+b,xmax],[x1+b,xmax]]
            c0=line(L0); l1=line(L1); l2=line(L2); l3=line(L3)
            tri=c0+l1+l3+l2
            g=g+tri
        else:
            den=(c*x1+d)**2 +c**2 *y1**2 
            x1_t=(a*c*(x1**2 +y1**2 )+(a*d+b*c)*x1+b*d)/den
            y1_t=y1/den
            den=(c*x2+d)**2 +c**2 *y2**2 
            x2_t=(a*c*(x2**2 +y2**2 )+(a*d+b*c)*x2+b*d)/den
            y2_t=y2/den
            inf_t=a/c
            c0=_geodesic_between_two_points(x1_t,y1_t,x2_t,y2_t)
            c1=_geodesic_between_two_points(x1_t,y1_t,inf_t,0. )
            c2=_geodesic_between_two_points(x2_t,y2_t,inf_t,0.0)
            tri=c0+c1+c2
            g=g+tri
    return g

def allbsd(line):
    r""" Parses one line from an allbsd file.  Returns the label and a
    dict containing fields with keys 'conductor', 'iso', 'number',
    'ainvs', 'rank', 'torsion', 'torsion_primes', 'tamagawa_product',
    'real_period', 'special_value', 'regulator', 'sha_an', 'sha',
    'sha_primes', all values being strings or floats or ints or lists
    of ints.

    Input line fields:

    conductor iso number ainvs rank torsion tamagawa_product real_period special_value regulator sha_an

    Sample input line:

    11 a 1 [0,-1,1,-10,-20] 0 5 5 1.2692093042795534217 0.25384186085591068434 1 1.00000000000000000000

    """
    data = split(line)
    label = data[0] + data[1] + data[2]
    ainvs = parse_ainvs(data[3])

    torsion = ZZ(data[5])
    sha_an = RR(data[10])
    sha = sha_an.round()
    sha_primes = sha.prime_divisors()
    torsion_primes = torsion.prime_divisors()

    data = {
        'conductor': int(data[0]),
        'iso': data[0] + data[1],
        'number': int(data[2]),
        'ainvs': ainvs,
        'rank': int(data[4]),
        'tamagawa_product': int(data[6]),
        'real_period': float(data[7]),
        'special_value': float(data[8]),
        'regulator': float(data[9]),
        'sha_an': float(sha_an),
        'sha':  int(sha),
        'sha_primes':  [int(p) for p in sha_primes],
        'torsion':  int(torsion),
        'torsion_primes':  [int(p) for p in torsion_primes]
        }

    return label, data

def _geodesic_between_two_points(x1, y1, x2, y2):
    r""" Geodesic path between two points hyperbolic upper half-plane

    INPUTS:

    - ''(x1,y1)'' -- starting point (0<y1<=infinity)
    - ''(x2,y2)'' -- ending point   (0<y2<=infinity)
    - ''z0''  -- (default I) the point in the upper corresponding
                 to the point 0 in the disc. I.e. the transform is
                 w -> (z-I)/(z+I)
    OUTPUT:

    - ''ca'' -- a polygonal approximation of a circular arc centered
    at c and radius r, starting at t0 and ending at t1


    EXAMPLES::


        sage: l=_geodesic_between_two_points(0.1,0.2,0.0,0.5)

    """
    pi = RR.pi()
    from sage.plot.plot import line
    from sage.functions.trig import arcsin

    # logging.debug("z1=%s,%s" % (x1,y1))
    # logging.debug("z2=%s,%s" % (x2,y2))
    if abs(x1 - x2) < 1e-10:
        # The line segment [x=x1, y0<= y <= y1]
        return line([[x1, y1], [x2, y2]])  # [0,0,x0,infinity]
    c = RR(y1 ** 2 - y2 ** 2 + x1 ** 2 - x2 ** 2) / RR(2 * (x1 - x2))
    r = RR(sqrt(y1 ** 2 + (x1 - c) ** 2))
    r1 = RR(y1 / r)
    r2 = RR(y2 / r)
    if abs(r1 - 1) < 1e-12:
        r1 = RR(1.0)
    elif abs(r2 + 1) < 1e-12:
        r2 = -RR(1.0)
    if abs(r2 - 1) < 1e-12:
        r2 = RR(1.0)
    elif abs(r2 + 1) < 1e-12:
        r2 = -RR(1.0)
    if x1 >= c:
        t1 = RR(arcsin(r1))
    else:
        t1 = RR(pi) - RR(arcsin(r1))
    if x2 >= c:
        t2 = RR(arcsin(r2))
    else:
        t2 = RR(pi) - arcsin(r2)
    # tmid = (t1 + t2) * RR(0.5)
    # a0 = min(t1, t2)
    # a1 = max(t1, t2)
    # logging.debug("c,r=%s,%s" % (c,r))
    # logging.debug("t1,t2=%s,%s"%(t1,t2))
    return _circ_arc(t1, t2, c, r)

def add_sha_tor_primes(N1,N2):
    """
    Add the 'sha', 'sha_primes', 'torsion_primes' fields to every
    curve in the database whose conductor is between N1 and N2
    inclusive.
    """
    query = {}
    query['conductor'] = { '$gte': int(N1), '$lte': int(N2) }
    res = curves.find(query)
    res = res.sort([('conductor', pymongo.ASCENDING)])
    n = 0
    for C in res:
        label = C['lmfdb_label']
        if n%1000==0: print label
        n += 1
        torsion = ZZ(C['torsion'])
        sha = RR(C['sha_an']).round()
        sha_primes = sha.prime_divisors()
        torsion_primes = torsion.prime_divisors()
        data = {}
        data['sha'] = int(sha)
        data['sha_primes'] = [int(p) for p in sha_primes]
        data['torsion_primes'] = [int(p) for p in torsion_primes]
        curves.update({'lmfdb_label': label}, {"$set": data}, upsert=True)

def _circ_arc(t0, t1, c, r, num_pts=500):
    r""" Circular arc
    INPUTS:
    - ''t0'' -- starting parameter
    - ''t1'' -- ending parameter
    - ''c''  -- center point of the circle
    - ''r''  -- radius of circle
    - ''num_pts''  -- (default 100) number of points on polygon
    OUTPUT:
    - ''ca'' -- a polygonal approximation of a circular arc centered
    at c and radius r, starting at t0 and ending at t1


    EXAMPLES::

        sage: ca=_circ_arc(0.1,0.2,0.0,1.0,100)

    """
    from sage.plot.plot import parametric_plot
    from sage.functions.trig import cos, sin
    from sage.all import var

    t00 = t0
    t11 = t1
    ## To make sure the line is correct we reduce all arguments to the same branch,
    ## e.g. [0,2pi]
    pi = RR.pi()
    while t00 < 0.0:
        t00 = t00 + RR(2.0 * pi)
    while t11 < 0:
        t11 = t11 + RR(2.0 * pi)
    while t00 > 2 * pi:
        t00 = t00 - RR(2.0 * pi)
    while t11 > 2 * pi:
        t11 = t11 - RR(2.0 * pi)

    xc = CC(c).real()
    yc = CC(c).imag()
    # L0 =
    # [[RR(r*cos(t00+i*(t11-t00)/num_pts))+xc,RR(r*sin(t00+i*(t11-t00)/num_pts))+yc]
    # for i in range(0 ,num_pts)]
    t = var("t")
    if t11 > t00:
        ca = parametric_plot((r * cos(t) + xc, r * sin(t) + yc), (t, t00, t11))
    else:
        ca = parametric_plot((r * cos(t) + xc, r * sin(t) + yc), (t, t11, t00))
    return ca

def draw_transformed_triangle_H(A,xmax=20):
    r"""
    Draw the modular triangle translated by A=[a,b,c,d]
    """
    #print "A=",A,type(A)
    pi=RR.pi()
    pi_3 = pi / RR(3.0)
    from sage.plot.plot import (Graphics,line)
    from sage.functions.trig import (cos,sin)
    x1=RR(-0.5) ; y1=RR(sqrt(3 )/2 )
    x2=RR(0.5) ; y2=RR(sqrt(3 )/2 )
    a,b,c,d = A #[0,0]; b=A[0,1]; c=A[1,0]; d=A[1,1]
    if a<0:
        a=RR(-a); b=RR(-b); c=RR(-c); d=RR(-d) 
    else:
        a=RR(a); b=RR(b); c=RR(c); d=RR(d) 
    if c==0: # then this is easier
        if a*d<>0:
            a=a/d; b=b/d; 
        L0 = [[a*cos(pi_3*RR(i/100.0))+b,a*sin(pi_3*RR(i/100.0))] for i in range(100 ,201 )]
        L1 = [[a*x1+b,a*y1],[a*x1+b,xmax]]
        L2 = [[a*x2+b,a*y2],[a*x2+b,xmax]]
        L3 = [[a*x2+b,xmax],[a*x1+b,xmax]]
        c0=line(L0); l1=line(L1); l2=line(L2); l3=line(L3)
        tri=c0+l1+l3+l2
    else:
        den=(c*x1+d)**2 +c**2 *y1**2 
        x1_t=(a*c*(x1**2 +y1**2 )+(a*d+b*c)*x1+b*d)/den
        y1_t=y1/den
        den=(c*x2+d)**2 +c**2 *y2**2 
        x2_t=(a*c*(x2**2 +y2**2 )+(a*d+b*c)*x2+b*d)/den
        y2_t=y2/den
        inf_t=a/c
        #print "A=",A
        #print "arg1=",x1_t,y1_t,x2_t,y2_t
        c0=_geodesic_between_two_points(x1_t,y1_t,x2_t,y2_t)
        #print "arg1=",x1_t,y1_t,inf_t
        c1=_geodesic_between_two_points(x1_t,y1_t,inf_t,0. )
        #print "arg1=",x2_t,y2_t,inf_t
        c2=_geodesic_between_two_points(x2_t,y2_t,inf_t,0.0)
        tri=c0+c1+c2
    return tri

def draw_funddom(coset_reps,format="S"):
    r""" Draw a fundamental domain for G.
    
    INPUT:
    
    - ``format``  -- (default 'Disp') How to present the f.d.
    -   ``S`` -- Display directly on the screen
    
    EXAMPLES::        


    sage: G=MySubgroup(Gamma0(3))
    sage: G._draw_funddom()
        
    """
    pi=RR.pi()
    pi_3 = pi / RR(3.0)
    from sage.plot.plot import (Graphics,line)
    from sage.functions.trig import (cos,sin)
    g=Graphics()
    x1=RR(-0.5) ; y1=RR(sqrt(3 )/2 )
    x2=RR(0.5) ; y2=RR(sqrt(3 )/2 )
    xmax=RR(20.0) 
    l1 = line([[x1,y1],[x1,xmax]])
    l2 = line([[x2,y2],[x2,xmax]])
    l3 = line([[x2,xmax],[x1,xmax]]) # This is added to make a closed contour
    c0=_circ_arc(RR(pi/3.0) ,RR(2.0*pi)/RR(3.0) ,0 ,1 ,100 )
    tri=c0+l1+l3+l2
    g=g+tri
    for A in coset_reps:
        if list(A)==[1,0,0,1]:
            continue

        tri=draw_transformed_triangle_H(A,xmax=xmax)
        g=g+tri
    return g

def kbarbar(weight):
    # The weight part of the analytic conductor
    return psi(RR(weight)/2).exp() / (2*RR.pi())

def _geodesic_between_two_points_d(x1,y1,x2,y2,z0=I):
    r""" Geodesic path between two points represented in the unit disc
         by the map w = (z-I)/(z+I)
    INPUTS:
    - ''(x1,y1)'' -- starting point (0<y1<=infinity)
    - ''(x2,y2)'' -- ending point   (0<y2<=infinity)
    - ''z0''  -- (default I) the point in the upper corresponding
                 to the point 0 in the disc. I.e. the transform is
                 w -> (z-I)/(z+I)
    OUTPUT:
    - ''ca'' -- a polygonal approximation of a circular arc centered
    at c and radius r, starting at t0 and ending at t1

    
    EXAMPLES::

        sage: l=_geodesic_between_two_points_d(0.1,0.2,0.0,0.5)
    
    """
    pi=RR.pi()
    from sage.plot.plot import line
    from sage.functions.trig import (cos,sin)    
    # First compute the points
    if(y1<0  or y2<0 ):
        raise ValueError,"Need points in the upper half-plane! Got y1=%s, y2=%s" %(y1,y2)
    if(y1==infinity):
        P1=CC(1 )
    else:
        P1=CC((x1+I*y1-z0)/(x1+I*y1-z0.conjugate()))
    if(y2==infinity):
        P2=CC(1 )
    else:
        P2=CC((x2+I*y2-z0)/(x2+I*y2-z0.conjugate()))
        # First find the endpoints of the completed geodesic in D
    if(x1==x2):
        a=CC((x1-z0)/(x1-z0.conjugate()))
        b=CC(1 )
    else:
        c=RR(y1**2 -y2**2 +x1**2 -x2**2 )/RR(2 *(x1-x2))
        r=RR(sqrt(y1**2 +(x1-c)**2 ))
        a=c-r
        b=c+r
        a=CC((a-z0)/(a-z0.conjugate()))
        b=CC((b-z0)/(b-z0.conjugate()))
    if( abs(a+b) < 1E-10 ): # On a diagonal
        return line([[P1.real(),P1.imag()],[P2.real(),P2.imag()]])
    th_a=a.argument()
    th_b=b.argument()
    # Compute the center of the circle in the disc model
    if( min(abs(b-1 ),abs(b+1 ))< 1E-10  and  min(abs(a-1 ),abs(a+1 ))>1E-10 ):
        c=b+I*(1 -b*cos(th_a))/sin(th_a)
    elif( min(abs(b-1 ),abs(b+1 ))> 1E-10  and  min(abs(a-1 ),abs(a+1 ))<1E-10 ):
        c=a+I*(1 -a*cos(th_b))/RR(sin(th_b))
    else:
        cx=(sin(th_b)-sin(th_a))/sin(th_b-th_a)
        c=cx+I*(1 -cx*cos(th_b))/RR(sin(th_b))
    # First find the endpoints of the completed geodesic
    r=abs(c-a)
    t1=CC(P1-c).argument()
    t2=CC(P2-c).argument()
    #print "t1,t2=",t1,t2
    return _circ_arc(t1,t2,c,r)