def symbolic_modulus_switching(n, alpha, q, h, m=None, epsilon=None):
    if epsilon is None:
        epsilon = var("epsilon")
        assume(epsilon>0)
        assume(epsilon<ZZ(1)/2)

    delta_0 = var("delta_0")
    assume(delta_0>=1.0)

    if m is None:
        m = sqrt(n*log(q)/log(delta_0))

    e = alpha*q/sqrt(2*pi)

    c = e * sqrt(m-n)/sqrt(h)

    v = delta_0**m * (q/c)**(n/m)  # norm of the vector
    v_ = v**2/m  # variance of each component

    v_r = (m-n) * e**2 *v_   # noise contribution
    v_l = h * v_ * c**2  # nose contribution of rounding noise

    # epsilon = exp(-pi*(|v|^2/q^2))

    f = log(1/epsilon)/pi == (v_l + v_r)/q**2

    # solve
    f = 2* q**2 * m * f * pi
    f = f.simplify_full()
    f = f.solve(delta_0**(2*m))[0]
    f = f.log().canonicalize_radical()
    f = f.solve(log(delta_0))[0]
    f = f.simplify_log()
    return f

def symbolic_sis(n, alpha, q, m=None, epsilon=None):
    if epsilon is None:
        epsilon = var("epsilon")
        assume(epsilon>0)
        assume(epsilon<ZZ(1)/2)
    delta_0 = var("delta_0")
    assume(delta_0>=1.0)

    e = alpha*q/sqrt(2*pi)

    if m is None:
        m = sqrt(n*log(q)/log(delta_0))

    v = e * delta_0**m * q**(n/m)  # norm of the vector

    # epsilon = exp(-pi*(|v|^2/q^2))
    f = log(1/epsilon)/pi == (v/q)**2

    # solve
    f = 2* q**2 * m * f * pi
    f = f.simplify_full()
    f = f.solve(delta_0**(2*m))[0]
    f = f.log().canonicalize_radical()
    f = f.solve(log(delta_0))[0]
    f = f.simplify_log()
    return f

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 quadratic_L_function__exact(n, d):
    r"""
    Returns the exact value of a quadratic twist of the Riemann Zeta function
    by `\chi_d(x) = \left(\frac{d}{x}\right)`.

    The input `n` must be a critical value.

    EXAMPLES::

        sage: quadratic_L_function__exact(1, -4)
        1/4*pi
        sage: quadratic_L_function__exact(-4, -4)
        5/2
        sage: quadratic_L_function__exact(2, 1)
        1/6*pi^2

    TESTS::

        sage: quadratic_L_function__exact(2, -4)
        Traceback (most recent call last):
        ...
        TypeError: n must be a critical value (i.e. odd > 0 or even <= 0)

    REFERENCES:

    - [Iwa1972]_, pp 16-17, Special values of `L(1-n, \chi)` and `L(n, \chi)`
    - [IR1990]_
    - [Was1997]_
    """
    from sage.all import SR, sqrt
    if n <= 0:
        return QuadraticBernoulliNumber(1-n,d)/(n-1)
    elif n >= 1:
        # Compute the kind of critical values (p10)
        if kronecker_symbol(fundamental_discriminant(d), -1) == 1:
            delta = 0
        else:
            delta = 1

        # Compute the positive special values (p17)
        if ((n - delta) % 2 == 0):
            f = abs(fundamental_discriminant(d))
            if delta == 0:
                GS = sqrt(f)
            else:
                GS = I * sqrt(f)
            ans = SR(ZZ(-1)**(1+(n-delta)/2))
            ans *= (2*pi/f)**n
            ans *= GS     # Evaluate the Gauss sum here! =0
            ans *= QQ.one()/(2 * I**delta)
            ans *= QuadraticBernoulliNumber(n,d)/factorial(n)
            return ans
        else:
            if delta == 0:
                raise TypeError("n must be a critical value (i.e. even > 0 or odd < 0)")
            if delta == 1:
                raise TypeError("n must be a critical value (i.e. odd > 0 or even <= 0)")

def crack_when_pq_close(n):
    t = Integer(ceil(sqrt(n)))
    while True:
        k = t**2 - n
        if k > 0:
            s = Integer(int(round(sqrt(t**2 - n))))
            if s**2 + n == t**2:
                return t + s, t - s
        t += 1

def x5_jacobi_pwsr(prec):
    mx = int(ceil(sqrt(8 * prec) / QQ(2)) + 1)
    mn = int(floor(-(sqrt(8 * prec) - 1) / QQ(2)))
    mx1 = int(ceil((sqrt(8 * prec + 1) - 1) / QQ(2)) + 1)
    mn1 = int(floor((-sqrt(8 * prec + 1) - 1) / QQ(2)))
    R = LaurentPolynomialRing(QQ, names="t")
    t = R.gens()[0]
    S = PowerSeriesRing(R, names="q1")
    q1 = S.gens()[0]
    eta_3 = sum([QQ(-1) ** n * (2 * n + 1) * q1 ** (n * (n + 1) // 2)
                 for n in range(mn1, mx1)]) + bigO(q1 ** (prec + 1))
    theta = sum([QQ(-1) ** n * q1 ** (((2 * n + 1) ** 2 - 1) // 8) * t ** (n + 1)
                 for n in range(mn, mx)])
    # ct = qexp_eta(ZZ[['q1']], prec + 1)
    return theta * eta_3 ** 3 * QQ(8) ** (-1)

def insert_EC_L_functions(start=1, end=100):
    curves = C.ellcurves.curves
    for N in range(start, end):
        print "Processing conductor", N
        sys.stdout.flush()
        query = curves.find({'conductor': N, 'number': 1})
        for curve in query:
            E = EllipticCurve([int(x) for x in curve['ainvs']])
            L = lc.Lfunction_from_elliptic_curve(E)
            first_zeros = L.find_zeros_via_N(curve['rank'] + 1)
            if len(first_zeros) > 1:
                if not first_zeros[-2] == 0:
                    print "problem"

            z = float(first_zeros[-1])

            Lfunction_data = {}
            Lfunction_data['first_zero'] = z
            Lfunction_data['description'] = 'Elliptic curve L-function for curve ' + str(curve['label'][:-1])
            Lfunction_data['degree'] = 2
            Lfunction_data['signature'] = [0, 1]
            Lfunction_data['eta'] = [(1.0, 0), ]

            Lfunction_data['level'] = N
            Lfunction_data['special'] = {'type': 'elliptic', 'label': curve['label'][:-1]}

            coeffs = []

            for k in range(1, 11):
                coeffs.append(CC(E.an(k) / sqrt(k)))

            Lfunction_data['coeffs'] = [(float(x.real()), float(x.imag())) for x in coeffs]

            Lfunctions.insert(Lfunction_data)

def is_quotient(M, sym, rank):
    symbol = GenusSymbol_global_ring(MatrixSpace(ZZ, rank, rank).one())
    symbol._local_symbols = [
        Genus_Symbol_p_adic_ring(p, syms) for p, syms in sym.iteritems()]
    s = get_symbol_string(symbol)
    print s
    N = FiniteQuadraticModule(s)
    t = N.order() / M.order()
    if not Integer(t).is_square():
        return False
    else:
        t = sqrt(t)
    for p in Integer(N.order()).prime_factors():
        if not N.signature(p) == M.signature(p):
            return False
    # if not N.signature() == M.signature():
    #    return False
    for G in N.subgroups():
        if G.order() == t and G.is_isotropic():
            Q = G.quotient()
            if Q.is_isomorphic(M):
                print Q
                return N
            else:
                del Q
    del N
    return False

def class_nr_pos_def_qf(D):
    r"""
    Compute the class number of positive definite quadratic forms.
    For fundamental discriminants this is the class number of Q(sqrt(D)),
    otherwise it is computed using: Cohen 'A course in Computational Algebraic Number Theory', p. 233
    """
    if D>0:
        return 0
    D4 = D % 4
    if D4 == 3 or D4==2:
        return 0
    K = QuadraticField(D)
    if is_fundamental_discriminant(D):
        return K.class_number()
    else:
        D0 = K.discriminant()
        Df = ZZ(D).divide_knowing_divisible_by(D0)
        if not is_square(Df):
            raise ArithmeticError,"DId not get a discrinimant * square! D={0} disc(D)={1}".format(D,D0)
        D2 = sqrt(Df)
        h0 = QuadraticField(D0).class_number()
        w0 = _get_w(D0)
        w = _get_w(D)
        #print "w,w0=",w,w0
        #print "h0=",h0
        h = 1
        for p in prime_divisors(D2):
            h = QQ(h)*(1-kronecker(D0,p)/QQ(p))
        #print "h=",h
        #print "fak=",
        h=QQ(h*h0*D2*w)/QQ(w0)
        return h

def _pell_solve_1(D,m): # m^2 < D
    root_d = Integer(floor(sqrt(D)))
    a = Integer(floor(root_d))
    P = Integer(0)
    Q = Integer(1)
    p = [Integer(1),Integer(a)]
    q = [Integer(0),Integer(1)]
    i = Integer(1)
    x0 = Integer(0)
    y0 = Integer(0)
    prim_sols = []
    test = Integer(0)
    while not (Q == 1 and i%2 == 1) or i == 1:
        test = p[i]**2 - D* (q[i]**2)
        if test == 1:
            x0 = p[i]
            y0 = q[i]
        test = (m/test)
        if is_square(test) and test >= 1:
            test = Integer(test)
            prim_sols.append((test*p[i],test*q[i]))
        i+=1
        P = a*Q - P
        Q = (D-P**2)/Q
        a = Integer(floor((P+root_d)/Q))
        p.append(a*p[i-1]+p[i-2])
        q.append(a*q[i-1]+q[i-2])
    return (x0,y0), prim_sols

def quadratic_L_function__exact(n, d):
    r"""
    Returns the exact value of a quadratic twist of the Riemann Zeta function
    by `\chi_d(x) = \left(\frac{d}{x}\right)`.  

    References:

    - Iwasawa's "Lectures on p-adic L-functions", p16-17, "Special values of
      `L(1-n, \chi)` and `L(n, \chi)`
    - Ireland and Rosen's "A Classical Introduction to Modern Number Theory"
    - Washington's "Cyclotomic Fields"

    EXAMPLES::

        sage: bool(quadratic_L_function__exact(1, -4) == pi/4)
        True

    """
    from sage.all import SR, sqrt
    if n<=0:
        k = 1-n
        return -QuadraticBernoulliNumber(k,d)/k
    elif n>=1:
        ## Compute the kind of critical values (p10)
        if kronecker_symbol(fundamental_discriminant(d), -1) == 1:
            delta = 0
        else:
            delta = 1        

        ## Compute the positive special values (p17)
        if ((n - delta) % 2 == 0):
            f = abs(fundamental_discriminant(d))
            if delta == 0:                            
                GS = sqrt(f)
            else:
                GS = I * sqrt(f)
            ans = SR(ZZ(-1)**(1+(n-delta)/2))
            ans *= (2*pi/f)**n
            ans *= GS     ## Evaluate the Gauss sum here! =0
            ans *= 1/(2 * I**delta)
            ans *= QuadraticBernoulliNumber(n,d)/factorial(n)
            return ans
        else:
            if delta == 0:
                raise TypeError, "n must be a critical value!\n" + "(I.e. even > 0 or odd < 0.)"
            if delta == 1:
                raise TypeError, "n must be a critical value!\n" + "(I.e. odd > 0 or even <= 0.)"

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 _pell_solve_2(D,m): # m^2 >= D
    prim_sols = []
    t,u = _pell_solve_1(D,1)[0]
    if m > 0:
        L = Integer(0)
        U = Integer(floor(sqrt(m*(t-1)/(2*D))))
    else:
        L = Integer(ceil(sqrt(-m/D)))
        U = Integer(floor(sqrt(-m*(t+1)/(2*D))))
    for y in range(L,U+1):
        y = Integer(y)
        x = (m + D*(y**2))
        if is_square(x):
            x = Integer(sqrt(x))
            prim_sols.append((x,y))
            if not ((-x*x - y*y*D) % m == 0 and (2*y*x) % m == 0): # (x,y) and (-x,y) are in different solution classes, so add both
                prim_sols.append((-x,y))
    return (t,u),prim_sols

def gamma__exact(n):
    """
    Evaluates the exact value of the gamma function at an integer or
    half-integer argument.

    EXAMPLES::

        sage: gamma__exact(4)
        6
        sage: gamma__exact(3)
        2
        sage: gamma__exact(2)
        1
        sage: gamma__exact(1)
        1

        sage: gamma__exact(1/2)
        sqrt(pi)
        sage: gamma__exact(3/2)
        1/2*sqrt(pi)
        sage: gamma__exact(5/2)
        3/4*sqrt(pi)
        sage: gamma__exact(7/2)
        15/8*sqrt(pi)

        sage: gamma__exact(-1/2)
        -2*sqrt(pi)
        sage: gamma__exact(-3/2)
        4/3*sqrt(pi)
        sage: gamma__exact(-5/2)
        -8/15*sqrt(pi)
        sage: gamma__exact(-7/2)
        16/105*sqrt(pi)

    """
    from sage.all import sqrt
    ## SANITY CHECK
    if (not n in QQ) or (denominator(n) > 2):
        raise TypeError, "Oops!  You much give an integer or half-integer argument."

    if (denominator(n) == 1):
        if n <= 0:
            return infinity
        if n > 0:
            return factorial(n-1) 
    else:
        ans = QQ(1)
        while (n != QQ(1)/2):
            if (n<0): 
                ans *= QQ(1)/n
                n = n + 1
            elif (n>0):
                n = n - 1
                ans *= n 
            
        ans *= sqrt(pi)
        return ans

 def renormalise_coefficients(self):
     """
     This turns a list of algebraically normalised coefficients
     as above into a list of automorphically normalised,
     i.e. s <-> 1-s
     """
     # this also turns them into floats and complex.
     for n in range(len(self.dirichlet_coefficients)):
         self.dirichlet_coefficients[n] /= sqrt(float(n+1)**self.motivic_weight)

def download_hecke_algebras_full_lists_gen(**args):
    label = str(args['orbit_label'])
    res = db.hecke_orbits.lucky({'orbit_label': label})
    mydate = time.strftime("%d %B %Y")
    if res is None:
        return "No generators available"
    lang = args['lang']
    c = download_comment_prefix[lang]
    mat_start = "Mat(" if lang == 'gp' else "Matrix("
    mat_end = "~)" if lang == 'gp' else ")"
    entry = lambda r: "".join([mat_start,str(r),mat_end])

    outstr = c + 'Hecke algebra for Gamma0(%s) and weight %s, orbit label %s. List of generators for the algebra. Downloaded from the LMFDB on %s. \n\n'%(res['level'], res['weight'], res['orbit_label'], mydate)
    outstr += download_assignment_start[lang] + '[\\\n'
    outstr += ",\\\n".join([entry([list(k) for k in matrix(sqrt(len(r)), sqrt(len(r)), r).rows()]) for r in [[int(i) for i in j] for j in res['Zbasis']] ])
    outstr += ']'
    outstr += download_assignment_end[lang]
    outstr += '\n'
    return outstr

def vect_to_sym(v):
    n = ZZ(round((-1+sqrt(1+8*len(v)))/2))
    M = matrix(n)
    k = 0
    for i in range(n):
        for j in range(i, n):
            M[i,j] = v[k]
            M[j,i] = v[k]
            k=k+1
    return [[int(M[i,j]) for i in range(n)] for j in range(n)]

    def init_dir_char(self, chi):
        """
        Initiate with a Web Dirichlet character.
        """
        self.original_object = [chi]
        chi = chi.chi.primitive_character()
        self.object_type = "dirichletcharacter"
        self.dim = 1
        self.motivic_weight = 0
        self.conductor = ZZ(chi.conductor())
        self.bad_semistable_primes = []
        self.bad_pot_good = self.conductor.prime_factors()
        if chi.is_odd():
            aa = 1
            bb = I
        else:
            aa = 0
            bb = 1
        self.sign = chi.gauss_sum_numerical() / (bb * float(sqrt(chi.modulus())) )
        # this has now type python complex. later we need a gp complex
        self.sign = ComplexField()(self.sign)
        self.mu_fe = [aa]
        self.nu_fe = []
        self.gammaV = [aa]
        self.langlands = True
        self.selfdual = (chi.multiplicative_order() <= 2)
        # rather than all(  abs(chi(m).imag) < 0.0001 for m in range(chi.modulus() ) )
        self.primitive = True
        self.set_dokchitser_Lfunction()
        self.set_number_of_coefficients()
        self.dirichlet_coefficients = [ chi(m) for m in range(self.numcoeff + 1) ]
        if self.selfdual:
            self.coefficient_type = 2
        else:
            self.coefficient_type = 3
        self.coefficient_period = chi.modulus()
        self.besancon_bound = 10000

        def eu(p):
            """
            local euler factor
            """
            if self.selfdual:
                K = QQ
            else:
                K = ComplexField()
            R = PolynomialRing(K, "T")
            T = R.gens()[0]
            if self.conductor % p != 0:
                return  1 - ComplexField()(chi(p)) * T
            else:
                return R(1)

        self.local_euler_factor = eu
        self.ld.gp().quit()

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