Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
NCTM or The Math Forum.



Re: PARI/gp wins my Riemann zeta speed contest, so far
Posted:
Feb 15, 2013 6:43 AM


On 02/07/2013 01:19 PM, David Bernier wrote: > On 01/21/2013 12:23 AM, David Bernier wrote: >> I had done my own Riemann zeta computations with Bernoulli numbers, >> computing in days a few 10,000 (say 55,000) decimals of >> the imaginary part of the first nontrivial zero >> 1/2 + i*14.13... using PARIgp. >> >> This time, I used PARIgp's own zeta(.): >> >> ? system(date);zeta(t);system(date) >> Sun Jan 20 03:36:28 EST 2013 >> Sun Jan 20 07:10:17 EST 2013 >> >> >> // 20,000 decimals precision >> // t is with 1/10^400 of first nontrivial zero. >> >> 3 hours and 34 minutes for 20017 significant digits >> near 1/2 + i*14.134725141734693790457251983562470270784257 >> >> >> >> ? \p >> realprecision = 20017 significant digits (20000 digits displayed) >> ? a=zeta(t); > [...] > > With realprecision \p set to 40,000 digits, something strange > happened with zeta(.): > > ? zeta(s) > *** at toplevel: zeta(s) > *** ^ > *** zeta: the PARI stack overflows ! > current stack size: 500000000 (476.837 Mbytes) > [hint] you can increase GP stack with allocatemem() > > *** Break loop: type 'break' to go back to GP > break> break > > ? allocatemem(2000000000) > *** Warning: new stack size = 2000000000 (1907.349 Mbytes). > ? zeta(s) > ^C *** at toplevel: zeta(s) > *** ^ > *** zeta: user interrupt after 38h, 47min, 59,768 ms. > > about 39 hours and has not completed. > 500 MB of stack needed. > > I was doublechecking my EulerMaclaurin computations. > > To be continued ... > > David Bernier > >
So, in one gp session, I have rho_1, the first nontrivial zero of the Riemann zeta function, to about 54,000 decimal digits accuracy.
I saved that to a file. In a second PARI/gp session, with 10,000 digits precision, I wanted to find the best N and nu for the EulerMacLaurin formula in section 6.4 of H. Edwards' book "Riemann's Zeta Function".
Best here means where the computation takes as little time as possible.
I think I'm quite close to minimum time for my setup:
? N %46 = 5500
? khyber %47 = 5300
? q=s/(2*N*exp(s*log(N)));
t=0;
q=s/(2*N*exp(s*log(N)));
for(X=1,khyber,t=t+bernfrac(2*X)*q;q = q*((s+2*X1)*(s+2*X)/(N*N))/((2*X+1)*(2*X+2)) ); (A)
t = t + exp((1s)*log(N))/(s1) + exp(s*log(N))/2 ;
g = sum(X=1,N1,exp(s*log(X))); (B)
z2 = t+g;
abs(z2)
%48 = 2.137363631927886601 E10002
time = 6min, 35,192 ms.
(A) is/does the the sum involving the evenindex Bernoulli numbers. 'q' is what gets multiplied by B_{2X}, so q is initialized before starting the computation in (A), and 'q' varies as X goes from 1 to 5300 (term up to B_10600 ).
(B) is/does the sum: sum_{k = 1, 5499} 1/(k^s).
s = 1/2 + i*14.13.... (first nontrivial zeta zero to 10,000 decimals).
With those paramters, the computation of zeta(s) took about 6 min 35 sec .
David Bernier
P.S.: I read SAGE has an improved (modular) Bernoulli number algorithm.
 dracut:/# lvm vgcfgrestore File descriptor 9 (/.console_lock) leaked on lvm invocation. Parent PID 993: sh Please specify a *single* volume group to restore.



