Learning with GP: Griffin,Ono,Rolen,Zaghier asymptotic formula for xi(s)

[Jacques Gélinas <jacquesg00@hotmail.com>]

The asymptotic formula Fsd for the Taylor coefficients of the Riemann xi function
at s=1/2 obtained recently in [1] gives quite a few correct leading digits !

1. The Riemann integral representation XiR for \Xi(it) is simpler than the usual
cosine transform reformulation of Jensen; \gamma(n) in [1, eq. [1], Table 2]
corresponds to our 8*n!*XikR(2*n). The first five Xi(x) MacLaurin coefficients
of x^[0,2,4,6,8] are computed via Double Exponential numerical integration and
compared to GP series values for equality of floating-point significant digits.

\\-------------------------------------------- Start cut
near(x,y,r=7/8) = exponent(normlp(x-y)/(normlp(x)+!x))/exponent(0.)>=r;
maclv(f,k) = Vec(f( x+O(x^(k+1)) ));

Xi(x) = my(s=1/2+x); Pi^(-s/2)*gamma(1+s/2)*(s-1)*zeta(s);

psik(x,k=0) = sumpos(n=1, exp(-n^2*Pi*x)*(-n^2*Pi)^k );
near(0, 1/2 + psik(1,0) + 4*psik(1,1) )    \\ from Jacobi theta identity

Itab = intnuminit(1,[oo,1]);               \\ reuse DE nodes and weights
XiR(t) = 1/2 - (t^2 + 1/4)*\
           intnum(x=1,[oo,1], psik(x) * x^(-3/4) * cos(t/2*log(x)), Itab);
F(n) = intnum(x=1,[oo,1], psik(x) * x^(-3/4) * log(x)^n, Itab );
XikR(k) = if(k%2, 0, if(!k, 1/2, F(k-2)/(k-2)!/2^(k-2)) - F(k)/k!/2^k/4);
near( maclv(Xi,8), apply(XikR,[0..8]) )
\\-------------------------------------------- End paste

2. Only the first term of psik(x) is kept to estimate F(n).
The modified integrand f(t) = exp(-Pi*t)*t^(-3/4)*log(t)^n
has its maximum at t=a, solution of n = (Pi*a + 3/4)*log(a).
Following the saddle point method, G and exp(G) are expanded
where G=log[f((1+y)*a)/f(a)] and z=1/log(a) is a parameter.
The coefficients A2,A3,A4,A5 are those of [1] with z=1/L.

\\-------------------------------------------- Start cut
La(n) = solve(x=0,10, n - (Pi*exp(x) + 3/4)*x );  \\L=log(a)

G(y)  = 3/4*(y - log(1+y)) + n*log(1+z*log(1+y)) - n*z*y;
       \\ G = [0=(n*z - (Pi*a+3/4))*y] + A2*y^2 + A3*y^3 ...
A2 = -n*z*Pol([1/2,1/2],z) + 3/4/2;
A3 =  n*z*Pol([1/3,1/2,1/3],z) - 3/4/3;
A4 = -n*z*Pol([1/4,1/2,11/24,1/4],z) + 3/4/4;
A5 =  n*z*Pol([1/5,1/2,7/12,5/12,1/5],z) - 3/4/5;
A6 = -n*z*Pol([1/6,1/2,17/24,5/8,137/360,1/6],z) + 3/4/6;
maclv(G,6) == [A2, A3, A4, A5, A6]

       \\ H = exp(A2*y^2)*(1 + B3*y^3 + B4*y^4 + ...)
H(y) = exp(G(y) - A2*y^2);
maclv(H,6) == [1, 0, 0, A3, A4, A5, A6+A3^2/2!]
\\-------------------------------------------- End paste

3. Now substitute x=(1+y)a in F(n) with infinite y-limits
so that the odd powers of y do not contribute anymore :
F(n) ~=~ a*f(a)*\int_(-oo,oo) exp(A2*y^2)*(1 + B4*y^4 ...)
         = ...( 1 + 3*B4/(2*A2)^2 - 15*B6/(2*A2)^3 ...)
         = ...( 1 + b1/n + b2/n^2 + ...)
Our value of b(1)[2] is identical to b1(L) in [1, Th. 9].

\\-------------------------------------------- Start cut
ff(k)   = prod(j=1,k,2*j-1);
b(M)  = my(Bv=maclv(H,6*M)); concat(1,maclv(x->subst(\
             sum(j=2,3*M,ff(j)*Bv[2*j+1]/(-2*A2)^j),n,1/x),M));
subst(b(1),z,1/L) == [ 1, Pol([2,9,16,6,2],L)/24/(L+1)^3 ]

                   \\ Main saddle point formula of order M [1, Th.9]
Fsd(n,M=1) = { my(L=La(n), z=1/L, bv=eval(b(M)));
                        sqrt(2*Pi) *  L^n / sqrt(n*z*(1+z)-3/4) *
                        exp(L/4-n/L+3/4) * subst(Polrev(bv,x),x,1/n);}

                   \\ Asymptotic MacLaurin coefficient formula of order M
Xika(k,M=1) = if(!k, Xi(0), if(k==2, 1/2, \
                        Fsd(k-2,M)/(k-2)!/2^(k-2)) - Fsd(k,M)/k!/2^(k+2) );
\\-------------------------------------------- End paste

4. Relative errors for Xika(k,M) decrease with orders M = 1,2,3
and this asymptotic formula underestimates exact values,
for all the tested coefficients. Output:
    [1.86E-6, 5.92E- 9, 1.71E-11]
    [3.50E-8, 1.64E-11, 7.01E-15]
    [1.03E-8, 2.65E-12, 6.23E-16]

\\-------------------------------------------- Start cut
rer(x,y) = 1 - y/x;
K  = 200;
Kv = K*[1..5]; 
Xv = apply(XikR,Kv);
vecmax(apply(k->rer(Xv[k/K],Xika(k,1)),Kv)) < 3E-5
vecmax(apply(k->rer(Xv[k/K],Xika(k,2)),Kv)) < 3E-7
vecmax(apply(k->rer(Xv[k/K],Xika(k,3)),Kv)) < 3E-9
Mv = [1..3];
apply(M->rer(1.788323906851465527E- 2331,Xika( 1000,M)),Mv)
apply(M->rer(6.634931785662221226E-31412,Xika(10000,M)),Mv)
apply(M->rer(1.546757970209095945E-67925,Xika(20000,M)),Mv)
\\-------------------------------------------- End paste

------------------------------ Features used [4,5]
Pi, gamma, zeta, exp, log, cos, sqrt, sum, sumpos, prod, subst, 
intnuminit, intnum, solve, O(), f(series), Pol, Polrev, Vec, 
normlp, exponent, !x, k!, vecmax, apply, concat

------------------------------ References
[1] Michael Griffin, Ken Ono, Larry Rolen, and Don Zagier 2019
    Jensen polynomials for the Riemann zeta function and other sequences
    PNAS 116, no. 23, 11103–11110 http://www.pnas.org/content/116/23/11103
    http://phys.org/news/2019-05-mathematicians-revive-abandoned-approach-riemann.html
[2] Michael Griffin, Ken Ono, Larry Rolen, and Don Zagier 2018
    Jensen-Polya program for the Riemann hypothesis and related problems
    http://people.oregonstate.edu/~petschec/ONTD/Talk1.pdf
[3] Enrico Bombieri 2019 New progress on the zeta function
     http://www.pnas.org/content/116/23/11085
[4] Tutorial:             http://pari.math.u-bordeaux.fr/pub/pari/manuals/2.11.1/tutorial.pdf
[5] Reference card: http://pari.math.u-bordeaux.fr/pub/pari/manuals/2.11.1/refcard.pdf

Jacques Gélinas