Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
Drexel University or The Math Forum.
|
|
|
|
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 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 .
z2 = t+g means the final sum is comprised of (A), (B) and
the expression:
exp((1-s)*log(N))/(s-1) + exp(-s*log(N))/2 .
I get 6 min 32.6 sec with N-1 = 5199 terms in (B)
[i.e., N = 5200] and khyber = 5538 (alias 'nu') even-indexed
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 E-41161
? abs(q)*abs(bernreal(2*khyber))
%227 = 3.6961736962624999220164898087836688496 E-10014
? 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.
|
|
|
|
