A fast amortized algorithm for computing quadratic Dirichlet L-functions

A fast amortized algorithm for computing quadratic Dirichlet L-functions Abstract An algorithm to compute Dirichlet L-functions for many quadratic characters is derived. The algorithm is optimal (up to logarithmic factors) provided that the conductors of the characters under consideration span a dyadic window. 1. Introduction Let s=σ+it, where σ and t are real, q>1 a positive integer, χ a primitive quadratic Dirichlet character modulo q, and L(s,χ) the corresponding Dirichlet L-function. In this article, we derive an algorithm to compute L(s,χ) for all quadratic characters χ with conductor q∈[Q,Q+Δ). The algorithm is essentially optimal provided that Δ is of comparable size to Q, and so can be viewed as an ‘Odlyzko–Schönhage algorithm’ in the q-aspect (rather than the t-aspect). This represents progress in amortized algorithms for L-functions. To state the main theorem, let d(q) denote the number of divisors of q. Theorem 1.1 Given any t, there are constants Aj≔Aj(t), j∈{1,2,3,4}, such that for every ϵ>0and 1≤Δ<Q/2, the function L(1/2+it,χ)can be numerically evaluated to within ϵ for any quadratic characters χ with conductor q∈[Q,Q+Δ)using ≤A1d(q)log(Q/ϵ)elementary arithmetic operations on numbers of ≤A2log(Q/ϵ)bits provided that a precomputation requiring ≤A3Qlog4(Q/ϵ)operations and ≤A4Qlog4(Q/ϵ)bits of storage is first performed. The average of the divisor function d(q) over q∈[Q,Q+Δ), Δ≍Q, is ≪logQ (see e.g. [15, Section 12.1]). Hence, Theorem 1.1 says that the average complexity of computing L(1/2+it,χ) is poly-log in Q and 1/ϵ. This is a large improvement over the approximate functional equation which, in comparison, requires Q1/2+oϵ(1) time per evaluation. (There is currently no algorithm for computing L(s,χ) at a single χ that improves on the approximate functional equation in the q-aspect except for an algorithm due to the author in the case χ is to a powerful modulus, but this case excludes quadratic characters.) The assumption in the theorem that Δ<Q/2 is not essential, and Δ≪Q still works simply by dividing into further dyadic sections if needed. Moreover, the upper bound Q/2 can be replaced by cQ for any positive constant c<1, though the dependence of the constants Aj on c would then be of the form 1/log(1/c). The theorem is stated for s on the critical line because our main interest is in computing low zeros of quadratic Dirichlet L-functions. But the theorem applies just as well to other values of s; e.g. to s=1. Given the wide availability of fast computer memory nowadays, the algorithm underlying Theorem 1.1 is practical. In fact, using a clever implementation one can hope to save at least one power of the logarithm (and probably two powers in the case of storage space) on powers of the logarithm that appear in the theorem, and to remove the dependence on the requested accuracy ϵ in some cases. Of course, Theorem 1.1 is not useful if one is interested in evaluating L(1/2+it,χ) for fewer than Q different χ, since in this case the approximate functional equation fares better. However, in some applications, such as the locating of low zeros of all quadratic Dirichlet L-functions with conductor in a dyadic section, our theorem is suitable. And even if one is interested in smaller sections, the theorem may still be of utility. For the proof of the theorem yields that if Δ=Q1−δ, then the average time to compute all L(s,χ) with conductor in [Q,Q+Δ) is Qδ+oϵ(1) operations. This is an improvement over the approximate functional equation if δ<1/2. The strategy behind Theorem 1.1 is roughly this. First, we use the well-known formula involving quadratic Gauss sums to express each quadratic character χ as a Fourier series. This representation of χ implicitly involves a single application of quadratic reciprocity, which is important, as without it our algorithm achieves no savings. The resulting main sum for L(s,χ) is then an exponential sum of length Q1+oϵ(1) terms. The advantage of having replaced χ by a quadratic Gauss sum is that now each χ corresponds to a point on the unit circle where we want to evaluate the exponential sum. So we proceed similarly to the Odlyzko–Schönhage algorithm from here, applying a fast multipole method to numerically evaluate the exponential sum at all the required points on the unit circle in (Δ+Q)1+oϵ(1) time and requiring (Δ+Q)1+oϵ(1) bits of storage space. Actually, our algorithm is slightly more complicated than this because there is dependence on q through a weighting function Vz(w) in the main sum; see formula (2.7). But this dependence is relatively weak and is easily removed via a Taylor expansion. Another complication is that representing χ by a quadratic Gauss sum when q is composite and varying in a range necessitates using an inclusion-exclusion formula; see formula (4.5). Interestingly, Theorem 1.1 might generalize to certain families of elliptic curve L-functions, though there are important differences with the present case. This will be the subject of future work. Sections 2–5 present the essential ingredients of the algorithm in Theorem 1.1. Section 6 contains a proof of the theorem. Section 7 discusses improvements to the algorithm. And Section 8 contains proofs of auxiliary lemmas. Remark. An algorithm for computing the ‘complete’ family of all Dirichlet L-functions L(s,χ) to primitive χ modulo a given q was recently derived by Platt [11]. This algorithm uses the Fast Fourier Transform in the expected way. 2. The main sum Let χ be a non-principal quadratic Dirichlet character of modulus q>1. The Dirichlet L-function associated with χ is   L(s,χ)≔∑n=1∞χ(n)ns,σ>0. (2.1)The completed L-function is Λ(s,χ)≔qs/2γ(s)L(s,χ), where   γ(z)≔π−z/2Γ(z+a2),a≔1−χ(−1)2∈{0,1}. (2.2)If χ is primitive, then the functional equation takes on the simple form Λ(s,χ)=Λ(1−s,χ). Using the method of the smoothed functional equation (see e.g. [12]), one obtains the following formula to compute L(s,χ). Let Γ(z,w) denote the incomplete gamma function,   Γ(z,w)=wz∫1∞e−wyyz−1dy,Re(w)>0. (2.3)And let   Vz(w)≔Γ(z/2,w)Γ(z/2),ρ(s)≔q1/2−sγ(1−s)γ(s). (2.4)Then we have   L(s,χ)=∑n=1∞χ(n)nsVs+a(πn2/q)+ρ(s)∑n=1∞χ(n)n1−sV1−s+a(πn2/q). (2.5)The weighting function V ensures the absolute convergence of the series in (2.5). So the series provides an analytic continuation of L(s,χ) to the entire complex plane. Henceforth, we restrict to the critical line s=1/2+it. As a consequence of the functional equation, Λ(1/2+it,χ) is real-valued for real t. This desirable feature of Λ simplifies locating roots on the critical line by looking for sign changes of Λ(1/2+it). However, unlike L(1/2+it,χ) itself, the function Λ(1/2+it,χ) experiences exponential decay with t due to the built-in γ factor. This makes Λ unsuitable for numerical computations when t is large; e.g. very high accuracy is needed to confirm sign changes. In view of this, one is motivated to introduce the functions   θ(t,a)≔arg[qit/2γ(1/2+it)],Z(t,χ)≔eiθ(t,a)L(1/2+it,χ). (2.6)So Z(t,χ) is of the same modulus as L(1/2+it,χ) while also analytic and real-valued for real t. Furthermore, the formula (5) now reads   Z(t,χ)=2Reeiθ(t,a)∑n=1∞χ(n)n1/2+itV1/2+it+a(πn2/q). (2.7) As pointed out in [12], there is still a problem of catastrophic cancellation in (2.7) because the terms in the series get very large and oscillatory if t is large. To fix the problem, one can use the same modification in [12], after which the algorithm that we describe here works even when t is large. Our focus, though, will be on computing low zeros of Z(t,χ). So we might imagine without loss of generality that t∈[0,1] say (or t∈[0,c1] for any fixed c1>0). Then the problem of catastrophic cancellation no longer arises, and the formula (2.7) is good enough as is. One could also use formula (2.5) and work with Λ(1/2+it,χ) directly. In either case, the dependence of the constants Aj on t is roughly of the form eπt/4. There is a sharp drop-off in the size of the weighting function V when n is around the square-root of the analytic conductor q≔q(s)=q∣s∣/2π, where s=1/2+it. The drop-off is so quick that one can obtain a truncation error of at most ϵ1 simply on stopping the series in (2.7) after ≈2qlog(q/ϵ1) terms. Hence, we can choose a uniform truncation point N≔NQ,Δ for the series in (2.7) for all quadratic conductors q∈[Q,Q+Δ). This yields the formula   Z(t,χ)=2Re[eiθ(t,a)F(t,χ)+R1(t,N,q)], (2.8)where   F(t,χ)=∑n=1Nχ(n)n1/2+itV1/2+it+a(πn2/q). (2.9) We now enforce the bound Δ<Q, as in Theorem 1.1, as well as the bound N>2Q/π, which our choice of N will ultimately satisfy if Q is large enough. Then, given 0<ϵ1<1, we can ensure via Lemma 8.1 in Section 8 that   R1(t,N,Q,Δ)≔maxq∈[Q,Q+Δ)∣R1(t,N,q)∣<ϵ1, (2.10)subject to   N≥(2Q/π)log(Q/ϵ1). (2.11)Note that this choice of N satisfies our earlier working assumption N>2Q/π, provided that Q>104 say. As for the rotation phase θ(t,a) in (2.7), its numerical evaluation of this can be done efficiently using known numerical algorithms for the Gamma function (e.g. [9, 13]), or using the asymptotic expansion   θ(t,a)∼t2log(qt2πe)+π(2a−1)8+148t+⋯,t→∞, (2.12)for large t, which is derived via the Stirling formula. In summary, to compute Z(t,χ) for a given t and χ with conductor q∈[Q,Q+Δ), it is enough to focus on the sum F(t,χ) in (2.9). The treatment of F(t,χ) will be separated into several cases. To begin, recall that a real primitive character exists to a prime power modulus if and only if the modulus is an odd prime p in which case the character is the quadratic symbol   (np)=((−1)(p−1)/2pn), (2.13)or the modulus is 4 in which case the character is χ4(n)=(−1)(n−1)/21nodd, or the modulus is 8 in which case there are two characters, χ8(n)=(−1)(n2−1)/81nodd and χ4(n)χ8(n)=(−1)(n2+4n−5)/81nodd;see e.g. [3]. Hence, the only quadratic conductors q are products such moduli, subject to the moduli being relatively prime. If q is not divisible by 8, then there is one primitive quadratic character with modulus q, otherwise there are two primitive characters. The cases of even and odd q will be handled separately by our algorithm. Furthermore, in order to fix the form of the functional equation of the L-function, we deal with the possibilities a=0 or a=1 separately also. Since our algorithm will apply to each situation analogously, we henceforth fix a=0 and q is odd. This corresponds to positive odd fundamental discriminants, which we denote by Qodd+. In particular, there is now a unique χq associated with each q∈Qodd+, and we may unambiguously write F(t,q) instead of F(t,χq) for the main sum. 3. The Taylor expansion We use the Taylor expansion to remove the dependence on q in the weighting function V. To this end, define   Gz(w)≔∫1∞e−wyyz−1dy. (3.1)By definition, we have   V1/2+it(πn2/q)n1/2+it=C(t,q)G1/4+it/2(πn2/q),C(t,q)≔(π/q)1/4+it/2Γ(1/4+it/2). (3.2)We apply the Taylor expansion to Gz(w) in the w variable to obtain   G1/4+it/2(πn2/q)=∑r=0∞G1/4+it/2(r)(πn2/Q)r!(πn2(Q−q)Qq)r, (3.3)where Gz(r)(w) is the rth derivative of G with respect to w. Let   cr(t,n)≔G1/4+it/2(r)(πn2/Q)r!(πn2Q)r, (3.4)and let   Fr(t,q)≔∑n=1Ncr(t,n)χq(n). (3.5)Truncating the series expansion in (3.3) after R terms then yields   F(t,q)=C(t,q)∑r=0R−1(Q−qq)rFr(t,q)+R2(t,N,q,R). (3.6)Furthermore, given 0<ϵ2<1, we can ensure via Lemma 8.2 in Section 8 that   R2(t,N,Q,Δ,R)≔maxq∈[Q,Q+Δ)∣R2(t,N,q,R)∣<ϵ2, (3.7)subject to   R≥log(N/ϵ2)log(Q/Δ). (3.8)Note that R depends only logarithmically on Q and ϵ2, and so grows slowly with these parameters. The importance of the Taylor expansion (3.6) is that it allows one to avoid having to re-calculate the weighting function V for each new q. This has a large impact on practical running time, as was observed in [14]. Indeed, using the recursions described in [14], the computation of the derivatives of Gz(w) can be reduced quickly to that of computing Gz(w) itself. In turn, the function Gz(w) can be computed efficiently using one of the many algorithms for the incomplete Gamma function; see [12] for a survey. Remark. We still remain within the target complexity of Theorem 1.1 even if each evaluation of G1/4+it/2(πn2/Q) takes as much as Δ/N time and space. This is because there are only N values of G1/4+it/2(πn2/Q) that we need to compute. In summary, we may assume that the coefficients cr(t,n) have been precomputed. 4. Exponential sums Consider the quadratic exponential sum (Gauss sum)   gq(n)≔∑ℓ=02n−1eπiqℓ2/n. (4.1)By results in e.g. [2, 3, 5] we have the Fourier expansion   χq(n)=g(q)gq(2n)n,g(q)≔eπiq(q−2)/8−πi/822, (4.2)subject to q being odd and (n,q)=1. Moving the factor g(q), which is independent of n, to the outside, we see that   Fr(t,q)=g(q)∑n=1(n,q)=1Ncr(t,n)gq(2n)n. (4.3) We wish to remove the summation condition (n,q)=1 in (4.3) because it will prevent us from applying the multi-evaluation method in the next section. With this in mind, note that if q∈Qodd+, then q=p1⋯pω where ph are distinct odd primes. So letting   S˜r(t,q,a)≔∑1≤m≤N/acr(t,am)gq(2am)am, (4.4) (q)={primeP:P∣qandP≤N}, and applying inclusion-exclusion, Fr(t,q) can be expressed as   Fr(t,q)=g(q)∑h=0∣(q)∣(−1)h∑{P1,…,Ph}⊆(q)S˜r(t,q,P1⋯Ph). (4.5)If P1⋯Ph>N, then S˜r(t,q,P1⋯Ph)=0 because the sum is empty. In particular, we may restrict the inner sum in (4.5) to subsets {P1,…,Ph} such that the product is ≤N. Moreover, using the notation a=P1⋯Ph and b=q/a, we have gq(2am)=agb(2m). So, if we let   Sr(t,a,b)≔∑1≤m≤N/acr(t,am)gb(2m)m, (4.6)then   S˜r(t,q,P1…Ph)=aSr(t,a,b). (4.7) Put together, and in light of the preceding observations, it suffices to compute the sum Sr(t,a,b) for all triples (r,a,b) such that Q/a≤b<Q/a+Δ/a, 1≤a≤N, and 0≤r<R. After that, Z(t,χq) can be recovered for all odd quadratic conductors q∈[Q,Q+Δ) by back-substitution into formulas (4.7), (4.5), (3.6) and (2.8). 5. Multiple evaluations algorithm Consider an exponential sum of the form   Z(h)≔∑k=1Kake2πiαkh,0≤αk<1,∣ak∣≤1. (5.1)The Odlyzko–Schönhage algorithm [10] provides an efficient method to compute Z(h) for all integers h∈[−H/2,H/2). Briefly, the problem is transformed via the fast Fourier transform to that of evaluating a sum of rational functions of the form   ∑k=1Kβkz−e−2πiαk,βk≔ake−2πiαk(e2πiαkH−1), (5.2)at all the Hth roots of unity. After that, either the Odlyzko–Schönhage rational function method, which is based on Taylor expansions, or the Greengard–Rokhlin algorithm [7], which additionally uses certain recursions and translations, can be used to perform the rational function evaluation. The Greengard–Rokhlin algorithm seems to perform better in practice; see [6] for an example implementation in the context of the zeta function. (There is vast literature on improvements and generalizations of the Greengard–Rokhlin algorithm; see e.g. [1] for a survey.) Either of these methods yields an essentially optimal average time (up to logarithmic factors) for computing Z(h). In fact, using the notation so far and letting M=max{H,K}, one easily deduces the following from [10, Section 3 & Theorem 4.1]. Theorem 5.1 (Odlyzko–Schönhage) The function Z(h)can be computed for all integers h∈[−H/2,H/2)with error <ϵ3using ≪Mlog2(M/ϵ3)arithmetic operations on numbers of ≪log(M/ϵ3)bits and ≪Mlog2(M/ϵ3)bits of storage. One can also use a non-uniform fast Fourier transform method instead, which avoids the intermediate step involving the rational functions (5.2); see [4, 8]. This is likely be more efficient in practice. 6. Proof of Theorem 1.1 Proof We assume that Q>104, otherwise the computation of F(t,q) can be done directly for each quadratic conductor q∈[Q,Q+Δ). In what follows, we will use the estimate (valid for w>0 and Re(z)≤1)   ∣Gz(r)(w)∣≤∫1∞e−wyyrdy≤w−r−1∫0∞e−uurdu≤r!wr+1. (6.1)This estimate implies, in particular, that ∣cr(t,n)∣≤Q/π. This bound is implicitly used when we apply Theorem 5.1 next. (Theorem 5.1 assumes that the coefficients of the exponential sum Z(h) are bounded by ≪1, and more generally the conclusion of the theorem still holds if the coefficients are bounded by ≪M, say.) We use Theorem 5.1 to bound the complexity of computing, to within ϵ3∈(0,1), the sum   Sr(t,a,b)=∑1≤m≤N/a0≤ℓ<4mcr(t,am)eπibℓ2/2mm (6.2)for all b∈[Q/a,Q/a+Δ/a), a∈[1,N], and r∈[0,R). Given a, the length of the sum Sr(t,a,b) is about 2N(N+a)/a2 terms, and the number of points b where we want to evaluate it is ≤Δ/a+1 points. Therefore, by Theorem 5.1, there exists an absolute constant A5 such that the total cost of computing Sr(t,a,b) for a given a and r, and for all b∈[Q/a,Q/a+Δ/a), is   ≤A5(Δ/a+N2/a2)log2((Δ+N)/ϵ3) (6.3)arithmetic operations on numbers of log((Δ+N)/ϵ3) bits and requiring a similar amount of storage space. Summing over a and r, this cost is   ≤A6R(ΔlogN+N2)log2((Δ+N)/ϵ3) (6.4)arithmetic operations, where A6 is an absolute constant (e.g. we can take A6=2A5). We choose R to be the ceiling of the lower bound in (3.8). In addition, we note that Δ<2πN2/logN since Q>104 by assumption. So there is a constant A7 such that (6.4) is bounded by   ≤A7N2log(N/ϵ2)log(Q/Δ)log2(N/ϵ3). (6.5)Choosing N to be the ceiling of the lower bound in (2.11), and substituting into (6.5) now gives that this is   ≤A8Qlog(Q/ϵ1)log(Q/(ϵ1ϵ2))log(Q/Δ)log2(Q/(ϵ1ϵ3)), (6.6)where A8 is an absolute constant. We choose ϵ1=ϵ2=ϵ/8, and ϵ3=ϵ/(2RNQ). This choice of ϵ3 ensures that S˜r(t,q,a) can be recovered from Sr(t,a,b) via formula (4.7) to within ϵ/(2RQ). In turn, F(t,q) can be recovered from S˜r(t,q,a) via formulas (4.5) and (3.6) to within ϵ/4. (Here, we used that ∣g(q)C(t,q)∣≤0.34 for t∈[0,1], and we bounded the length of the sum in (4.5) by d(q)≤Q, which is very generous.) Overall, the recovery process thus described takes at most an additional ≤A9d(q)R operation for some absolute constant A9. Finally, the value of Z(t,χq) can be computed using formula (2.8), and the accumulated (roundoff and truncation) error in computing Z(t,χq) this way is <2ϵ1+2ϵ2+RNQϵ3≤ϵ. Note that that to execute formula (4.5), we needed information about the factorization of q. Acquiring this information causes little difficulty, as one can determine the factorization of each q∈[Q,Q+Δ) in poly-log time on average using an obvious generalization of the sieve of Eratosthenes.□ 7. Improvements The most significant improvement to Theorem 1.1 would be to reduce the window size Δ that is required for optimal performance, currently Δ≍Q. One might even hope for optimal performance over windows as small as Δ≍Q. Unfortunately, it is not clear how to do this, though a key step in that direction seems to involve repeated applications of the reciprocity law for the quadratic Gauss sum gb(m) in order reduce the length of the main sum. One can use various identities for gb(m) to obtain immediate savings. For example,   gb(2m)=4∑ℓ=0m−1eπibℓ2/2m+2eπim/2−1ifb≡1(mod4),e3πim/2−1ifb≡3(mod4). (7.1)Using this identity enables cutting the length of the sum Sr(t,a,b) from about 2(N/a)2 terms to about 12(N/a)2 terms, a speed-up by a factor of 4. There may be other identities and symmetries of gb(m) that can be similarly exploited. Another potentially interesting improvement is to obtain a ‘hybrid algorithm’ that works uniformly in the t and q aspects. One difficulty here is that the dependence on t in the main sum occurs implicitly, through the weighting function V. However, as suggested in [12], one can express V using a Riemann sum and then interchange the order of summation and integration. This could produce an exponential sum more amenable to fast evaluation via the algorithm underlying Theorem 5.1. 8. Lemmas Lemma 8.1 Assume that Q>104and Q>Δ. If N≥(2Q/π)log(Q/ϵ1), then ∣R1(t,N,Q,Δ)∣<ϵ1. Proof First, we recall the definition (2.4) of the weighting function Vz(w), the integral representation (2.3) of Γ(z,w), and the lower bound mint∈[0,1]∣Γ(1/4+it/2)∣>1. Using this, we deduce that the size of nth term in the tail of the series (2.7) is   ≤1n(πn2q)a/2+1/4∫1∞e−πn2y/qya/2−3/4dy<1n(qπn2)3/4−a/2e−πn2/q, (8.1)where we arrived at the last inequality on using that a≤1 together with the trivial bound ya/2−3/4≤1 for y≥1. Thus, summing the terms in the series (2.7) over n>N, bounding each term using (8.1), and estimating the sum by an integral (cf. [16, p. 320]) gives   ∣R1(t,N,Q,Δ)∣≤12(qπ)7/4e−πN2/qN2. (8.2)Taking N as in the lemma and using the bound q<2Q yields the desired bound.□ Lemma 8.2 Assume Q>104and Q>Δ. If R≥log(N/ϵ2)/log(Q/Δ), then ∣R2(t,N,Q,Δ,R)∣<ϵ2. Proof If we truncate the Taylor series (3.3) at r=R, then the remainder can be expressed using the well-known integral form for the Taylor remainder. By [14, Lemma 2.1], this remainder, call it υ(q,r), is of size   ∣υ(q,r)∣≤(q/Q−1)RR!Γ(R,πn2/q)≤(q/Q−1)RR. (8.3)Therefore, recalling (3.6) and using that ∣C(t,q)∣≤(π/q)1/4 for t∈[0,1], and ∑1≤n≤N1/n≤2N,we obtain   ∣R2(t,N,q,R)∣<2N(πq)1/4(q/Q−1)RR. (8.4)Since Q≤q<Q+Δ, this gives   ∣R2(t,N,Q,Δ,R)∣<2NR(πQ)1/4(ΔQ)R. (8.5)So, taking R as in the lemma more than suffices to give the desired bound.□ Funding Preparation of this material is partially supported by the National Science Foundation under agreements No. DMS-1406190. The author is pleased to thank ICERM where part of this work was conducted. Acknowledgements This paper was motivated by a question of Brian Conrey about computing with various families of L-functions during the ICERM special semester in Fall 2015. I would like to thank Leslie Greengard for pointing out references about the non-uniform fast Fourier transform. References 1 R. Beatson and L. Greengard, A short course on fast multipole methods, Wavelets, multilevel methods and elliptic PDEs (Leicester, 1996), Numer. Math. Sci. Comput. , Oxford University Press, New York, 1997, 1– 37. 2 B. C. Berndt, R. J. Evans and K. S. Williams, Gauss and Jacobi sums, Canadian Mathematical Society Series of Monographs and Advanced Texts , John Wiley & Sons, Inc., New York, 1998). A Wiley-Interscience Publication. 3 H. Davenport, Multiplicative number theory, Graduate Texts in Mathematics  vol. 74, third edn, Springer-Verlag, New York, 2000). Revised and with a preface by Hugh L. Montgomery. 4 A. Dutt and V. Rokhlin, Fast Fourier transforms for nonequispaced data, SIAM J. Sci. Comput.  14 ( 1993), 1368– 1393. Google Scholar CrossRef Search ADS   5 H. Fiedler, W. Jurkat and O. Körner, Asymptotic expansions of finite theta series, Acta Arith.  32 ( 1977), 129– 146. Google Scholar CrossRef Search ADS   6 X. Gourdon, The 1013 first zeros of the Riemann zeta function and zero computation at very large heights, Online manuscript ( 2004). 7 L. Greengard and V. Rokhlin, A fast algorithm for particle simulations, J. Comput. Phys.  73 ( 1987), 325– 348. Google Scholar CrossRef Search ADS   8 L. Greengard and J.-Y. Lee, Accelerating the nonuniform fast Fourier transform, SIAM Rev.  46 ( 2004), 443– 454. Google Scholar CrossRef Search ADS   9 C. Lanczos, A precision approximation of the gamma function, J. Soc. Indust. Appl. Math. Ser. B Numer. Anal.  1 ( 1964), 86– 96. Google Scholar CrossRef Search ADS   10 A. M. Odlyzko and A. Schönhage, Fast algorithms for multiple evaluations of the Riemann zeta function, Trans. Amer. Math. Soc.  309 ( 1988), 797– 809. Google Scholar CrossRef Search ADS   11 D. J. Platt, Numerical computations concerning the GRH, Math, Comp.  85 ( 2016), 3009– 3027. 12 M. Rubinstein, Computational methods and experiments in analytic number theory, Recent perspectives in random matrix theory and number theory, London Math. Soc. Lecture Note Ser.  vol. 322, Cambridge University Press, Cambridge, 2005, 425– 506. Google Scholar CrossRef Search ADS   13 J. L. Spouge, Computation of the gamma, digamma, and trigamma functions, SIAM J. Numer. Anal.  31 ( 1994), 931– 944. Google Scholar CrossRef Search ADS   14 J. Stopple, The quadratic character experiment, Experiment. Math.  18 ( 2009), 193– 200. Google Scholar CrossRef Search ADS   15 E. C. Titchmarsh, The theory of the Riemann zeta-function , second edn, The Clarendon Press Oxford University Press, New York, 1986), Edited and with a preface by D. R. Heath-Brown. 16 P. J. Weinberger, On small zeros of Dirichlet L-functions, Math. Comp.  29 ( 1975), 319– 328. Collection of articles dedicated to Derrick Henry Lehmer on the occasion of his seventieth birthday. © 2017. Published by Oxford University Press. All rights reserved. For permissions, please email: journals.permissions@oup.com http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png The Quarterly Journal of Mathematics Oxford University Press

A fast amortized algorithm for computing quadratic Dirichlet L-functions

Loading next page...
 
/lp/ou_press/a-fast-amortized-algorithm-for-computing-quadratic-dirichlet-l-DtFhyv10mo
Publisher
Oxford University Press
Copyright
© 2017. Published by Oxford University Press. All rights reserved. For permissions, please email: journals.permissions@oup.com
ISSN
0033-5606
eISSN
1464-3847
D.O.I.
10.1093/qmath/hax029
Publisher site
See Article on Publisher Site

Abstract

Abstract An algorithm to compute Dirichlet L-functions for many quadratic characters is derived. The algorithm is optimal (up to logarithmic factors) provided that the conductors of the characters under consideration span a dyadic window. 1. Introduction Let s=σ+it, where σ and t are real, q>1 a positive integer, χ a primitive quadratic Dirichlet character modulo q, and L(s,χ) the corresponding Dirichlet L-function. In this article, we derive an algorithm to compute L(s,χ) for all quadratic characters χ with conductor q∈[Q,Q+Δ). The algorithm is essentially optimal provided that Δ is of comparable size to Q, and so can be viewed as an ‘Odlyzko–Schönhage algorithm’ in the q-aspect (rather than the t-aspect). This represents progress in amortized algorithms for L-functions. To state the main theorem, let d(q) denote the number of divisors of q. Theorem 1.1 Given any t, there are constants Aj≔Aj(t), j∈{1,2,3,4}, such that for every ϵ>0and 1≤Δ<Q/2, the function L(1/2+it,χ)can be numerically evaluated to within ϵ for any quadratic characters χ with conductor q∈[Q,Q+Δ)using ≤A1d(q)log(Q/ϵ)elementary arithmetic operations on numbers of ≤A2log(Q/ϵ)bits provided that a precomputation requiring ≤A3Qlog4(Q/ϵ)operations and ≤A4Qlog4(Q/ϵ)bits of storage is first performed. The average of the divisor function d(q) over q∈[Q,Q+Δ), Δ≍Q, is ≪logQ (see e.g. [15, Section 12.1]). Hence, Theorem 1.1 says that the average complexity of computing L(1/2+it,χ) is poly-log in Q and 1/ϵ. This is a large improvement over the approximate functional equation which, in comparison, requires Q1/2+oϵ(1) time per evaluation. (There is currently no algorithm for computing L(s,χ) at a single χ that improves on the approximate functional equation in the q-aspect except for an algorithm due to the author in the case χ is to a powerful modulus, but this case excludes quadratic characters.) The assumption in the theorem that Δ<Q/2 is not essential, and Δ≪Q still works simply by dividing into further dyadic sections if needed. Moreover, the upper bound Q/2 can be replaced by cQ for any positive constant c<1, though the dependence of the constants Aj on c would then be of the form 1/log(1/c). The theorem is stated for s on the critical line because our main interest is in computing low zeros of quadratic Dirichlet L-functions. But the theorem applies just as well to other values of s; e.g. to s=1. Given the wide availability of fast computer memory nowadays, the algorithm underlying Theorem 1.1 is practical. In fact, using a clever implementation one can hope to save at least one power of the logarithm (and probably two powers in the case of storage space) on powers of the logarithm that appear in the theorem, and to remove the dependence on the requested accuracy ϵ in some cases. Of course, Theorem 1.1 is not useful if one is interested in evaluating L(1/2+it,χ) for fewer than Q different χ, since in this case the approximate functional equation fares better. However, in some applications, such as the locating of low zeros of all quadratic Dirichlet L-functions with conductor in a dyadic section, our theorem is suitable. And even if one is interested in smaller sections, the theorem may still be of utility. For the proof of the theorem yields that if Δ=Q1−δ, then the average time to compute all L(s,χ) with conductor in [Q,Q+Δ) is Qδ+oϵ(1) operations. This is an improvement over the approximate functional equation if δ<1/2. The strategy behind Theorem 1.1 is roughly this. First, we use the well-known formula involving quadratic Gauss sums to express each quadratic character χ as a Fourier series. This representation of χ implicitly involves a single application of quadratic reciprocity, which is important, as without it our algorithm achieves no savings. The resulting main sum for L(s,χ) is then an exponential sum of length Q1+oϵ(1) terms. The advantage of having replaced χ by a quadratic Gauss sum is that now each χ corresponds to a point on the unit circle where we want to evaluate the exponential sum. So we proceed similarly to the Odlyzko–Schönhage algorithm from here, applying a fast multipole method to numerically evaluate the exponential sum at all the required points on the unit circle in (Δ+Q)1+oϵ(1) time and requiring (Δ+Q)1+oϵ(1) bits of storage space. Actually, our algorithm is slightly more complicated than this because there is dependence on q through a weighting function Vz(w) in the main sum; see formula (2.7). But this dependence is relatively weak and is easily removed via a Taylor expansion. Another complication is that representing χ by a quadratic Gauss sum when q is composite and varying in a range necessitates using an inclusion-exclusion formula; see formula (4.5). Interestingly, Theorem 1.1 might generalize to certain families of elliptic curve L-functions, though there are important differences with the present case. This will be the subject of future work. Sections 2–5 present the essential ingredients of the algorithm in Theorem 1.1. Section 6 contains a proof of the theorem. Section 7 discusses improvements to the algorithm. And Section 8 contains proofs of auxiliary lemmas. Remark. An algorithm for computing the ‘complete’ family of all Dirichlet L-functions L(s,χ) to primitive χ modulo a given q was recently derived by Platt [11]. This algorithm uses the Fast Fourier Transform in the expected way. 2. The main sum Let χ be a non-principal quadratic Dirichlet character of modulus q>1. The Dirichlet L-function associated with χ is   L(s,χ)≔∑n=1∞χ(n)ns,σ>0. (2.1)The completed L-function is Λ(s,χ)≔qs/2γ(s)L(s,χ), where   γ(z)≔π−z/2Γ(z+a2),a≔1−χ(−1)2∈{0,1}. (2.2)If χ is primitive, then the functional equation takes on the simple form Λ(s,χ)=Λ(1−s,χ). Using the method of the smoothed functional equation (see e.g. [12]), one obtains the following formula to compute L(s,χ). Let Γ(z,w) denote the incomplete gamma function,   Γ(z,w)=wz∫1∞e−wyyz−1dy,Re(w)>0. (2.3)And let   Vz(w)≔Γ(z/2,w)Γ(z/2),ρ(s)≔q1/2−sγ(1−s)γ(s). (2.4)Then we have   L(s,χ)=∑n=1∞χ(n)nsVs+a(πn2/q)+ρ(s)∑n=1∞χ(n)n1−sV1−s+a(πn2/q). (2.5)The weighting function V ensures the absolute convergence of the series in (2.5). So the series provides an analytic continuation of L(s,χ) to the entire complex plane. Henceforth, we restrict to the critical line s=1/2+it. As a consequence of the functional equation, Λ(1/2+it,χ) is real-valued for real t. This desirable feature of Λ simplifies locating roots on the critical line by looking for sign changes of Λ(1/2+it). However, unlike L(1/2+it,χ) itself, the function Λ(1/2+it,χ) experiences exponential decay with t due to the built-in γ factor. This makes Λ unsuitable for numerical computations when t is large; e.g. very high accuracy is needed to confirm sign changes. In view of this, one is motivated to introduce the functions   θ(t,a)≔arg[qit/2γ(1/2+it)],Z(t,χ)≔eiθ(t,a)L(1/2+it,χ). (2.6)So Z(t,χ) is of the same modulus as L(1/2+it,χ) while also analytic and real-valued for real t. Furthermore, the formula (5) now reads   Z(t,χ)=2Reeiθ(t,a)∑n=1∞χ(n)n1/2+itV1/2+it+a(πn2/q). (2.7) As pointed out in [12], there is still a problem of catastrophic cancellation in (2.7) because the terms in the series get very large and oscillatory if t is large. To fix the problem, one can use the same modification in [12], after which the algorithm that we describe here works even when t is large. Our focus, though, will be on computing low zeros of Z(t,χ). So we might imagine without loss of generality that t∈[0,1] say (or t∈[0,c1] for any fixed c1>0). Then the problem of catastrophic cancellation no longer arises, and the formula (2.7) is good enough as is. One could also use formula (2.5) and work with Λ(1/2+it,χ) directly. In either case, the dependence of the constants Aj on t is roughly of the form eπt/4. There is a sharp drop-off in the size of the weighting function V when n is around the square-root of the analytic conductor q≔q(s)=q∣s∣/2π, where s=1/2+it. The drop-off is so quick that one can obtain a truncation error of at most ϵ1 simply on stopping the series in (2.7) after ≈2qlog(q/ϵ1) terms. Hence, we can choose a uniform truncation point N≔NQ,Δ for the series in (2.7) for all quadratic conductors q∈[Q,Q+Δ). This yields the formula   Z(t,χ)=2Re[eiθ(t,a)F(t,χ)+R1(t,N,q)], (2.8)where   F(t,χ)=∑n=1Nχ(n)n1/2+itV1/2+it+a(πn2/q). (2.9) We now enforce the bound Δ<Q, as in Theorem 1.1, as well as the bound N>2Q/π, which our choice of N will ultimately satisfy if Q is large enough. Then, given 0<ϵ1<1, we can ensure via Lemma 8.1 in Section 8 that   R1(t,N,Q,Δ)≔maxq∈[Q,Q+Δ)∣R1(t,N,q)∣<ϵ1, (2.10)subject to   N≥(2Q/π)log(Q/ϵ1). (2.11)Note that this choice of N satisfies our earlier working assumption N>2Q/π, provided that Q>104 say. As for the rotation phase θ(t,a) in (2.7), its numerical evaluation of this can be done efficiently using known numerical algorithms for the Gamma function (e.g. [9, 13]), or using the asymptotic expansion   θ(t,a)∼t2log(qt2πe)+π(2a−1)8+148t+⋯,t→∞, (2.12)for large t, which is derived via the Stirling formula. In summary, to compute Z(t,χ) for a given t and χ with conductor q∈[Q,Q+Δ), it is enough to focus on the sum F(t,χ) in (2.9). The treatment of F(t,χ) will be separated into several cases. To begin, recall that a real primitive character exists to a prime power modulus if and only if the modulus is an odd prime p in which case the character is the quadratic symbol   (np)=((−1)(p−1)/2pn), (2.13)or the modulus is 4 in which case the character is χ4(n)=(−1)(n−1)/21nodd, or the modulus is 8 in which case there are two characters, χ8(n)=(−1)(n2−1)/81nodd and χ4(n)χ8(n)=(−1)(n2+4n−5)/81nodd;see e.g. [3]. Hence, the only quadratic conductors q are products such moduli, subject to the moduli being relatively prime. If q is not divisible by 8, then there is one primitive quadratic character with modulus q, otherwise there are two primitive characters. The cases of even and odd q will be handled separately by our algorithm. Furthermore, in order to fix the form of the functional equation of the L-function, we deal with the possibilities a=0 or a=1 separately also. Since our algorithm will apply to each situation analogously, we henceforth fix a=0 and q is odd. This corresponds to positive odd fundamental discriminants, which we denote by Qodd+. In particular, there is now a unique χq associated with each q∈Qodd+, and we may unambiguously write F(t,q) instead of F(t,χq) for the main sum. 3. The Taylor expansion We use the Taylor expansion to remove the dependence on q in the weighting function V. To this end, define   Gz(w)≔∫1∞e−wyyz−1dy. (3.1)By definition, we have   V1/2+it(πn2/q)n1/2+it=C(t,q)G1/4+it/2(πn2/q),C(t,q)≔(π/q)1/4+it/2Γ(1/4+it/2). (3.2)We apply the Taylor expansion to Gz(w) in the w variable to obtain   G1/4+it/2(πn2/q)=∑r=0∞G1/4+it/2(r)(πn2/Q)r!(πn2(Q−q)Qq)r, (3.3)where Gz(r)(w) is the rth derivative of G with respect to w. Let   cr(t,n)≔G1/4+it/2(r)(πn2/Q)r!(πn2Q)r, (3.4)and let   Fr(t,q)≔∑n=1Ncr(t,n)χq(n). (3.5)Truncating the series expansion in (3.3) after R terms then yields   F(t,q)=C(t,q)∑r=0R−1(Q−qq)rFr(t,q)+R2(t,N,q,R). (3.6)Furthermore, given 0<ϵ2<1, we can ensure via Lemma 8.2 in Section 8 that   R2(t,N,Q,Δ,R)≔maxq∈[Q,Q+Δ)∣R2(t,N,q,R)∣<ϵ2, (3.7)subject to   R≥log(N/ϵ2)log(Q/Δ). (3.8)Note that R depends only logarithmically on Q and ϵ2, and so grows slowly with these parameters. The importance of the Taylor expansion (3.6) is that it allows one to avoid having to re-calculate the weighting function V for each new q. This has a large impact on practical running time, as was observed in [14]. Indeed, using the recursions described in [14], the computation of the derivatives of Gz(w) can be reduced quickly to that of computing Gz(w) itself. In turn, the function Gz(w) can be computed efficiently using one of the many algorithms for the incomplete Gamma function; see [12] for a survey. Remark. We still remain within the target complexity of Theorem 1.1 even if each evaluation of G1/4+it/2(πn2/Q) takes as much as Δ/N time and space. This is because there are only N values of G1/4+it/2(πn2/Q) that we need to compute. In summary, we may assume that the coefficients cr(t,n) have been precomputed. 4. Exponential sums Consider the quadratic exponential sum (Gauss sum)   gq(n)≔∑ℓ=02n−1eπiqℓ2/n. (4.1)By results in e.g. [2, 3, 5] we have the Fourier expansion   χq(n)=g(q)gq(2n)n,g(q)≔eπiq(q−2)/8−πi/822, (4.2)subject to q being odd and (n,q)=1. Moving the factor g(q), which is independent of n, to the outside, we see that   Fr(t,q)=g(q)∑n=1(n,q)=1Ncr(t,n)gq(2n)n. (4.3) We wish to remove the summation condition (n,q)=1 in (4.3) because it will prevent us from applying the multi-evaluation method in the next section. With this in mind, note that if q∈Qodd+, then q=p1⋯pω where ph are distinct odd primes. So letting   S˜r(t,q,a)≔∑1≤m≤N/acr(t,am)gq(2am)am, (4.4) (q)={primeP:P∣qandP≤N}, and applying inclusion-exclusion, Fr(t,q) can be expressed as   Fr(t,q)=g(q)∑h=0∣(q)∣(−1)h∑{P1,…,Ph}⊆(q)S˜r(t,q,P1⋯Ph). (4.5)If P1⋯Ph>N, then S˜r(t,q,P1⋯Ph)=0 because the sum is empty. In particular, we may restrict the inner sum in (4.5) to subsets {P1,…,Ph} such that the product is ≤N. Moreover, using the notation a=P1⋯Ph and b=q/a, we have gq(2am)=agb(2m). So, if we let   Sr(t,a,b)≔∑1≤m≤N/acr(t,am)gb(2m)m, (4.6)then   S˜r(t,q,P1…Ph)=aSr(t,a,b). (4.7) Put together, and in light of the preceding observations, it suffices to compute the sum Sr(t,a,b) for all triples (r,a,b) such that Q/a≤b<Q/a+Δ/a, 1≤a≤N, and 0≤r<R. After that, Z(t,χq) can be recovered for all odd quadratic conductors q∈[Q,Q+Δ) by back-substitution into formulas (4.7), (4.5), (3.6) and (2.8). 5. Multiple evaluations algorithm Consider an exponential sum of the form   Z(h)≔∑k=1Kake2πiαkh,0≤αk<1,∣ak∣≤1. (5.1)The Odlyzko–Schönhage algorithm [10] provides an efficient method to compute Z(h) for all integers h∈[−H/2,H/2). Briefly, the problem is transformed via the fast Fourier transform to that of evaluating a sum of rational functions of the form   ∑k=1Kβkz−e−2πiαk,βk≔ake−2πiαk(e2πiαkH−1), (5.2)at all the Hth roots of unity. After that, either the Odlyzko–Schönhage rational function method, which is based on Taylor expansions, or the Greengard–Rokhlin algorithm [7], which additionally uses certain recursions and translations, can be used to perform the rational function evaluation. The Greengard–Rokhlin algorithm seems to perform better in practice; see [6] for an example implementation in the context of the zeta function. (There is vast literature on improvements and generalizations of the Greengard–Rokhlin algorithm; see e.g. [1] for a survey.) Either of these methods yields an essentially optimal average time (up to logarithmic factors) for computing Z(h). In fact, using the notation so far and letting M=max{H,K}, one easily deduces the following from [10, Section 3 & Theorem 4.1]. Theorem 5.1 (Odlyzko–Schönhage) The function Z(h)can be computed for all integers h∈[−H/2,H/2)with error <ϵ3using ≪Mlog2(M/ϵ3)arithmetic operations on numbers of ≪log(M/ϵ3)bits and ≪Mlog2(M/ϵ3)bits of storage. One can also use a non-uniform fast Fourier transform method instead, which avoids the intermediate step involving the rational functions (5.2); see [4, 8]. This is likely be more efficient in practice. 6. Proof of Theorem 1.1 Proof We assume that Q>104, otherwise the computation of F(t,q) can be done directly for each quadratic conductor q∈[Q,Q+Δ). In what follows, we will use the estimate (valid for w>0 and Re(z)≤1)   ∣Gz(r)(w)∣≤∫1∞e−wyyrdy≤w−r−1∫0∞e−uurdu≤r!wr+1. (6.1)This estimate implies, in particular, that ∣cr(t,n)∣≤Q/π. This bound is implicitly used when we apply Theorem 5.1 next. (Theorem 5.1 assumes that the coefficients of the exponential sum Z(h) are bounded by ≪1, and more generally the conclusion of the theorem still holds if the coefficients are bounded by ≪M, say.) We use Theorem 5.1 to bound the complexity of computing, to within ϵ3∈(0,1), the sum   Sr(t,a,b)=∑1≤m≤N/a0≤ℓ<4mcr(t,am)eπibℓ2/2mm (6.2)for all b∈[Q/a,Q/a+Δ/a), a∈[1,N], and r∈[0,R). Given a, the length of the sum Sr(t,a,b) is about 2N(N+a)/a2 terms, and the number of points b where we want to evaluate it is ≤Δ/a+1 points. Therefore, by Theorem 5.1, there exists an absolute constant A5 such that the total cost of computing Sr(t,a,b) for a given a and r, and for all b∈[Q/a,Q/a+Δ/a), is   ≤A5(Δ/a+N2/a2)log2((Δ+N)/ϵ3) (6.3)arithmetic operations on numbers of log((Δ+N)/ϵ3) bits and requiring a similar amount of storage space. Summing over a and r, this cost is   ≤A6R(ΔlogN+N2)log2((Δ+N)/ϵ3) (6.4)arithmetic operations, where A6 is an absolute constant (e.g. we can take A6=2A5). We choose R to be the ceiling of the lower bound in (3.8). In addition, we note that Δ<2πN2/logN since Q>104 by assumption. So there is a constant A7 such that (6.4) is bounded by   ≤A7N2log(N/ϵ2)log(Q/Δ)log2(N/ϵ3). (6.5)Choosing N to be the ceiling of the lower bound in (2.11), and substituting into (6.5) now gives that this is   ≤A8Qlog(Q/ϵ1)log(Q/(ϵ1ϵ2))log(Q/Δ)log2(Q/(ϵ1ϵ3)), (6.6)where A8 is an absolute constant. We choose ϵ1=ϵ2=ϵ/8, and ϵ3=ϵ/(2RNQ). This choice of ϵ3 ensures that S˜r(t,q,a) can be recovered from Sr(t,a,b) via formula (4.7) to within ϵ/(2RQ). In turn, F(t,q) can be recovered from S˜r(t,q,a) via formulas (4.5) and (3.6) to within ϵ/4. (Here, we used that ∣g(q)C(t,q)∣≤0.34 for t∈[0,1], and we bounded the length of the sum in (4.5) by d(q)≤Q, which is very generous.) Overall, the recovery process thus described takes at most an additional ≤A9d(q)R operation for some absolute constant A9. Finally, the value of Z(t,χq) can be computed using formula (2.8), and the accumulated (roundoff and truncation) error in computing Z(t,χq) this way is <2ϵ1+2ϵ2+RNQϵ3≤ϵ. Note that that to execute formula (4.5), we needed information about the factorization of q. Acquiring this information causes little difficulty, as one can determine the factorization of each q∈[Q,Q+Δ) in poly-log time on average using an obvious generalization of the sieve of Eratosthenes.□ 7. Improvements The most significant improvement to Theorem 1.1 would be to reduce the window size Δ that is required for optimal performance, currently Δ≍Q. One might even hope for optimal performance over windows as small as Δ≍Q. Unfortunately, it is not clear how to do this, though a key step in that direction seems to involve repeated applications of the reciprocity law for the quadratic Gauss sum gb(m) in order reduce the length of the main sum. One can use various identities for gb(m) to obtain immediate savings. For example,   gb(2m)=4∑ℓ=0m−1eπibℓ2/2m+2eπim/2−1ifb≡1(mod4),e3πim/2−1ifb≡3(mod4). (7.1)Using this identity enables cutting the length of the sum Sr(t,a,b) from about 2(N/a)2 terms to about 12(N/a)2 terms, a speed-up by a factor of 4. There may be other identities and symmetries of gb(m) that can be similarly exploited. Another potentially interesting improvement is to obtain a ‘hybrid algorithm’ that works uniformly in the t and q aspects. One difficulty here is that the dependence on t in the main sum occurs implicitly, through the weighting function V. However, as suggested in [12], one can express V using a Riemann sum and then interchange the order of summation and integration. This could produce an exponential sum more amenable to fast evaluation via the algorithm underlying Theorem 5.1. 8. Lemmas Lemma 8.1 Assume that Q>104and Q>Δ. If N≥(2Q/π)log(Q/ϵ1), then ∣R1(t,N,Q,Δ)∣<ϵ1. Proof First, we recall the definition (2.4) of the weighting function Vz(w), the integral representation (2.3) of Γ(z,w), and the lower bound mint∈[0,1]∣Γ(1/4+it/2)∣>1. Using this, we deduce that the size of nth term in the tail of the series (2.7) is   ≤1n(πn2q)a/2+1/4∫1∞e−πn2y/qya/2−3/4dy<1n(qπn2)3/4−a/2e−πn2/q, (8.1)where we arrived at the last inequality on using that a≤1 together with the trivial bound ya/2−3/4≤1 for y≥1. Thus, summing the terms in the series (2.7) over n>N, bounding each term using (8.1), and estimating the sum by an integral (cf. [16, p. 320]) gives   ∣R1(t,N,Q,Δ)∣≤12(qπ)7/4e−πN2/qN2. (8.2)Taking N as in the lemma and using the bound q<2Q yields the desired bound.□ Lemma 8.2 Assume Q>104and Q>Δ. If R≥log(N/ϵ2)/log(Q/Δ), then ∣R2(t,N,Q,Δ,R)∣<ϵ2. Proof If we truncate the Taylor series (3.3) at r=R, then the remainder can be expressed using the well-known integral form for the Taylor remainder. By [14, Lemma 2.1], this remainder, call it υ(q,r), is of size   ∣υ(q,r)∣≤(q/Q−1)RR!Γ(R,πn2/q)≤(q/Q−1)RR. (8.3)Therefore, recalling (3.6) and using that ∣C(t,q)∣≤(π/q)1/4 for t∈[0,1], and ∑1≤n≤N1/n≤2N,we obtain   ∣R2(t,N,q,R)∣<2N(πq)1/4(q/Q−1)RR. (8.4)Since Q≤q<Q+Δ, this gives   ∣R2(t,N,Q,Δ,R)∣<2NR(πQ)1/4(ΔQ)R. (8.5)So, taking R as in the lemma more than suffices to give the desired bound.□ Funding Preparation of this material is partially supported by the National Science Foundation under agreements No. DMS-1406190. The author is pleased to thank ICERM where part of this work was conducted. Acknowledgements This paper was motivated by a question of Brian Conrey about computing with various families of L-functions during the ICERM special semester in Fall 2015. I would like to thank Leslie Greengard for pointing out references about the non-uniform fast Fourier transform. References 1 R. Beatson and L. Greengard, A short course on fast multipole methods, Wavelets, multilevel methods and elliptic PDEs (Leicester, 1996), Numer. Math. Sci. Comput. , Oxford University Press, New York, 1997, 1– 37. 2 B. C. Berndt, R. J. Evans and K. S. Williams, Gauss and Jacobi sums, Canadian Mathematical Society Series of Monographs and Advanced Texts , John Wiley & Sons, Inc., New York, 1998). A Wiley-Interscience Publication. 3 H. Davenport, Multiplicative number theory, Graduate Texts in Mathematics  vol. 74, third edn, Springer-Verlag, New York, 2000). Revised and with a preface by Hugh L. Montgomery. 4 A. Dutt and V. Rokhlin, Fast Fourier transforms for nonequispaced data, SIAM J. Sci. Comput.  14 ( 1993), 1368– 1393. Google Scholar CrossRef Search ADS   5 H. Fiedler, W. Jurkat and O. Körner, Asymptotic expansions of finite theta series, Acta Arith.  32 ( 1977), 129– 146. Google Scholar CrossRef Search ADS   6 X. Gourdon, The 1013 first zeros of the Riemann zeta function and zero computation at very large heights, Online manuscript ( 2004). 7 L. Greengard and V. Rokhlin, A fast algorithm for particle simulations, J. Comput. Phys.  73 ( 1987), 325– 348. Google Scholar CrossRef Search ADS   8 L. Greengard and J.-Y. Lee, Accelerating the nonuniform fast Fourier transform, SIAM Rev.  46 ( 2004), 443– 454. Google Scholar CrossRef Search ADS   9 C. Lanczos, A precision approximation of the gamma function, J. Soc. Indust. Appl. Math. Ser. B Numer. Anal.  1 ( 1964), 86– 96. Google Scholar CrossRef Search ADS   10 A. M. Odlyzko and A. Schönhage, Fast algorithms for multiple evaluations of the Riemann zeta function, Trans. Amer. Math. Soc.  309 ( 1988), 797– 809. Google Scholar CrossRef Search ADS   11 D. J. Platt, Numerical computations concerning the GRH, Math, Comp.  85 ( 2016), 3009– 3027. 12 M. Rubinstein, Computational methods and experiments in analytic number theory, Recent perspectives in random matrix theory and number theory, London Math. Soc. Lecture Note Ser.  vol. 322, Cambridge University Press, Cambridge, 2005, 425– 506. Google Scholar CrossRef Search ADS   13 J. L. Spouge, Computation of the gamma, digamma, and trigamma functions, SIAM J. Numer. Anal.  31 ( 1994), 931– 944. Google Scholar CrossRef Search ADS   14 J. Stopple, The quadratic character experiment, Experiment. Math.  18 ( 2009), 193– 200. Google Scholar CrossRef Search ADS   15 E. C. Titchmarsh, The theory of the Riemann zeta-function , second edn, The Clarendon Press Oxford University Press, New York, 1986), Edited and with a preface by D. R. Heath-Brown. 16 P. J. Weinberger, On small zeros of Dirichlet L-functions, Math. Comp.  29 ( 1975), 319– 328. Collection of articles dedicated to Derrick Henry Lehmer on the occasion of his seventieth birthday. © 2017. Published by Oxford University Press. All rights reserved. For permissions, please email: journals.permissions@oup.com

Journal

The Quarterly Journal of MathematicsOxford University Press

Published: Mar 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Unlimited reading

Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere.

Stay up to date

Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates.

Organize your research

It’s easy to organize your research with our built-in tools.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

Monthly Plan

  • Read unlimited articles
  • Personalized recommendations
  • No expiration
  • Print 20 pages per month
  • 20% off on PDF purchases
  • Organize your research
  • Get updates on your journals and topic searches

$49/month

Start Free Trial

14-day Free Trial

Best Deal — 39% off

Annual Plan

  • All the features of the Professional Plan, but for 39% off!
  • Billed annually
  • No expiration
  • For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588

$360/year

billed annually
Start Free Trial

14-day Free Trial