## Experimenting with Euler-McLaurin summation for zeta, Update

When using the Euler-MacLaurin formula to compute the zeta function along the critical line to very high precision, I’ve found that using the analytic expression for the Bernoulli numbers does save time for the term involving the Bernoulli numbers. So, in my computations of zeta near the first non-trivial zeta zero, the lion’s share of the time is now spent in computing the partial sums of the zeta function, following H. Edwards expression (1) in Section 6.4 of his book “Riemann’s Zeta Function”. I use `N’ for Edwards’s `N’, and `khyber’ for Edwards’s `nu’, the Greek letter. In PARI/gp, I read a file containing an integer close to t_0 * 10^90000, where t_0 is the imaginary part of the first non-trivial zeta zero. This allows me to set-up `s’ in “exact” rational arithmetic, viz.:

? \p 90000

realprecision = 90010 significant digits (90000 digits displayed)

? read(round90kay)

%207 1413…… [90k digits]

<snip>

? bb=%207

? tt = bb/(10^90000);

? s = 1/2 + I*tt;

This completes the set-up. Both the real and imaginary parts of `s’ are now rational numbers.

? N

%214 = 35000

So, N = 35,000.

The partial sum of zeta to 35000 terms takes almost 20 hours to compute:

? g = sum(X=1,N-1,exp(-s*log(X)));

time = 19h, 45min, 50,152 ms.

? \p

realprecision = 90010 significant digits (90000 digits displayed)

? N

%214 = 35000

On the other hand, the terms involving the even Bernoulli numbers take a fraction of that time.

My latest run is with nu = 60000, and it took about 6 hours:

? t=0; q=s/(2*N*exp(s*log(N))); q=q/(Pi*Pi);

for(X=1,khyber,t=t+ (((-1)^(X-1))*zeta(2*X))*q;

q = q*((s+2*X-1)*(s+2*X)/(N*N)); q=q/(4*Pi*Pi) );

t = t + exp((1-s)*log(N))/(s-1) + exp(-s*log(N))/2 ;

z2 = t+g;

—

I left space around the even-index Bernoulli number sum for greater legibility. Once again, it is:

for(X=1,khyber,t=t+ (((-1)^(X-1))*zeta(2*X))*q;

q = q*((s+2*X-1)*(s+2*X)/(N*N)); q=q/(4*Pi*Pi) );

This uses an exact analytic formula for the Bernoulli numbers that can be found at Mathworld in their article on Bernoulli numbers. For example, zeta(1000), B_1000 and Pi are all inter-related.

time = 5h, 54min, 10,238 ms.

? floor(log(abs(z2))/log(10))

%223 = -83677

This shows that combining the partial sums in `g’ with the Bernoulli numbers terms stored in `t’, we get 83,677 significant digits.

At this value of `nu’ (khyber), increasing `nu’ by 1000 delivers roughly 560 extra digits of accuracy. I set my next run to:

nu = 71000 .

For a fixed N, as nu increases, eventually there is no gain from futher increasing nu. Edwards mention the Backlund estimate on the error term. One has to know the approximate size of the Bernoulli numbers, but that is quite straightforward.

Fredrik Johansson has an article at arxiv about software he wrote to compute the Hurwitz zeta function to very high accuracy and related topics, such as the Stieltjes constants, “

# Rigorous high-precision computation of the Hurwitz zeta function and its derivatives” .

time = 5h, 54min, 10,238 ms. [ for nu = 60000 ].