Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

Views expressed in these public forums are not endorsed by Drexel University or The Math Forum.


Math Forum » Discussions » sci.math.* » sci.math

Topic: PARI/gp wins my Riemann zeta speed contest, so far
Replies: 9   Last Post: Feb 15, 2013 10:43 PM

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
David Bernier

Posts: 3,401
Registered: 12/13/04
Re: PARI/gp wins my Riemann zeta speed contest, so far
Posted: Feb 15, 2013 6:43 AM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

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 non-trivial zero
>> 1/2 + i*14.13... using PARI-gp.
>>
>> This time, I used PARI-gp'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 non-trivial 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 top-level: 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 top-level: 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 double-checking my Euler-Maclaurin computations.
>
> To be continued ...
>
> David Bernier
>
>


So, in one gp session, I have rho_1, the first non-trivial
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 Euler-MacLaurin 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 set-up:

? 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*X-1)*(s+2*X)/(N*N))/((2*X+1)*(2*X+2)) ); (A)

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

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

z2 = t+g;

abs(z2)


%48 = 2.137363631927886601 E-10002

time = 6min, 35,192 ms.

(A) is/does the the sum involving the even-index 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 non-trivial 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.



Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.