Ticket #8135: riemann_pi.gp

File riemann_pi.gp, 4.9 KB (added by kevin.stueve, 9 years ago)

pari-gp code from Tomás Oliveira e Silva

Line 
1/*
2** Riemann's formula for pi(x)
3**
4** T. Oliveira e Silva, January 2010
5*/
6
7/*
8** enable the computation of pi(x) for x up to 10^8
9*/
10default(primelimit,10^8);
11
12/*
13** computation of li(x^{1/2+I t})+li(x^{1/2-I t})
14**
15** li(x^\rho) = x^\rho/u*(1+1/u+2!/u^2+3!/u^3+...),
16** with \rho=1/2+it and u=\rho\log(x)
17*/
18li(x,t)=
19{ local(u,s0,s1,s2,lx,tol,k);
20  u=(1/2+I*t)*log(x); /* log(x^{1/2+I t}) */
21  s0=2*x^(1/2+I*t)/u;
22  tol=1e-9/abs(s0); /* we desire an error of less than 10^{-9} */
23  s1=s2=1;
24  k=1;
25  while(1,
26    if(k>100,
27      error("li: unable to attain the desired precision");
28    );
29    s2*=k/u;
30    s1+=s2;
31    if(abs(s2)<tol,
32      break();
33    );
34    k++;
35  );
36  return(real(s0*s1));
37}
38
39/*
40** the imaginary part of the first 100 zeros of the zeta function on the critical line
41*/
42zz=vector(100);
43zz[  1]= 14.13472514173469379046;
44zz[  2]= 21.02203963877155499263;
45zz[  3]= 25.01085758014568876321;
46zz[  4]= 30.42487612585951321031;
47zz[  5]= 32.93506158773918969066;
48zz[  6]= 37.58617815882567125722;
49zz[  7]= 40.91871901214749518740;
50zz[  8]= 43.32707328091499951950;
51zz[  9]= 48.00515088116715972794;
52zz[ 10]= 49.77383247767230218192;
53zz[ 11]= 52.97032147771446064415;
54zz[ 12]= 56.44624769706339480437;
55zz[ 13]= 59.34704400260235307965;
56zz[ 14]= 60.83177852460980984426;
57zz[ 15]= 65.11254404808160666088;
58zz[ 16]= 67.07981052949417371448;
59zz[ 17]= 69.54640171117397925293;
60zz[ 18]= 72.06715767448190758252;
61zz[ 19]= 75.70469069908393316833;
62zz[ 20]= 77.14484006887480537268;
63zz[ 21]= 79.33737502024936792276;
64zz[ 22]= 82.91038085408603018316;
65zz[ 23]= 84.73549298051705010574;
66zz[ 24]= 87.42527461312522940653;
67zz[ 25]= 88.80911120763446542368;
68zz[ 26]= 92.49189927055848429626;
69zz[ 27]= 94.65134404051988696660;
70zz[ 28]= 95.87063422824530975874;
71zz[ 29]= 98.83119421819369223332;
72zz[ 30]=101.31785100573139122879;
73zz[ 31]=103.72553804047833941640;
74zz[ 32]=105.44662305232609449367;
75zz[ 33]=107.16861118427640751512;
76zz[ 34]=111.02953554316967452466;
77zz[ 35]=111.87465917699263708561;
78zz[ 36]=114.32022091545271276589;
79zz[ 37]=116.22668032085755438216;
80zz[ 38]=118.79078286597621732298;
81zz[ 39]=121.37012500242064591895;
82zz[ 40]=122.94682929355258820082;
83zz[ 41]=124.25681855434576718473;
84zz[ 42]=127.51668387959649512428;
85zz[ 43]=129.57870419995605098577;
86zz[ 44]=131.08768853093265672357;
87zz[ 45]=133.49773720299758645013;
88zz[ 46]=134.75650975337387133133;
89zz[ 47]=138.11604205453344320019;
90zz[ 48]=139.73620895212138895045;
91zz[ 49]=141.12370740402112376194;
92zz[ 50]=143.11184580762063273941;
93zz[ 51]=146.00098248676551854740;
94zz[ 52]=147.42276534255960204952;
95zz[ 53]=150.05352042078488035143;
96zz[ 54]=150.92525761224146676185;
97zz[ 55]=153.02469381119889619826;
98zz[ 56]=156.11290929423786756975;
99zz[ 57]=157.59759181759405988753;
100zz[ 58]=158.84998817142049872417;
101zz[ 59]=161.18896413759602751944;
102zz[ 60]=163.03070968718198724331;
103zz[ 61]=165.53706918790041883004;
104zz[ 62]=167.18443997817451344096;
105zz[ 63]=169.09451541556882148951;
106zz[ 64]=169.91197647941169896670;
107zz[ 65]=173.41153651959155295985;
108zz[ 66]=174.75419152336572581338;
109zz[ 67]=176.44143429771041888889;
110zz[ 68]=178.37740777609997728583;
111zz[ 69]=179.91648402025699613934;
112zz[ 70]=182.20707848436646191541;
113zz[ 71]=184.87446784838750880096;
114zz[ 72]=185.59878367770747146653;
115zz[ 73]=187.22892258350185199164;
116zz[ 74]=189.41615865601693708485;
117zz[ 75]=192.02665636071378654728;
118zz[ 76]=193.07972660384570404740;
119zz[ 77]=195.26539667952923532146;
120zz[ 78]=196.87648184095831694862;
121zz[ 79]=198.01530967625191242492;
122zz[ 80]=201.26475194370378873302;
123zz[ 81]=202.49359451414053427769;
124zz[ 82]=204.18967180310455433072;
125zz[ 83]=205.39469720216328602521;
126zz[ 84]=207.90625888780620986150;
127zz[ 85]=209.57650971685625985284;
128zz[ 86]=211.69086259536530756391;
129zz[ 87]=213.34791935971266619064;
130zz[ 88]=214.54704478349142322294;
131zz[ 89]=216.16953850826370026587;
132zz[ 90]=219.06759634902137898568;
133zz[ 91]=220.71491883931400336912;
134zz[ 92]=221.43070555469333873210;
135zz[ 93]=224.00700025460433521173;
136zz[ 94]=224.98332466958228750378;
137zz[ 95]=227.42144427967929131046;
138zz[ 96]=229.33741330552534810776;
139zz[ 97]=231.25018870049916477381;
140zz[ 98]=231.98723525318024860377;
141zz[ 99]=233.69340417890830064070;
142zz[100]=236.52422966581620580248;
143
144
145/*
146** pi(x)=li(x)
147**      -sum_{\rho} li(x^{\rho})
148**      -log(2)
149**      -sum_{k=2}^{\infty} pi(x^{1/k})
150**      +\int_x^\infty dt/(t*(t^2-1)*log(t))
151**
152** given the above prime limit, works for x up to 10^{16}
153*/
154riemann_pi(x)=
155{ local(k,y,z);
156  y=-eint1(-log(x)); /* $li(x)=int_0^x du/\log u$ */
157  y-=log(2);
158  k=2;
159  /* the pi(x^(1/k)) terms */
160  while(1,
161    z=floor((x+0.1)^(1/k)); /* add 0.1 to avoid possible roundoff errors */
162    if(z<2,
163      break();
164    );
165    y-=primepi(z)/k;
166    k++;
167  );
168  /* \int_x^\infty dt/(t*(t^2-1)*log(t)) */
169  /* this term is always smaller than 1/(2*x^2*log(x)), and so can usually be ignored */
170  k=1;
171  while(1,
172    z=incgam(0,2*k*log(x));
173    y+=z;
174    k++;
175    if(z<1e-9,
176      break;
177    );
178  );
179  /* the zeta zeros */
180  for(k=1,100,
181    y-=li(x,zz[k]);
182  );
183  /* done */
184  return(y);
185}