Just another site

Experimenting with Euler-McLaurin summation for zeta function

If we look at MathWorld at the page for the Bernoulli numbers, we find one analytic formula for the even Bernoulli numbers (eqn. 41 or 42).

I’ve been using this analytic expression in formula 1 of Section 6.4 of Edwards’ book on the Riemann zeta function.  I can do this with PARI/gp.  One term that requires a lot of time to

compute is sum_{n = 1, … N-1} n^(-s)  ;  the analytic formula for the Bernoulli numbers doesn’t help with that “partial sum” term.  Where the analytic formula will help is in the terms involving the Bernoulli numbers.  For increased precision, one needs both more terms in the “partial sum” and more Bernoulli-number terms. So I’m experimenting. “khyber” stands for Edwards’ “nu”, roughly the number of even Bernoulli numbers included.


? khyber

%185 = 37000
? N=N+1000

%186 = 16300   // value of N

? t=0;

  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 ; 
  g = sum(X=1,N-1,exp(-s*log(X))); 
  z2 = t+g;

time = 3h, 22min, 34,486 ms.
? floor(log(abs(z2))/log(10))

%188 = -42576   // 42,576 decimals Ok …

\precision : 50,000 digits [ in computations ]


Written by meditationatae

January 1, 2014 at 3:11 pm

Posted in History

Tagged with , ,

4 Responses

Subscribe to comments with RSS.

  1. The analytic formula for the even Bernoulli numbers is at Equation (41) on the MathWorld page linked to just here:


    January 1, 2014 at 3:47 pm

    • 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) );
      /\/\/\ Even-index Bernoulli number term.

      t = t + exp((1-s)*log(N))/(s-1) + exp(-s*log(N))/2 ;
      /\/\/\ term requiring very little work.

      g = sum(X=1,N-1,exp(-s*log(X)));
      /\/\/\ partial sum , i.e. sum_{n=1 … N-1} n^(-s) .

      z2 = t+g;
      /\/\/\ Euler-MacLaurin estimate or formula for zeta(s), typically
      Re(s) = 1/2 in my computations, and in fact s is very very close
      (to within 1/10^50000, say) of the first non-trivial zeta zero with Im(s)>0 ,
      s ~= 1/2 + i*14.1etc .


      January 1, 2014 at 4:01 pm

    • I should say that Simon Plouffe at one time used zeta at 2, 4, 6, 8, 10 approximated closely enough, in order to compute the even-indexed Bernoulli numbers, the B_{2k}, for k = 1, 2, 3, …
      Based on his idea, I used “low-precision” zeta at the {2k}’s to approximate the B_{2k}’s as much as needed in , say, 50,000 or 100000- digit precison values of rho_1, the imaginary part of the first non-trivial zeta zero, in the upper-half of the complex plane.
      In the Euler-MacLaurin approximation to high-high-precision zeta,
      as the {2k}’s increase to 30,000 or even say 60,000 ,
      the multiplicand (?) in Euler-MacLaurin (ref: Edwards, Riemann Zeta Function)
      get “extremely” small;
      For that reason 100,000-decimals precision in e.g. B_{60000} is not needed:
      it gets multiplied by a very very very small number.
      This was borne-out in later posts and rho_1 of zeta computations done in early to mid 2014.

      It remains that, to my knowledge, a faster method to compute zeta (special-purpose) was one of the topics of an article by a computational programmer from Europe, known to Doctor or Professor Luschny of UBC or Simon Fraser University in Vancouver, Canada.

      My own (perhaps un-published) value for rho_1 to 100,000Decimals did agree with that of the European programmer;


      Please look up:

      “Fredrik Johansson Arb C library ”

      [end insert]

      but his (Fredrik Johansson’s ) program, or rather the C Arb library,
      is really fast, and therefore for “just zeta” (high-precision),
      I think it’s recommendable. (but it can do e.g. Hurwitz (?) zeta functions and, I think, a lot more as well).


      September 23, 2014 at 5:18 pm

  2. Reblogged this on meditationatae and commented:

    testing 123..


    May 29, 2014 at 7:46 am

Comments are closed.

%d bloggers like this: