To numerically verify R.H. up to some maximum height T, people often use the Riemann-Siegel formula
Wikipedia article on R.S. or Z function
for at least some of the Zeta/Z computations on zeros.
This method is faster than using a method called Euler-MacLaurin summation, although the latter is preferred for dozens or hundreds of digits accuracy.
I’m in the process of duplicating the results of R.P. Brent at ANU in Canberra, Australia.
It was published in Math. Comp. in 1979, vol. 33, no. 148 October 1979:
link to AMS page on the 1979 paper
It explains the details on Gram points, good and bad, the Gram blocks (as they have become known), Rosser’s rule and Gram blocks where Rosser’s rule fails with “missing zeros”, and so on.
My program is written in the C programming language, and is roughly 1200 lines long.
Sample output:
We didn’t find all the missing zeros in the next Gram block
index1 = 69784846
index2 = 69784847
success = 0
We now look in the preceding Gram Block
Gram point guess = 30461845.1127759444134426
Gram point guess = 30461845.5209311383532622
We found the 2 missing zeros in the preceding Gram block
index1 = 69784843
index2 = 69784844
The updated zeros count is: 75000000
real 667m28.360s // wall clock time.
user 667m15.456s
sys 0m3.475s
As I recall, I used the T value of 32,585,736.4 which corresponds to the location of Gram point g_{74,999,999};
Convention is that g_0 ~= 18.0 , so that [g_0, g_1)
Nearest integer to the n-th Gram point
contains one zeta zero, the second, with the first at about 14.13 . On average, one has k+1 zeta zeros from height 0 to the height of g_k . I’m confident my program is working quite well, as results are corroborated by the work of others, especially Brent and some co-authors on larger computations that were done later, in the 1980s .
I’m running a second job covering the Gram points g_{74,999,999} , where the zero count reached 75 million, to the Gram point G_{199,999,999} where the zeros count should be 200,000,000 +/- 1 at most.
cf.: file riemannsiegel3.gp in /home/david
https://www.ams.org/journals/mcom/1982-39-160/S0025-5718-1982-0669660-1/
The authors set T to be Gram point 200,000,000 and find
200,000,001 zeros up to height g_{200,000,000}.
For record verifications of RH (say to g_{10^13} ), it’s unwieldy to use nothing more sophisticated than the Riemann-Siegel forumula. Xavier Gourdon used the Odlydzko-Schonhage method in his unpublished (or not peer-reviewed) paper/work , as well as the fast multipole (or Greengard-Rokhlin method). He also used the FFT and other devices.
So, my program’s goal is to give exact counts of zeta zeros on the critical line in 0< Im(s) < T, for T from 10^7 to 10^9, and beyond, if I can understand and implement the O.S. method or the Greengard-Rokhlin method.
email contact:
Below, I post user-defined functions in PARI/gp:
PARI/GP Development Headquarters
source: the file
N =
(X)->floor(u(X)^2)
Psi =
(Z)->cos(2*Pi*(Z^2-Z-1.0/16))/cos(2*Pi*Z)
R1 =
(X)->(-1)^(N(X)-1)*Psi(p(X))/u(X)
R2 =
(X)->(-1)^(N(X)-1)*(Psi(p(X))/u(X)-(1/(96*Pi^2))*d3psi(p(X))/(u(X)^3))
R3 =
(X)->(-1)^(N(X)-1)*(Psi(p(X))/u(X)-(1/(96*Pi^2))*d3psi(p(X))/(u(X)^3)+(d2psi(p(X))/(64*(Pi^2))+d6psi(p(X))/(18432*(Pi^4)))/(u(X)^5))
R4 =
(X)->(-1)^(N(X)-1)*(Psi(p(X))/u(X)-(1/(96*Pi^2))*d3psi(p(X))/(u(X)^3)+(d2psi(p(X))/(64*(Pi^2))+d6psi(p(X))/(18432*(Pi^4)))/(u(X)^5)-(d1psi(p(X))/(64*(Pi^2))+d5psi(p(X))/(3840*(Pi^4))+d9psi(p(X))/(5308416*(Pi^6)))/(u(X)^7))
Z =
(X)->2*sum(Y=1,floor(sqrt(X/(2*Pi))),cos(rstheta(X)-X*log(Y))/sqrt(Y))
d1psi =
(Z)->(-0.5*Psi(Z-epsilon)+0.5*Psi(Z+epsilon))/epsilon
d2psi =
(Z)->(Psi(Z-epsilon)-2*Psi(Z)+Psi(Z+epsilon))/(epsilon^2)
d3psi =
(Z)->(-0.5*Psi(Z-epsilon*2)+Psi(Z-epsilon)-Psi(Z+epsilon)+0.5*Psi(Z+epsilon*2))/(epsilon^3)
d5psi =
(Z)->(-0.5*Psi(Z-3*epsilon)+2*Psi(Z-2*epsilon)-2.5*Psi(Z-epsilon)+2.5*Psi(Z+epsilon)-2*Psi(Z+2*epsilon)+0.5*Psi(Z+3*epsilon))/(epsilon^5)
d6psi =
(Z)->(Psi(Z-3*epsilon)-6*Psi(Z-2*epsilon)+15*Psi(Z-epsilon)-20*Psi(Z)+15*Psi(Z+epsilon)-6*Psi(Z+2*epsilon)+Psi(Z+3*epsilon))/(epsilon^6)
d9psi =
(Z)->(-Psi(Z-5*epsilon)+8*Psi(Z-4*epsilon)-27*Psi(Z-3*epsilon)+48*Psi(Z-2*epsilon)-42*Psi(Z-epsilon)+42*Psi(Z+epsilon)-48*Psi(Z+2*epsilon)+27*Psi(Z+3*epsilon)-8*Psi(Z+4*epsilon)+Psi(Z+5*epsilon))/(2*(epsilon^9))
p =
(X)->u(X)^2-N(X)
rs1 =
(X)->Z(X)+R1(X)
rs2 =
(X)->Z(X)+R2(X)
rs3 =
(X)->Z(X)+R3(X)
rs4 =
(X)->Z(X)+R4(X)
rstheta =
(X)->(X/2.0)*log(X/(2.0*Pi))-X/2.0-Pi/8.0+1/(48.0*X)+7.0/(5760.0*X^3)+31.0/(80640*X^5)+381.0/(1290240*X^7)
u =
(X)->(X/(2*Pi))^(0.25)
?
====================================
? epsilon
%405 = 1.00000000000000000000000000000000000000000000 E-5
