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


Math Forum
»
Discussions
»
sci.math.*
»
sci.math
Notice: We are no longer accepting new posts, but the forums will continue to be readable.
Topic:
PARI/gp wins my Riemann zeta speed contest, so far
Replies:
9
Last Post:
Feb 15, 2013 10:43 PM




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


On 02/15/2013 06:43 AM, David Bernier wrote: > 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 .
z2 = t+g means the final sum is comprised of (A), (B) and the expression: exp((1s)*log(N))/(s1) + exp(s*log(N))/2 .
I get 6 min 32.6 sec with N1 = 5199 terms in (B) [i.e., N = 5200] and khyber = 5538 (alias 'nu') evenindexed Bernoulli numbers. [ Edwards' 'nu' in section 6.4 of his book].
For large nu, and s very close to 1/2 + i*14.13, the error term is bounded in z2 = t+g (within a small factor) by  B_2nu * q_{2nu}  , where B_2nu * q_{2nu} is the term with the Bernoulli number B_2nu in (A).
I get: B_{2*5538} = B_11076 ~= 10^31147 ,
so how many significant digits of B_11076 are needed in a 10,000 decimal precision value for zeta(1/2 + i*14.13...) ?
It seems it would be very few, if I've got the bounds correct (say 1 to 20 significant digits of B_11076 ... )
? abs(bernreal(2*khyber)) %225 = 2.6420009078322364970094272369181821282 E31147
? abs(q) %226 = 1.3990054603331733520906578588928549721 E41161
? abs(q)*abs(bernreal(2*khyber)) %227 = 3.6961736962624999220164898087836688496 E10014
? N %228 = 5200
? khyber %229 = 5538
David Bernier  dracut:/# lvm vgcfgrestore File descriptor 9 (/.console_lock) leaked on lvm invocation. Parent PID 993: sh Please specify a *single* volume group to restore.



