def number_of_coefficients_needed(Q, kappa_fe, lambda_fe, max_t):
    # TODO: This doesn't work. Trouble when computing t0
    # We completely mimic what lcalc does when it decides whether
    # to print a warning.

    DIGITS = 14    # These are names of lcalc parameters, and we are
    DIGITS2 = 2    # mimicking them.

    logger.debug("Start NOC")
    theta = sum(kappa_fe)
    c = DIGITS2 * log(10.0)
    a = len(kappa_fe)

    c1 = 0.0
    for j in range(a):
        logger.debug("In loop NOC")
        t0 = kappa_fe[j] * max_t + complex(lambda_fe[j]).imag()
        logger.debug("In loop 2 NOC")
        if abs(t0) < 2 * c / (math.pi * a):
            logger.debug("In loop 3_1 NOC")
            c1 += kappa_fe[j] * pi / 2.0
        else:
            c1 += kappa_fe[j] * abs(c / (t0 * a))
            logger.debug("In loop 3_2 NOC")

    return int(round(Q * exp(log(2.3 * DIGITS * theta / c1) * theta) + 10))

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 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 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 central_char_function(self):
        dim = self.dimension()
        dfactor = (-1) ** dim
        # doubling insures integers below
        # we could test for when we need it, but then we carry the "if"
        # throughout
        charf = 2 * self.character_field()
        localfactors = self.local_factors_table()
        bad = [0 if dim + 1 > len(z) else 1 for z in localfactors]
        localfactors = [self.from_conjugacy_class_index_to_polynomial(j + 1) for j in range(len(localfactors))]
        localfactors = [z.leading_coefficient() * dfactor for z in localfactors]
        # Now take logs to figure out what power these are
        mypi = RealField(100)(pi)
        localfactors = [charf * log(z) / (2 * I * mypi) for z in localfactors]
        localfactorsa = [z.real().round() % charf for z in localfactors]
        # Test to see if we are ok?
        localfactorsa = [localfactorsa[j] if bad[j] > 0 else -1 for j in range(len(localfactorsa))]

        def myfunc(inp, n):
            fn = list(factor(inp))
            pvals = [[localfactorsa[self.any_prime_to_cc_index(z[0]) - 1], z[1]] for z in fn]
            # -1 is the marker that the prime divides the conductor
            for j in range(len(pvals)):
                if pvals[j][0] < 0:
                    return -1
            pvals = sum([z[0] * z[1] for z in pvals])
            return pvals % n

        return myfunc

def _number_of_bits(n):
    """
    Description:
        
        Returns the number of bits in the binary representation of n
    
    Input:
    
        n - integer
    
    Output:
    
        num_bits - number of bits
    
    """    
    if n == 0:
        return 1
    else:
        return floor(log(n).n()/log(2).n()) + 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 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

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 dimension(self,k,ignore=False, debug = 0):
        if k < 2 and not ignore:
            raise NotImplementedError("k has to >= 2")
        s = self._signature
        if not (2*k in ZZ):
            raise ValueError("k has to be integral or half-integral")
        if (2*k+s)%4 != 0 and not ignore:
            raise NotImplementedError("2k has to be congruent to -signature mod 4")
        if self._v2.has_key(0):
            v2 = self._v2[0]
        else:
            v2 = 1

        if self._g != None:
            if not self._aniso_formula:
                vals = self._g.values()
            #else:
                #print "using aniso_formula"
            M = self._g
        else:
            vals = self._M.values()
            M = self._M
            
        prec = ceil(max(log(M.order(),2),52)+1)+17
        #print prec
        RR = RealField(prec)
        CC = ComplexField(prec)
        d = self._d
        m = self._m
        if debug > 0: print d,m
            
        if self._alpha3 == None:
            if self._aniso_formula:
                self._alpha4 = 1
                self._alpha3 = -sum([BB(a)*mm for a,mm in self._v2.iteritems() if a != 0])
                #print self._alpha3
                self._alpha3 += Integer(d) - Integer(1) - self._g.beta_formula()
                #print self._alpha3, self._g.a5prime_formula()
                self._alpha3 = self._alpha3/RR(2)
            else:
                self._alpha3 = sum([(1-a)*mm for a,mm in self._v2.iteritems() if a != 0])
                #print self._alpha3
                self._alpha3 += sum([(1-a)*mm for a,mm in vals.iteritems() if a != 0])
                #print self._alpha3
                self._alpha3 = self._alpha3 / Integer(2)
                self._alpha4 = 1/Integer(2)*(vals[0]+v2) # the codimension of SkL in MkL
        alpha3 = self._alpha3
        alpha4 = self._alpha4
        if debug > 0: print alpha3, alpha4
        g1=M.char_invariant(1)
        g1=CC(g1[0]*g1[1])
        #print g1
        g2=M.char_invariant(2)
        g2=RR(real(g2[0]*g2[1]))
        if debug > 0: print g2, g2.parent()
        g3=M.char_invariant(-3)
        g3=CC(g3[0]*g3[1])
        if debug > 0: print RR(d) / RR(4), sqrt(RR(m)) / RR(4), CC(exp(2 * CC.pi() * CC(0,1) * (2 * k + s) / Integer(8)))
        alpha1 = RR(d) / RR(4) - (sqrt(RR(m)) / RR(4)  * CC(exp(2 * CC.pi() * CC(0,1) * (2 * k + s) / Integer(8))) * g2)
        if debug > 0: print alpha1
        alpha2 = RR(d) / RR(3) + sqrt(RR(m)) / (3 * sqrt(RR(3))) * real(exp(CC(2 * CC.pi() * CC(0,1) * (4 * k + 3 * s - 10) / 24)) * (g1+g3))
        if debug > 0: print alpha1, alpha2, g1, g2, g3, d, k, s
        dim = real(d + (d * k / Integer(12)) - alpha1 - alpha2 - alpha3)
        if debug > 0:
            print "dimension:", dim
        if abs(dim-round(dim)) > 1e-6:
            raise RuntimeError("Error ({0}) too large in dimension formula for {1} and k={2}".format(abs(dim-round(dim)), self._M if self._M is not None else self._g, k))
        dimr = dim
        dim = Integer(round(dim))
        if k >=2 and dim < 0:
            raise RuntimeError("Negative dimension (= {0}, {1})!".format(dim, dimr))
        return dim

def silke(A, c, beta, h, m=None, scale=1, float_type="double"):
    """

    :param A:    LWE matrix
    :param c:    LWE vector
    :param beta: BKW block size
    :param m:    number of samples to consider
    :param scale: scale rhs of lattice by this factor

    """
    from fpylll import BKZ, IntegerMatrix, LLL, GSO
    from fpylll.algorithms.bkz2 import BKZReduction as BKZ2

    if m is None:
        m = A.nrows()

    L = dual_instance1(A, scale=scale)
    L = IntegerMatrix.from_matrix(L)
    L = LLL.reduction(L, flags=LLL.VERBOSE)
    M = GSO.Mat(L, float_type=float_type)
    bkz = BKZ2(M)
    t = 0.0
    param = BKZ.Param(block_size=beta,
                      strategies=BKZ.DEFAULT_STRATEGY,
                      auto_abort=True,
                      max_loops=16,
                      flags=BKZ.VERBOSE|BKZ.AUTO_ABORT|BKZ.MAX_LOOPS)
    bkz(param)
    t += bkz.stats.total_time

    H = copy(L)

    import pickle
    pickle.dump(L, open("L-%d-%d.sobj"%(L.nrows, beta), "wb"))

    E = []
    Y = set()
    V = set()
    y_i = vector(ZZ, tuple(L[0]))
    Y.add(tuple(y_i))
    E.append(apply_short1(y_i, A, c, scale=scale)[1])

    v = L[0].norm()
    v_ = v/sqrt(L.ncols)
    v_r = 3.2*sqrt(L.ncols - A.ncols())*v_/scale
    v_l = sqrt(h)*v_

    fmt = u"{\"t\": %5.1fs, \"log(sigma)\": %5.1f, \"log(|y|)\": %5.1f, \"log(E[sigma]):\" %5.1f}"

    print
    print fmt%(t,
               log(abs(E[-1]), 2),
               log(L[0].norm(), 2),
               log(sqrt(v_r**2 + v_l**2), 2))
    print
    for i in range(m):
        t = cputime()
        M = GSO.Mat(L, float_type=float_type)
        bkz = BKZ2(M)
        t = cputime()
        bkz.randomize_block(0, L.nrows, stats=None, density=3)
        LLL.reduction(L)
        y_i = vector(ZZ, tuple(L[0]))
        l_n = L[0].norm()
        if L[0].norm() > H[0].norm():
            L = copy(H)
        t = cputime(t)

        Y.add(tuple(y_i))
        V.add(y_i.norm())
        E.append(apply_short1(y_i, A, c, scale=scale)[1])
        if len(V) >= 2:
            fmt =  u"{\"i\": %4d, \"t\": %5.1fs, \"log(|e_i|)\": %5.1f, \"log(|y_i|)\": %5.1f,"
            fmt += u"\"log(sigma)\": (%5.1f,%5.1f), \"log(|y|)\": (%5.1f,%5.1f), |Y|: %5d}"
            print fmt%(i+2, t, log(abs(E[-1]), 2), log(l_n, 2), log_mean(E), log_var(E), log_mean(V), log_var(V), len(Y))

    return E

def log_mean(X):
    return log(mean([abs(x) for x in X]), 2)

    def dimension(self,k,ignore=False, debug = 0):
        if k < 2 and not ignore:
            raise NotImplementedError("k has to >= 2")
        s = self._signature
        if not (2*k in ZZ):
            raise ValueError("k has to be integral or half-integral")
        if (2*k+s)%2 != 0:
            return 0
        m = self._m
        n2 = self._n2
        if self._v2.has_key(0):
            v2 = self._v2[0]
        else:
            v2 = 1

        if self._g != None:
            if not self._aniso_formula:
                vals = self._g.values()
            #else:
                #print "using aniso_formula"
            M = self._g
        else:
            vals = self._M.values()
            M = self._M

        if (2*k+s)%4 == 0:
            d = Integer(1)/Integer(2)*(m+n2) # |dimension of the Weil representation on even functions|
            self._d = d
            self._alpha4 = 1/Integer(2)*(vals[0]+v2) # the codimension of SkL in MkL
        else:
            d = Integer(1)/Integer(2)*(m-n2) # |dimension of the Weil representation on odd functions|
            self._d = d
            self._alpha4 = 1/Integer(2)*(vals[0]-v2) # the codimension of SkL in MkL
            
        prec = ceil(max(log(M.order(),2),52)+1)+17
        #print prec
        RR = RealField(prec)
        CC = ComplexField(prec)
        
        if debug > 0: print "d, m = {0}, {1}".format(d,m)
        eps = exp( 2 * CC.pi() * CC(0,1) * (s + 2*k) / Integer(4) )
        eps = round(real(eps))        
        if self._alpha3 is None or self._last_eps != eps:
            self._last_eps = eps
            if self._aniso_formula:
                self._alpha4 = 1
                self._alpha3 = -sum([BB(a)*mm for a,mm in self._v2.iteritems() if a != 0])
                #print self._alpha3
                self._alpha3 += Integer(d) - Integer(1) - self._g.beta_formula()
                #print self._alpha3, self._g.a5prime_formula()
                self._alpha3 = self._alpha3/RR(2)
            else:
                self._alpha3 = eps*sum([(1-a)*mm for a,mm in self._v2.iteritems() if a != 0])
                if debug>0: print "alpha3t = ", self._alpha3
                self._alpha3 += sum([(1-a)*mm for a,mm in vals.iteritems() if a != 0])
                #print self._alpha3
                self._alpha3 = self._alpha3 / Integer(2)
        alpha3 = self._alpha3
        alpha4 = self._alpha4
        if debug > 0: print alpha3, alpha4
        g1=M.char_invariant(1)
        g1=CC(g1[0]*g1[1])
        #print g1
        g2=M.char_invariant(2)
        g2=CC(g2[0]*g2[1])
        if debug > 0: print g2, g2.parent()
        g3=M.char_invariant(-3)
        g3=CC(g3[0]*g3[1])
        if debug > 0: print "eps = {0}".format(eps)
        if debug > 0: print "d/4 = {0}, m/4 = {1}, e^(2pi i (2k+s)/8) = {2}".format(RR(d) / RR(4), sqrt(RR(m)) / RR(4), CC(exp(2 * CC.pi() * CC(0,1) * (2 * k + s) / Integer(8))))
        if eps == 1:
            g2_2 = real(g2)
        else:
            g2_2 = imag(g2)*CC(0,1)
        alpha1 = RR(d) / RR(4) - sqrt(RR(m)) / RR(4)  * CC(exp(2 * CC.pi() * CC(0,1) * (2 * k + s) / Integer(8)) * g2_2)
        if debug > 0: print alpha1
        alpha2 = RR(d) / RR(3) + sqrt(RR(m)) / (3 * sqrt(RR(3))) * real(exp(CC(2 * CC.pi() * CC(0,1) * (4 * k + 3 * s - 10) / 24)) * (g1 + eps*g3))
        if debug > 0: print "alpha1 = {0}, alpha2 = {1}, alpha3 = {2}, g1 = {3}, g2 = {4}, g3 = {5}, d = {6}, k = {7}, s = {8}".format(alpha1, alpha2, alpha3, g1, g2, g3, d, k, s)
        dim = real(d + (d * k / Integer(12)) - alpha1 - alpha2 - alpha3)
        if debug > 0:
            print "dimension:", dim
        if abs(dim-round(dim)) > 1e-6:
            raise RuntimeError("Error ({0}) too large in dimension formula for {1} and k={2}".format(abs(dim-round(dim)), self._M if self._M is not None else self._g, k))
        dimr = dim
        dim = Integer(round(dim))
        if k >=2 and dim < 0:
            raise RuntimeError("Negative dimension (= {0}, {1})!".format(dim, dimr))
        return dim

def affine_global_height(affine_point, abs_val=lambda x: x.abs()):
    """The affine logarithmic height function (works only over the rationals?)."""
    # use .numerical_approx() if we want a float
    return log(affine_height(affine_point, abs_val))

 def _sage_(self):
     import sage.all as sage
     return sage.log(self.args[0]._sage_())

def log_var(X):
    return log(variance(X).sqrt(), 2)