Yasuaki Honda at google dot com
[1] Wikipedia von Mangoldt function, https://en.wikipedia.org/wiki/Von_Mangoldt_function
[2] Conrey, The Riemann Hypothesis https://www.ams.org/notices/200303/fea-conrey-web.pdf
[4] The LMFDB Collaboration, The L-functions and Modular Forms Database, home page of the Zeros of zeta(s), https://www.lmfdb.org/zeros/zeta/?limit=200&N=1, 2020 , [Online; accessed 10 October 2020].
scale:1000000;
mangoldt(n)='if n=1 then 1 elseif n=p^k then log(p) else 0;
mangoldt_c(n):=if floor(n)<2 then 0 else block([ifc:ifactors(floor(n))],if length(ifc)=1 then log(ifc[1][1]) else 0)$
(matchdeclare(n,numberp),tellsimp(mangoldt(n),mangoldt_c(n)))$
'psi(x)=sum(mangoldt(n),n,1,floor(x));
kill(psi_c)$
psi_c[n]:=if n=1 then float(mangoldt(1)) else psi_c[n-1]+float(mangoldt(n))$
psi(x):=psi_c[floor(x)]$
for i:1 thru scale do psi(i)$
'integrate(('psi(x)-x)*x^(-1/2+%i*t),x,0,'scale);
xres:0.002$
xlist:exp(makelist(i,i,0,float(log(scale)),xres))$
tmax:60;
tres:0.1;
psilist: map(psi,xlist)-xlist$
innerproduct(list1,list2):=if not(length(list1)=length(list2)) then error("length of arg1 and arg2 does not match.") else
block([c:length(list1),total:0,t1,t2],
for i:1 thru c do (t1:pop(list1),t2:pop(list2),total:total+t1*t2),
return(total))$
xxlist(t):=imagpart(xlist^(-0.5+%i*t))$
primepeek():=block([res:[],c:0],for t:0 step tres thru tmax do (
if mod(c,100)=0 then print(t), c:c+1,
res:append(res,[innerproduct(xxlist(t),psilist)])),return(res))$
compile(all);
warning: encountered undefined variable psi_c in translation. warning: encountered undefined variable xlist in translation. warning: encountered undefined variable tres in translation. warning: encountered undefined variable tmax in translation. warning: encountered undefined variable psilist in translation.
errlist:primepeek()$
0 9.99999999999998 20.00000000000001 30.00000000000016 40.0000000000003 50.00000000000044
plot2d([discrete,makelist(i*tres,i,0,length(errlist)-1),errlist/length(xlist)]);
/*
The LMFDB Collaboration, The L-functions and Modular Forms Database,
home page of the Zeros of zeta(s),
https://www.lmfdb.org/zeros/zeta/?limit=200&N=1, 2020 , [Online; accessed 10 October 2020].
*/
img_rho:[
14.1347251417346937904572519835625,
21.0220396387715549926284795938969,
25.0108575801456887632137909925628,
30.4248761258595132103118975305840,
32.9350615877391896906623689640747,
37.5861781588256712572177634807053,
40.9187190121474951873981269146334,
43.3270732809149995194961221654068,
48.0051508811671597279424727494277,
49.7738324776723021819167846785638,
52.9703214777144606441472966088808,
56.4462476970633948043677594767060,
59.3470440026023530796536486749922
]$
/* The rest 187 zeros are omitted as we draw graph up to 1/2+60*%i. */
plot2d([[discrete,makelist(i*tres,i,0,length(errlist)-1),errlist/length(xlist)],[discrete,makelist([img_rho[i],0],i,1,13)]],
[style,[lines],[points,1]],[legend,false]);