Riele, H.J.J. (Herman) te

English

CWI's Institutional Repository

Numer. Math. 58, 661-667 (1991) Numerische Mathematik © Springer-Verlag 1991 A new lower bound for the de Bruijn-Newman constant H.J.J. te Riele Centre for Mathematics and Computer Science, P.O. Box 4079, 1009 AB Amsterdam, The Netherlands Received August 20, 1989 /July 30, 1990 Summary. Strong numerical evidence is presented for a new lower bound for the so-called de Bruijn-Newman constant. This constant is related to the Rie-mann hypothesis. The new bound, - 5, is suggested by high-precision floating-point computations, with a mantissa of 250 decimal digits, of i) the coefficients of a so-called Jensen polynomial of degree 406, ii) the so-called Sturm sequence corresponding to this polynomial which implies that it has two complex zeros, and iii) the two complex zeros of this polynomial. A proof of the new bound could be given if one would repeat the computations i) and iii) with a floating-point accuracy of at least 2600 decimal digits. Subject classification: AMS(MOS): 65E05, 30DXX, 11M26; CR: Gl.5. 1 Introduction Recently, Csordas et al. [2] have introduced the so-called de Bruijn-Newman constant A as follows. Let the function H,_(x), AE9f', be defined by "' (1) H,_(x):= J eAl2 cJi(t) cos(xt)dt, 0 where "' (2) tl>(t)= L (2n4 n2 e9 t-3n2 ne 5 ') exp(-n2 ne4 r). n=l The function <P satisfies the following properties: i) tl>(z) is analytic in the strip -n/8<'Xz<n/8; ii) '1>(t)='1>(-t), and '1>(t)>O(tE9£>); iii) for any t:>O, lim <t><n>(t) exp[(n-c:) e4 i]=0, for each n=O, 1, ... t~oo 662 H.J .J. te Riele The function H ,_ is an entire function of order one, and H ,,(x) is real for real x. From results of de Bruijn [1] it follows that if the Riemann hypothesis is true, then H;.(x) must possess only real zeros for any ).~0. C.M. Newman has shown [6] that there exists a real number A, -oo <A ~-t, such that i) H,.(x) has only real zeros when A.~A, and ii) H,.(x) has some non-real zeros when ).<A. This number A was baptized the de Bruijn-Newman constant in [2]. The truth of the Riemann hypothesis would imply that A~ 0, whereas Newman [6] conjec-tures that A~O. In [2] it was proved that A> -50. In this note we will describe high-precision computations which provide strong numerical evidence for the new bound A> - 5. Moreover, our computa-tions show that trying to prove this result, or improve upon it, would be a formidable task, unless the algorithm used could be improved substantially. The computations were carried out on the CDC Cyber 995 (about 2 h CPU time for testing), and on the CDC Cyber 205 (about 30 h CPU time for 'produc-tion') of SARA (The Academic Computer Centre Amsterdam). Brent's MP pack-age was an indispensable tool for the high-precision floating-point computations. Since this package has not been vectorized, we used the Cyber 205 just as an extremely fast scalar machine. This note will rely heavily on [2]. We assume the reader to have a copy of [2] at hand (slight change: in the present paper we write bm and /Jm instead of 6m and Pm)· 2 Algorithm and results If we expand the cosine in (1) in its Taylor series, we obtain (3) H ( )= ~ (- lrbm(A.) x 2 m ;. X L.., (2 ) I ' m=O m · where 00 bm(A.) = J t2 meM2 <P(t) d t, 0 m=O, 1, ... ; AElR. The n-th degree Jensen polynomial Gn(t; A.) associated with H ,_ is defined by (4) G ( . A.)•=~ (n) k!bk(A.) k n t, k:--0 k (2k)! t ' and it is shown in [2] that if there exists a positive integer m and a real number X such that Gm(t; X) possesses a non-real zero, then X<A. The problem is to find m, given X. In [2] a lower bound for A was constructively obtained as follows. For A.= -50 and n= 16 the moments bk(A.), k=O, ... , n were computed with a known precision, by means of Romberg quadrature. The approximate Jensen polyno-A new lower bound for the de Bruijn-Newman constant 663 mial which we obtain by using this numerical approximation of bk().) in (4) is denoted by gn(t; A). The so-called Jenkins algorithm was used to compute all the zeros of gn (including one complex zero and its complex conjugate). A theorem of Ostrowski was then invoked to find an upper bound for the distance of each of these zeros to the corresponding zeros of Gn. This error bound was small enough to guarantee that the complex zero of g 11,(t; -50) found by high-precision computation indeed is an approximation of a complex zero of the Jensen polynomial G16 (t; -50). The sensitivity of the zeros of polynomials to errors in their coefficients required that the computations were performed in very high precision. Csordas et al. [2] used 110 digits of precision for their proof that - 50 <A. As a partial check, we repeated their computations in double precision on a CDC Cyber 995 (which means an accuracy of about 28 decimal digits) and could reproduce the complex zero of g16 (t; -50) with an accuracy of about 20 decimal digits. This illustrates the large amount of extra work needed to provide a proof of the existence of complex zeros of the Jensen polynomial Gn(t; Jc). In order to improve the result of Csordas et al., we realized that the degree of the Jensen polynomial Gn(t; },) which possesses complex zeros, might grow very fast with A. Consequently, finding all the zeros of G"' 11= 1, 2, ... (in order to prove the existence of complex ones) might become very expensive. Therefore, we decided to use so-called Sturm sequences [ 4] to get an indication of the existence of any complex zeros of the given Jensen polynomial. The computation of a Sturm sequence is much simpler than computing all the zeros of a given polynomial. A Sturm sequence associated with a given polynomial p0 (x) of degree m is a sequence of polynomials p0 (x), p1 (x), ... of strictly decreasing degree which can be defined as follows: P1 (x)==p~(x), Pi- i(x) ==qi(x) Pi(x)-Pi+ 1 (x), i = 1, 2, .. ., where qi(x) is found by the Euclidean algorithm, such that the degree of P; + 1 (x) is less than the degree of P;(x). If p0 (x) has only simple zeros, P;(x) has degree m-i, and the Sturm sequence consists of m+ 1 polynomials p0 (x), .. ., Pm(x). Let v(a) be the number of sign changes in the sequence {p;(a)}i=o (where zero values are skipped). Then v(a)-v(b) is the number of real zeros of the polynomial p0 (x) on the interval [a, b]. Note that v(± co) can be determined by inspection of the signs of the highest degree coefficients of the polynomials in the Sturm sequence. Our algorithm now works as follows. Suppose we know Ao and m=m().0 ), which is the smallest value for which gm(t; ),0 ) has complex zeros (to start with, we take Ao= - 50 and m = 16 from [2]). Then for a new value of }, which is somewhat larger than Ao we compute /J;(J,), i=O, 1, .. ., and for each new /3; we compute the coefficients of the associated Jensen polynomial, and by means of the associated Sturm sequence, its number of real zeros on the interval (- oo, 0]. This is continued until we have found 11 for which gn(t; A) should have complex zeros. Then we try to compute a complex zero of this polynomial by means of the Newton process, where the starting value is chosen as follows. Let z = z(Ao) be the known complex zero of gm(t; ).0 ). We tabulate the values 664 H.J .J. te Riele Table I. Minimal degrees m(i.) of Jensen polynomials with complex zeros }. m().) Complex zeros of gm(t; ).) Accuracy used 9i ±'! -50 16 -220.9191117 7.092565255 28D -49 16 -217.9076244 5.773253615 28D -48 16 -214.9084360 4.111013736 28D -47 16 -211.9217860 1.006843660 28D -46 17 -202.21965S3 5.677704348 28D -45 17 -199.3211883 3.991036911 28D -44 17 - 196.4360833 0.462709708 28D -43 18 -187.4386728 4.830351149 28D -42 18 -184.6425759 2.749091911 28D -41 19 -176.228937S 4.96997S476 28D -40 19 -173.5216696 3.024436421 28D and 40D -30 27 -116.8258164 2.400S95686 28D and SOD -20 41 - 111.0654985 1.322239430 SOD -10 97 -45.53019819 0.156978360 75D -5 406 - 24.34071458 0.031926616 250D of g.(t; ).) and its derivative, for some real values of t around ~(z), and we look for a local positive minimum, or a local negative maximum. In our experience, such a minimum, or maximum, is easy to find if A. is not too far away from ),0• Then we take c+di as starting value for the Newton process where c is the value oft for which g.(t; A.) assumes its local minimum or maximum, and where d = J:(z). In this way we found complex zeros of g.(t; A.) for A.= - 50(1) -40, - 30, - 20, -10, - 5. Table 1 presents the values for which we have determined the polynomial g.(t; ),) of smallest degree with complex zeros by means of the asso-ciated Sturm sequence. This degree is denoted by m=m(A.). In all cases this Jensen polynomial has m(A.)- 2 real zeros. Table 1 also lists the complex zeros found, truncated to 10 decimal digits, and the accuracy used. For A. close to - 50, the degree of the Jensen polynomial with complex zeros does not increase too quickly with A.. However, from A.~ - 20 this pattern changes drastically, as Table 1 shows. As ). increases, the imaginary parts of the complex zeros found seem to tend to zero. Our computations do not provide a mathematical proof of the existence of complex zeros of G.(t; A.), although there is strong numerical evidence. A proof of the new bound - 5 <A along the lines of Csordas et al. would require an extension of the accuracy we used (250 decimal digits) to at least 2600 decimal digits. The following simple error analysis, carried out for the case A.= - 5, may 406 help to convince the reader. If we define g406 (t; - 5)=: I ak t 406 -\ then, since k=O the ak are numerical approximations of the coefficients of G406 , computed with an accuracy of 250 decimal digits, we may write G 406 (t) = g406 (t) + e h(t), A new lower bound for the de Bruijn-Newman constant 665 406 where c=l0- 250 and h(t)= L 6kakt406 -k with 16kl<l. If 0 is a simple zero k=O of g4o6' and 0, the corresponding zero of G406 , then we have, as an approxima-tion of the first order in r, ( cf. [7, formula ( 5.8.1)]): Hence, h(O) o, = o- 8 -,-(O). g406 IO -OI <c lh(O)I < g4o6(101) ' I g~o6 (OJI 8 I g~o6 ( 0)1 ' since the coefficients of g406 are positive. For 0= -24.3407 ... +0.03192 ... i we found g406 (i01)=4.837 ... x10 12 and lg~06 (0)1=7.824 ... x10- 17, so that 10, - 01<10- 221 and 8, is indeed a complex zero of G406 (t; - 5). 3 Computational details In this section we shall explain the main details of how we computed fim(),) and the Sturm sequences of gn(t; A.). We write bm(A.) as the sum a co (5) bm(A)= Jt 2 me)·'2 <P(t)dt+ J t 2 rneM'<P(t)dt 0 a (Csordas et al. used a= 1). An upper bound for the second integral of (5) is found as follows. The function t 2 me''2 has maximum value (for t=(-m/A.)t), so that I t2 me.<i'<P(t)dt<exp [m (-1-log _:1,i.)] ! <P(t) dt <;exp [ m (-1 +log _:1)J+5a-ne4a] (cf. [3], in. (3.7)). This bound is used, for given ), and m, to choose a such that the contribution of the second integral in (5) to the value of bm(A.) is negligi-ble, in view of the precision used. E.g. for A.= - 5, we chose a= 1.65. Form= 406, this yields an upper bound of 10-4oo on the value of the second integral in (5). The smallest bm(-5) we computed is b344(-5)=1.46822 ... x10- 73 . Since we worked with a precision of 250 decimal digits, it follows that the contribution of the second integral in (5) to { bm ( - 5) }!~60 is indeed negligible. Let <PN(t) denote the sum of the first N terms of (2); then we have (cf. [3], eq. (4.6)) 0 < <P(t)- <PN(t) < nN3 exp(5t-nN 2 e4') (0 ~ t <co). 666 H.J .J. te Riele Given t, the number N is chosen such that the right hand side is less than 10-A where A is the number of decimal digits of precision employed in the computations plus{ -log10 (first term of cJ>(t))} (since this first term determines the size of cP(t)). Since the N in the exp is dominating, it is suffi~ient (most of the time) to choose N to be the smallest integer larger than i/(e- 4 r(A log 10+5t)/rr). Using the same notation as in [2], we now have to compute the integral (6) a b;,;l(A.)= J t 2 me-lr2 <PN(t) dt 0 to sufficient accuracy. In [2] this was done by Romberg quadrature. However, by inspecting the Romberg table for b~>(A_), we noticed that when going from left to right, i.e., when comparing T;j with T;,j+ 1, the accuracy did decrease (rather than increase, as one would expect: cf., e.g., [7, p. 141]). Moreover, the most accurate results were found in the first column of the Romberg table (just the trapezoidal rule results for step a, a/2, a/4, ... ), and the convergence in this column was much faster than quadratic. An explanation is given by the fact that the integrand in bm(A.) is an even function, and under certain condi-tions given in Theorem 2.2 of [5] the convergence of the trapezoidal rule for such functions is exponential. The integrand bm(A.) happens to satisfy these condi-tions. Therefore, it is unnecessary to apply Romberg quadrature. We just applied the composite trapezoidal rule with step a, a/2, a/4, ... , until a sufficiently small correction was obtained. For the computation of /Jm( - 5) we never needed to work with a step less than a/1024. Before applying the trapezoidal rule, a table of values of e"'r2 <P(t) was precomputed for t = j a/1024, j = 0, ... , 1024, since (a selection of) these values are needed for each /Jm(A.). In the final steps, we always observed an approximate doubling of the number of correct digits upon halving the step. The Sturm sequence associated with the polynomial gm(t; A.) was computed as follows. Let p0 (x)==gm(x; A.) and let m-i ( ) '\' m-i-j Pi x == L, c;j x , i=O, ... , m. j=O The coefficients c0 j of p0 (x) are computed by means of the relation (which follows from (4)): ml Coj= j!(2 (m~j))l /Jm-j(A), j=O, ... , m. Since p1 (x)=p~(x) we have clj=(m-j) c0 j,j=O, ... , m-1. Let qi(x)==qi 0 x+q;1, i = 1, .. ., m - 1. Then, by applying the definition of a Sturm sequence given in § 2, we find qio and q; 1, i= 1, ... , m-l, from, qio ciO -ci-1.0 =0, qio Ci1 +qi1ciO-ci-1,1 =0, A new lower bound for the de Bruijn-Newman constant 667 and C;+ 1,i,j=O, ... , m-i-2 from and Ci+l,m-i-l from C;+ l,m-i-1 =qil Ci,m-i-Ci- J,m-i+ 1· The number of negative real zeros of gm(t; J.) is just the number of sign changes in the sequence {(-1r-ic; 0 }f= 0 minus the number of sign changes in the sequence {c;,m-Ji'=o· In some instances, when going from gm to gm+l• this difference dropped down from m sharply. It turned out that this was caused always by insufficient precision used in the computation of the Sturm sequence associated with gm+ 1. The normal pattern (i.e., finding v ( - ro )-v (0) = m + 1, or m-1) could be restored easily, by increasing the accuracy. We have used the Sturm sequence technique to get an easy numerical indica-tion for the existence of complex zeros of gm(t; J.). The actual computation of these complex zeros by means of the Newton process provides the numerical evidence for their existence. The error analysis given at the end of§ 2 is meant to convince the reader who would rightly notice that these numerical computa-tions do not provide a rigorous mathematical proof. Acknowledgements. We like to thank Dr. N.M. Temme for his help and, in particular, for pointing out [5] to us. We also acknowledge the provision of computer time on the Cyber 205 by the Dutch Working Group on the Use of Supercomputers. References !. de Bruijn, N.G.: The roots of trigonometric integrals. Duke Math. J. 17, 197-226 (1950) 2. Csordas, G., Norfolk, T.S., Varga, R.S.: A lower bound for the de Bruijn-Newman constant A*. Numer. Math. 52, 483-497 (1988) 3. Csordas, G., Norfolk, T.S., Varga, R.S.: The Riemann hypothesis and the Turan inequalities. Trans. Amer. Math. Soc. 296, 521-541 (1986) 4. Henrici, P.: Applied and computational complex analysis, vol. !. New York: Wiley 1977 5. Kress, R.: On the general Hermite cardinal interpolation. Math. Comput. 26, 925-933 (1972) 6. Newman, C.M.: Fourier transforms with only real zeros. Proc. Am. Math. Soc. 61, 245- 251 (1976) 7. Stoer, J., Bulirsch, R.: Introduction to numerical analysis, Third correct printing. Berlin, Heidelberg, New York: Springer 1983 

A new lower bound for the De Bruijn-Newman constant

Abstract

Similar works

Full text

Available Versions

CWI's Institutional Repository