def saw_tooth_fn(x):
    if floor(x) == ceil(x):
        return 0
    elif x in QQ:
        return QQ(x)-QQ(floor(x))-QQ(1)/QQ(2)
    else:
        return x-floor(x)-0.5

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 tuples_even_wt_modular_forms(wt):
    '''
    Returns the list of tuples (p, q, r, s) such that
    4p + 6q + 10r +12s = wt.
    '''
    if wt < 0 or wt % 2 == 1:
        return []
    w = wt / 2
    return [(p, q, r, s) for p in range(0, floor(w / 2) + 1)
            for q in range(0, floor(w / 3) + 1)
            for r in range(0, floor(w / 5) + 1)
            for s in range(0, floor(w / 6) + 1)
            if 2 * p + 3 * q + 5 * r + 6 * s == w]

def make_curve(num_bits, num_curves=1): 
    """
    Description:
    
        Finds num_curves Barreto-Naehrig curves with a prime order that is at least 2^num_bits.
    
    Input:
    
        num_bits - number of bits for the prime order of the curve
        num_curves - number of curves to find
    
    Output:
    
        curves - list of the first num_curves BN curves each of prime order at least 2^num_bits;
                 each curve is represented as a tuple (q,t,r,k,D)
    
    """
    def P(y):
        x = Integer(y)
        return 36*pow(x,4) + 36*pow(x,3) + 24*pow(x,2) + 6*x + 1
    x = Integer(floor(pow(2, (num_bits)/4.0)/(sqrt(6))))
    q = 0
    r = 0
    t = 0
    curve_num = 0
    curves = []
    while curve_num < num_curves or (log(q).n()/log(2).n() < 2*num_bits and not (utils.is_suitable_q(q) and utils.is_suitable_r(r) and utils.is_suitable_curve(q,t,r,12,-3,num_bits))):
        t = Integer(6*pow(x,2) + 1)
        q = P(-x)
        r = q + 1 - t
        b = utils.is_suitable_q(q) and utils.is_suitable_r(r) and utils.is_suitable_curve(q,t,r,12,-3,num_bits)
        if b:
            try:
                assert floor(log(r)/log(2)) + 1 >= num_bits, 'Subgroup not large enough'  
                curves.append((q,t,r,12,-3))
                curve_num += 1
            except AssertionError as e:
                pass
        if curve_num < num_curves or not b:
            q = P(x)
            r = q+1-t
            if (utils.is_suitable_q(q) and utils.is_suitable_r(r) and utils.is_suitable_curve(q,t,r,12,-3,num_bits)):
                try:
                    assert floor(log(r)/log(2)) + 1 >= num_bits, 'Subgroup not large enough'  
                    curves.append((q,t,r,12,-3))
                    curve_num += 1
                except AssertionError as e:
                    pass  
        x += 1
    return curves

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 exact_cm_at_i_level_1(self, N=10,insert_in_db=True):
        r"""
        Use formula by Zagier (taken from pari implementation by H. Cohen) to compute the geodesic expansion of self at i
        and evaluate the constant term.

        INPUT:
        -''N'' -- integer, the length of the expansion to use.
        """
        try:
            [poldeg, monomials, X] = self.as_polynomial_in_E4_and_E6()
        except:
            return ""
        k = self.weight()
        tab = dict()
        QQ['x']
        tab[0] = 0 * x ** 0
        tab[1] = X[0] * x ** poldeg
        for ix in range(1, len(X)):
            tab[1] = tab[1] + QQ(X[ix]) * x ** monomials[ix][1]
        for n in range(1, N + 1):
            tmp = -QQ(k + 2 * n - 2) / QQ(12) * x * tab[n] + (x ** 2 - QQ(1)) / QQ(2) * ((tab[
                                                                                          n]).derivative())
            tab[n + 1] = tmp - QQ((n - 1) * (n + k - 2)) / QQ(144) * tab[n - 1]
        res = 0
        for n in range(1, N + 1):
            term = (tab[n](x=0)) * 12 ** (floor(QQ(n - 1) / QQ(2))) * x ** (n - 1) / factorial(n - 1)
            res = res + term
        
        return res

    def dirichlet_series_coeffs(self, prec, eps=1e-10):
        """
        Return the coefficients of the Dirichlet series representation
        of self, up to the given precision.

        INPUT:
           - prec -- positive integer
           - eps -- None or a positive real; any coefficient with absolute
             value less than eps is set to 0.
        """
        # Use multiplicativity to compute the Dirichlet series
        # coefficients, then make a DirichletSeries object.
        zero = RDF(0)
        coeffs = [RDF(0),RDF(1)] + [None]*(prec-2)

        from sage.all import log, floor   # TODO: slow
        
        # prime-power indexed coefficients
        for p in prime_range(2, prec):
            B = floor(log(prec, p)) + 1
            series = self._local_series(p, B)
            p_pow = p
            for i in range(1, B):
                coeffs[p_pow] = series[i] if (eps is None or abs(series[i])>eps) else zero
                p_pow *= p

        # non-prime-powers
        from sage.all import factor
        for n in range(2, prec):
            if coeffs[n] is None:
                a = prod(coeffs[p**e] for p, e in factor(n))
                coeffs[n] = a if (eps is None or abs(a) > eps) else zero

        return coeffs

 def create_small_record(self, min_prec=10, want_prec=100, max_length = 5242880, max_height_qexp = default_max_height):
     ### creates a duplicate record (fs) of this webnewform
     ### with lower precision to load faster on the web
     ### we aim to have at most max_length bytes
     ### but at least min_prec coefficients and we desire to have want_prec
     if min_prec>=self.prec:
         raise ValueError("Need higher precision, self.prec = {}".format(self.prec))
     if not hasattr(self, '_file_record_length'):
         self.update_from_db()
     l = self._file_record_length
         
     if l > max_length or self.prec > want_prec:
         nl = float(l)/float(self.prec)*float(want_prec)
         if nl > max_length:
             prec = max([floor(float(self.prec)/float(l)*float(max_length)), min_prec])
         else:
             prec = want_prec
         emf_logger.debug("Creating a new record with prec = {}".format(prec))
         self.prec = prec
         include_coeffs = self.complexity_of_first_nonvanishing_coefficients() <= default_max_height
         if include_coeffs:
             self.q_expansion = self.q_expansion.truncate_powerseries(prec)
             self._coefficients = {n:c for n,c in self._coefficients.iteritems() if n<prec}
         else:
             self.q_expansion = self.q_expansion.truncate_powerseries(1)
             self._coefficients = {}
             self.prec = 0
             self.coefficient_field = NumberField(self.absolute_polynomial, names=str(self.coefficient_field.gen()))
         self._embeddings['values'] = {n:c for n,c in self._embeddings['values'].iteritems() if n<prec}
         self._embeddings['prec'] = prec
         self.save_to_db()

 def cutout_digits(elt):
     digits = 1 if elt == 0 else floor(RR(abs(elt)).log(10)) + 1
     if digits > bigint_cutoff:
         # a large number would be replaced by ab...cd
         return digits - 7
     else:
         return 0

    def __init__(self,G,k=QQ(1)/QQ(2),number=0,ch=None,dual=False,version=1,dimension=1,**kwargs):
        r"""
        Initialize the Eta multiplier system: $\nu_{\eta}^{2(k+r)}$.
        INPUT:

        - G -- Group
        - ch -- character
        - dual -- if we have the dual (in this case conjugate)
        - weight -- Weight (recall that eta has weight 1/2 and eta**2k has weight k. If weight<>k we adjust the power accordingly.
        - number -- we consider eta^power (here power should be an integer so as not to change the weight...)
                
        """
        self._weight=QQ(k)
        if floor(self._weight-QQ(1)/QQ(2))==ceil(self._weight-QQ(1)/QQ(2)):
            self._half_integral_weight=1
        else:
            self._half_integral_weight=0
        MultiplierSystem.__init__(self,G,character=ch,dual=dual,dimension=dimension)
        number = number % 12
        if not is_even(number):
            raise ValueError,"Need to have v_eta^(2(k+r)) with r even!"
        self._pow=QQ((self._weight+number)) ## k+r
        self._k_den=self._pow.denominator()
        self._k_num=self._pow.numerator()
        self._K = CyclotomicField(12*self._k_den)
        self._z = self._K.gen()**self._k_num
        self._i = CyclotomicField(4).gen()
        self._fak = CyclotomicField(2*self._k_den).gen()**-self._k_num
        self._version = version
        
        self.is_consistent(k) # test consistency

def get_all_combis(g,n):
    dim = 3*g-3 + n
    reducible_boundaries = 0
    marks = range(1,n+1)
    if n != 0:
        first_mark_list = [marks.pop()]            
        for g1 in range(0, g + 1):
            for p in subsets(marks):
                r_marks = set(first_mark_list + p)
                if 3*g1 - 3 + len(r_marks) + 1 >= 0 and 3*(g-g1) - 3 + n - len(r_marks) + 1 >= 0:
                    reducible_boundaries+=1  
    
    else: #self.n == 0
        for g1 in range(1, floor(g/2.0)+1):
            reducible_boundaries+=1

    #print "computed red bound"
    indexes = range(1,n+dim+1) + range(n+dim+g+1, n+dim+g+reducible_boundaries + 2)
    codims = [1]*n + range(1,dim+1) + [1]*(reducible_boundaries +1)
    
    for w in WeightedIntegerVectors(dim,codims):
        combi = []
        #print w
        for index, wi in zip(indexes,w):
            combi += [index]*wi
        yield combi

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 as_polynomial_in_E4_and_E6(self,insert_in_db=True):
     r"""
     If self is on the full modular group writes self as a polynomial in E_4 and E_6.
     OUTPUT:
     -''X'' -- vector (x_1,...,x_n)
     with f = Sum_{i=0}^{k/6} x_(n-i) E_6^i * E_4^{k/4-i}
     i.e. x_i is the coefficient of E_6^(k/6-i)*
     """
     if(self.level() != 1):
         raise NotImplementedError("Only implemented for SL(2,Z). Need more generators in general.")
     if(self._as_polynomial_in_E4_and_E6 is not None and self._as_polynomial_in_E4_and_E6 != ''):
         return self._as_polynomial_in_E4_and_E6
     d = self._parent.dimension_modular_forms()  # dimension of space of modular forms
     k = self.weight()
     K = self.base_ring()
     l = list()
     # for n in range(d+1):
     #    l.append(self._f.q_expansion(d+2)[n])
     # v=vector(l) # (self._f.coefficients(d+1))
     v = vector(self.coefficients(range(d),insert_in_db=insert_in_db))
     d = dimension_modular_forms(1, k)
     lv = len(v)
     if(lv < d):
         raise ArithmeticError("not enough Fourier coeffs")
     e4 = EisensteinForms(1, 4).basis()[0].q_expansion(lv + 2)
     e6 = EisensteinForms(1, 6).basis()[0].q_expansion(lv + 2)
     m = Matrix(K, lv, d)
     lima = floor(k / 6)  # lima=k\6;
     if((lima - (k / 2)) % 2 == 1):
         lima = lima - 1
     poldeg = lima
     col = 0
     monomials = dict()
     while(lima >= 0):
         deg6 = ZZ(lima)
         deg4 = (ZZ((ZZ(k / 2) - 3 * lima) / 2))
         e6p = (e6 ** deg6)
         e4p = (e4 ** deg4)
         monomials[col] = [deg4, deg6]
         eis = e6p * e4p
         for i in range(1, lv + 1):
             m[i - 1, col] = eis.coefficients()[i - 1]
         lima = lima - 2
         col = col + 1
     if (col != d):
         raise ArithmeticError("bug dimension")
     # return [m,v]
     if self._verbose > 0:
         wmf_logger.debug("m={0}".format(m, type(m)))
         wmf_logger.debug("v={0}".format(v, type(v)))
     try:
         X = m.solve_right(v)
     except:
         return ""
     self._as_polynomial_in_E4_and_E6 = [poldeg, monomials, X]
     return [poldeg, monomials, X]

def orbit_label(j):
    x = AlphabeticStrings().gens()
    if(j < 26):
        label = str(x[j]).lower()
    else:
        j1 = j % 26
        j2 = floor(QQ(j) / QQ(26))
        label = str(x[j1]).lower()
        label = label + str(j2)
    return label

def euler_p_factor(root_list, PREC):
    ''' computes the coefficients of the pth Euler factor expanded as a geometric series
      ax^n is the Dirichlet series coefficient p^(-ns)
    '''
    PREC = floor(PREC)
    # return satake_list
    R = LaurentSeriesRing(CF, 'x')
    x = R.gens()[0]
    ep = prod([1 / (1 - a * x) for a in root_list])
    return ep + O(x ** (PREC + 1))

    def MgnLb_class(self,index):
        """
        Returns the class corresponding to the index from Carl Faber's ``MgnLb`` Maple program.
        This is useful for testing purposes. ::

            sage: from strataalgebra import *
            sage: s = StrataAlgebra(QQ,1,(1,2))
            sage: s.MgnLb_class(1)
            ps1
            sage: s.MgnLb_class(4)
            ka2
            sage: s.MgnLb_class(6)
            1/2*D_irr
            sage: s.MgnLb_class(2)
            ps2

        """
        #print "making classes again!"
        if index <= len(self.markings):
            return self.psi(index)
        index -= len(self.markings)
        if index <= self.moduli_dim:
            return self.kappa(index)
        index -= self.moduli_dim
        if index <= self.g:
            raise Exception("We don't do the ch classes!")
        index -= self.g
        if index == 1:
            return self.irr()
        index -=1
        
        marks = set(self.markings)
        reducible_boundaries = []
        if len(self.markings) != 0:
            first_mark_list = [marks.pop()]            
            for g1 in range(0, self.g + 1):
                for p in subsets(marks):
                    r_marks = set(first_mark_list + p)
                    if 3*g1 - 3 + len(r_marks) + 1 >= 0 and 3*(self.g-g1) - 3 + len(self.markings) - len(r_marks) + 1 >= 0:
                        reducible_boundaries.append( (g1, r_marks) )  
                        
            reducible_boundaries.sort(key = lambda b: sorted(list(b[1])))
            reducible_boundaries.sort(key = lambda b: len(b[1]))
            reducible_boundaries.sort(key = lambda b: b[0])
        
        else: #self.n == 0
            for g1 in range(1, floor(self.g/2.0)+1):
                reducible_boundaries.append( (g1, [])) 
            
        return self.boundary(*reducible_boundaries[index-1])   

def run(num_bits,k):
    """
    Description:
    
        Runs the Dupont-Enge-Morain method multiple times until a valid curve is found
    
    Input:
    
        num_bits - number of bits
        k - an embedding degree
    
    Output:
    
        (q,t,r,k,D) - an elliptic curve;
                      if no curve is found, the algorithm returns (0,0,0,k,0)
    
    """
    j,r,q,t = 0,0,0,0
    num_generates = 512
    h = num_bits/(euler_phi(k))
    tried = [(0,0)] # keep track of random values tried for efficiency
    for i in range(0,num_generates):
        D = 0
        y = 0
        while (D,y) in tried: # find a pair that we have not tried
            D = -randint(1, 1024) # pick a small D so that the CM method is fast
            D = fundamental_discriminant(D)
            m = 0.5*(h - log(-D).n()/(2*log(2)).n())
            if m < 1:
                m = 1
            y = randint(floor(2**(m-1)), floor(2**m))
        tried.append((D,y))
        q,t,r,k,D = method(num_bits,k,D,y) # run DEM
        if q != 0 and t != 0 and r != 0 and k != 0 and D != 0: # found an answer, so output it
            assert is_valid_curve(q,t,r,k,D), 'Invalid output'
            return q,t,r,k,D
    return 0,0,0,k,0 # found nothing

def HasFinitePointAt(F,p,c):
    # Tests whether y²=c*F(x) has a finite Qp-point with x and y both in Zp,
    # assuming that deg F = 6 and F integral
    Fp = GF(p)
    if p > 2 and Fp(F.leading_coefficient()): # Tests to accelerate case of large p
        # Write F(x) = c*lc(F)*R(x)²*S(x) mod p, R as big as possible
        R = F.parent()(1)
        S = F.parent()(1)
        for X in (F.base_extend(Fp)/Fp(F.leading_coefficient())).squarefree_decomposition():
            [G,v] = X # Term of the form G(x)^v in the factorisation of F(x) mod p
            S *= G**(v%2)
            R *= G**(v//2)
        r = R.degree()
        s = S.degree()
        if s == 0: # F(x) = C*R(x)² mod p, C = c*lc(F) constant
            if IsSquareInQp(c*F.leading_coefficient(),p):
                if p>r:# Then there is x s.t. R(x) nonzero and C is a square
                    return true
            #else: # C nonsquare, so if we have a Zp-point it must have R(x) = 0 mod p
            #Z = R.roots()
            ##TODO
        else:
            g = S.degree()//2 - 1 # genus of the curve y²=C*S(x)
            B = floor(p-1-2*g*sqrt(p)) # lower bound on number of points on y²=C*S(x) not at infty
            if B > r+s: # Then there is a point on y²=C*S(x) not at infty and with R(x) and S(x) nonzero
                return true
    #Now p is small, we can run a naive search
    q = p
    Z = []
    if p == 2:
        q = 8
    for x in range(q):
        y = F(x)
        # If we have found a root, save it (take care of the case p=2!)
        if (p > 2 or x < 2) and Fp(y) == 0:
            Z.append(x)
        # If we have a mod p point with y nonzero mod p, then it lifts, so we're done
        if IsSquareInQp(c*y, p):
            return true
    #So now, if we have a Qp-point, then its y-coordinate must be 0 mod p
    t = F.variables()[0]
    for z in Z:
        F1 = F(z+p*t)
        c1 = F1.content()
        F1 //= c1
        if HasFinitePointAt(F1,p,(c*c1).squarefree_part()):
            return true
    return false

 def __init__(self,group,dchar=(0,0),dual=False,is_trivial=False,dimension=1,**kwargs):
     if not ZZ(4).divides(group.level()):
         raise ValueError," Need level divisible by 4. Got:%s " % self._group.level()
     MultiplierSystem.__init__(self,group,dchar=dchar,dual=dual,is_trivial=is_trivial,dimension=dimension,**kwargs)
     self._i = CyclotomicField(4).gen()
     self._one = self._i**4
     self._weight= QQ(kwargs.get("weight",QQ(1)/QQ(2)))
     ## We have to make sure that we have the correct multiplier & character 
     ## for the desired weight
     if self._weight<>None:
         if floor(2*self._weight)<>ceil(2*self._weight):
             raise ValueError," Use ThetaMultiplier for half integral or integral weight only!"
         t = self.is_consistent(self._weight)
         if not t:
             self.set_dual()
             t1 = self.is_consistent(self._weight)    
             if not t1:
                 raise ArithmeticError,"Could not find consistent theta multiplier! Try to add a character."

def list_zeros(N=None,
               t=None,
               limit=None,
               fmt=None,
               download=None):
    if N is None:
        N = request.args.get("N", None, int)
    if t is None:
        t = request.args.get("t", 0, float)
    if limit is None:
        limit = request.args.get("limit", 100, int)
    if fmt is None:
        fmt = request.args.get("format", "plain")
    if download is None:
        download = request.args.get("download", "no")

    if limit < 0:
        limit = 100
    if N is not None:  # None is < 0!! WHAT THE WHAT!
        if N < 0:
            N = 0
    if t < 0:
        t = 0

    if limit > 100000:
        # limit = 100000
        #
        bread = [("L-functions", url_for("l_functions.l_function_top_page")),("Zeros of $\zeta(s)$", url_for(".zetazeros"))]
        return render_template('single.html', title="Too many zeros", bread=bread, kid = "dq.zeros.zeta.toomany")

    if N is not None:
        zeros = zeros_starting_at_N(N, limit)
    else:
        zeros = zeros_starting_at_t(t, limit)

    if fmt == 'plain':
        response = flask.Response(("%d %s\n" % (n, nstr(z,31+floor(log(z,10))+1,strip_zeros=False,min_fixed=-inf,max_fixed=+inf)) for (n, z) in zeros))
        response.headers['content-type'] = 'text/plain'
        if download == "yes":
            response.headers['content-disposition'] = 'attachment; filename=zetazeros'
    else:
        response = str(list(zeros))

    return response