#
# at most 600 nontrivial interpolated "zeros" (only 200 actual zeros)
# 08/27/19 original code by Daniel Hutama (computes Riemann's R(x)),
# modified by Darrell Cox to support linearly interpolated Riemann spectrum
#
import numpy as np
import cmath
import matplotlib.pyplot as plt
import math
import cmath
import mpmath as mp
from numpy import polyfit
plt.rcParams['text.usetex'] = True
plt.rcParams['text.latex.unicode'] = True
import scipy
from sage.functions.spike_function import SpikeFunction
from scipy.interpolate import interp1d
from scipy import stats
import scipy.integrate as integrate
import scipy.special as special
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
from sage.libs.lcalc.lcalc_Lfunction import *
# These formulas are the components of the explicit formula.
x = var('x')
def Riemann_main(x, numbersums, numberzeros, numbertrivials):
return sum([moebius(n)/n*Ei((1/n)*log(x)) for n in xrange(1, numbersums+1)])
def zeta_nontrivial_correction(x, numbersums, numberzeros, numbertrivials):
# zz = zeta_zeros()[0:numberzeros]
# linearly interpolated zeros (two values between successive zeros)
zz=[14.1347251417347000, 16.4304966407469880, 18.7262681397592740, 21.0220396387715600,
22.3516456192296040, 23.6812515996876480, 25.0108575801456890, 26.8155304287169610,
28.6202032772882330, 30.4248761258595090, 31.2616046131527380, 32.0983331004459630,
32.9350615877391920, 34.4854337781013530, 36.0358059684635140, 37.5861781588256820,
38.6970251099329520, 39.8078720610402210, 40.9187190121474980, 41.7215037684033310,
42.5242885246591630, 43.3270732809150020, 44.8864324809990550, 46.4457916810831080,
48.0051508811671610, 48.5947114133355380, 49.1842719455039160, 49.7738324776723000,
50.8393288110196830, 51.9048251443670670, 52.9703214777144570, 54.1289635508307680,
55.2876056239470800, 56.4462476970633920, 57.4131797989097090, 58.3801119007560270,
59.3470440026023520, 59.8419555099381740, 60.3368670172739950, 60.8317785246098100,
62.2587003657670760, 63.6856222069243430, 65.1125440480816020, 65.7682995418857960,
66.4240550356899890, 67.0798105294941680, 67.9020075900541120, 68.7242046506140550,
69.5464017111739850, 70.3866536989432920, 71.2269056867125980, 72.0671576744819050,
73.2796686826825830, 74.4921796908832620, 75.7046906990839260, 76.1847404890142170,
76.6647902789445080, 77.1448400688747990, 77.8756850526663270, 78.6065300364578550,
79.3373750202493680, 80.5283769648615930, 81.7193789094738180, 82.9103808540860290,
83.5187515628963690, 84.1271222717067100, 84.7354929805170510, 85.6320868580531140,
86.5286807355891770, 87.4252746131252250, 87.8865534779616410, 88.3478323427980570,
88.8091112076344590, 90.0367072286091310, 91.2643032495838040, 92.4918992705584910,
93.2117141938789620, 93.9315291171994320, 94.6513440405198880, 95.0577741030950280,
95.4642041656701680, 95.8706342282453080, 96.8574875582281010, 97.8443408882108940,
98.8311942181936870, 99.6600798140395910, 100.4889654098854900, 101.3178510057313800,
102.1204133506470400, 102.9229756955627000, 103.7255380404783400, 104.2992330444276000,
104.8729280483768500, 105.4466230523260900, 106.0206190963095300, 106.5946151402929700,
107.1686111842764000, 108.4555859705741600, 109.7425607568719200, 111.0295355431696700,
111.3112434211106600, 111.5929512990516500, 111.8746591769926400, 112.6898464231459900,
113.5050336692993500, 114.3202209154527100, 114.9557073839209900, 115.5911938523892800,
116.2266803208575500, 117.0813811692304400, 117.9360820176033300, 118.7907828659762100,
119.6505635781243600, 120.5103442902725100, 121.3701250024206500, 121.8956930994646200,
122.4212611965086000, 122.9468292935525800, 123.3834923804836400, 123.8201554674147100,
124.2568185543457700, 125.3434403294293400, 126.4300621045129100, 127.5166838795965000,
128.2040239863830200, 128.8913640931695600, 129.5787041999560600, 130.0816989769482700,
130.5846937539404800, 131.0876885309326700, 131.8910380882876400, 132.6943876456426200,
133.4977372029976000, 133.9173280531230100, 134.3369189032484300, 134.7565097533738800,
135.8763538537604100, 136.9961979541469400, 138.1160420545334400, 138.6560976870627500,
139.1961533195920700, 139.7362089521213900, 140.1987084360879700, 140.6612079200545500,
141.1237074040211300, 141.7864202052209600, 142.4491330064207900, 143.1118458076206300,
144.0748913673356000, 145.0379369270505700, 146.0009824867655100, 146.4749101053635300,
146.9488377239615600, 147.4227653425596100, 148.2996837019680400, 149.1766020613764600,
150.0535204207848800, 150.3440994846037300, 150.6346785484225800, 150.9252576122414700,
151.6250696785606100, 152.3248817448797500, 153.0246938111988900, 154.0540989722118800,
155.0835041332248600, 156.1129092942378800, 156.6078034686899500, 157.1026976431420200,
157.5975918175940600, 158.0150572688695500, 158.4325227201450400, 158.8499881714205100,
159.6296468268123500, 160.4093054822041900, 161.1889641375960300, 161.8028793207913600,
162.4167945039866900, 163.0307096871820000, 163.8661628540881500, 164.7016160209942900,
165.5370691879004100, 166.0861927846584400, 166.6353163814164600, 167.1844399781745100,
167.8211317906392700, 168.4578236031040300, 169.0945154155688200, 169.3670024368497900,
169.6394894581307500, 169.9119764794116900, 171.0784964928049900, 172.2450165061982800,
173.4115365195915500, 173.8590881875162800, 174.3066398554410100, 174.7541915233657300,
175.3166057814806300, 175.8790200395955300, 176.4414342977104300, 177.0867587905069300,
177.7320832833034400, 178.3774077760999700, 178.8904331908189800, 179.4034586055379900,
179.9164840202570000, 180.6800155082934800, 181.4435469963299600, 182.2070784843664600,
183.0962082723734700, 183.9853380603804900, 184.8744678483875000, 185.1159064581608200,
185.3573450679341500, 185.5987836777074700, 186.1421633129722600, 186.6855429482370400,
187.2289225835018600, 187.9580012743402100, 188.6870799651785500, 189.4161586560169300,
190.2863245575825600, 191.1564904591481900, 192.0266563607137900, 192.3776797750910800,
192.7287031894683800, 193.0797266038457000, 193.8082832957402100, 194.5368399876347200,
195.2653966795292300, 195.8024250666722700, 196.3394534538153100, 196.8764818409583200,
197.2560911193895100, 197.6357003978207000, 198.0153096762519200, 199.0984570987358800,
200.1816045212198400, 201.2647519437038000, 201.6743661338493900, 202.0839803239949800,
202.4935945141405400, 203.0589536104618800, 203.6243127067832300, 204.1896718031045500,
204.5913469361241300, 204.9930220691437100, 205.3946972021632900, 206.2318844307109200,
207.0690716592585500, 207.9062588878062200, 208.4630091641562300, 209.0197594405062500,
209.5765097168562600, 210.2812940096926000, 210.9860783025289400, 211.6908625953653000,
212.2432148501477500, 212.7955671049302000, 213.3479193597126800, 213.7476278343056000,
214.1473363088985300, 214.5470447834914300, 215.0878760250821800, 215.6287072666729300,
216.1695385082637100, 217.1355577885162800, 218.1015770687688500, 219.0675963490213900,
219.6167038457855900, 220.1658113425498000, 220.7149188393140100, 220.9535144111071100,
221.1921099829002100, 221.4307055546933300, 222.2894704546636700, 223.1482353546340100,
224.0070002546043200, 224.3324417262636500, 224.6578831979229800, 224.9833246695822900,
225.7960312062813000, 226.6087377429803100, 227.4214442796792900, 228.0601006216279800,
228.6987569635766700, 229.3374133055253600, 229.9750051038499500, 230.6125969021745400,
231.2501887004991700, 231.4958708847262000, 231.7415530689532400, 231.9872352531802400,
232.5559582284229200, 233.1246812036656000, 233.6934041789083100, 234.6370126745442800,
235.5806211701802500, 236.5242296658161900, 236.9394266041858500, 237.3546235425555200,
237.7698204809252000, 238.3650395117260200, 238.9602585425268400, 239.5554775733276400,
240.0533709809572800, 240.5512643885869300, 241.0491577962165800, 241.6405291755519000,
242.2319005548872300, 242.8232719342225900, 243.2391474551744500, 243.6550229761263000,
244.0708984970781600, 245.0929290230179400, 246.1149595489577300, 247.1369900748975100,
247.4586567366478400, 247.7803233983981700, 248.1019900601484700, 248.5925565883347000,
249.0831231165209400, 249.5736896447072000, 250.0541090281434700, 250.5345284115797400,
251.0149477950160100, 251.6999607793438300, 252.3849737636716400, 253.0699867479994800,
253.8154099836376600, 254.5608332192758500, 255.3062564549140300, 255.6644088680875200,
256.0225612812610100, 256.3807136944344600, 257.1239556268001300, 257.8671975591657900,
258.6104394915313900, 259.0317619909135900, 259.4530844902957900, 259.8744069896779900,
260.1846328279842700, 260.4948586662905500, 260.8050845045968900, 261.7280209713546200,
262.6509574381123600, 263.5738939048701500, 264.2352132162055200, 264.8965325275409000,
265.5578518388763300, 265.9102258197512500, 266.2625998006261600, 266.6149737815010800,
267.0506208819420500, 267.4862679823830300, 267.9219150828240600, 268.6047597298818900,
269.2876043769397300, 269.9704490239976200, 270.4783178965467400, 270.9861867690958700,
271.4940556416449900, 272.1492401572310800, 272.8044246728171600, 273.4596091884033100,
274.1689036753834900, 274.8781981623636700, 275.5874926493438600, 275.8756782672735400,
276.1638638852032300, 276.4520495031329100, 277.0516141787026100, 277.6511788542723000,
278.2507435298419400, 278.5769126624763300, 278.9030817951107200, 279.2292509277451700,
280.3078722068474900, 281.3864934859498100, 282.4651147650520800, 282.7138050877793500,
282.9624954105066200, 283.2111857332338900, 283.7527784824575300, 284.2943712316811700,
284.8359639809047500, 285.4464577749374700, 286.0569515689701900, 286.6674453630029100,
287.0822704091426700, 287.4970954552824200, 287.9119205014221700, 288.4678986440210300,
289.0238767866198900, 289.5798549292188100, 290.3353337291683300, 291.0908125291178400,
291.8462913290674100, 292.4170055991637000, 292.9877198692599900, 293.5584341393562800,
294.0274126326593800, 294.4963911259624800, 294.9653696192655200, 295.1679980391630800,
295.3706264590606500, 295.5732548789582700, 296.3752622732866400, 297.1772696676150100,
297.9792770619434400, 298.5996267258693700, 299.2199763897953100, 299.8403260537213000,
300.4433258565455800, 301.0463256593698700, 301.6493254621941600, 301.9984668379984200,
302.3476082138026900, 302.6967495896069000, 303.4192901733570100, 304.1418307571071300,
304.8643713408573000, 305.1525517612504900, 305.4407321816436800, 305.7289126020368100,
306.2257737774145900, 306.7226349527923600, 307.2194961281700800, 308.1828184676807000,
309.1461408071913300, 310.1094631467019000, 310.4613559412532600, 310.8132487358046200,
311.1651415303559800, 311.5860280804376400, 312.0069146305193000, 312.4278011806009100,
312.9469626974535600, 313.4661242143062100, 313.9852857311589100, 314.4820625172645200,
314.9788393033701400, 315.4756160894757500, 316.2286793737739000, 316.9817426580720500,
317.7348059423702000, 318.1075720470190000, 318.4803381516678100, 318.8531042563166100,
319.6221142739156000, 320.3911242915145900]
n_dict = {n:[] for n in xrange(1, numbersums+1)}
for n in xrange(1, numbersums+1):
if moebius(n)==1:
n_dict[n] = -sum([(1/n)*2*real(Ei((0.5+I*gamma)/n*log(x))) for gamma in zz])
elif moebius(n)==0:
n_dict[n] = 0
elif moebius(n)==-1:
n_dict[n] = sum([(1/n)*2*real(Ei((0.5+I*gamma)/n*log(x))) for gamma in zz])
return sum(n_dict.values())
def trivial_zeros_correction(x, numbersums, numberzeros, numbertrivials):
n_dict = {n:[] for n in xrange(1, numbersums + 1)}
for n in xrange(1, numbersums+1):
if moebius(n)==1:
n_dict[n] = -sum([(1/n)*Ei((-2*m/n)*log(x)) for m in xrange(1, numbertrivials+1)])
elif moebius(n)==0:
n_dict[n] = 0
elif moebius(n)==-1:
n_dict[n] = sum([(1/n)*Ei((-2*m/n)*log(x)) for m in xrange(1, numbertrivials+1)])
return sum(n_dict.values())
def pi_0(x, numbersums, numberzeros, numbertrivials):
return sum([Riemann_main(x, numbersums, numberzeros, numbertrivials), zeta_nontrivial_correction(x, numbersums, numberzeros, numbertrivials), trivial_zeros_correction(x, numbersums, numberzeros, numbertrivials)])
mindistance=2
maxdistance=10000
nontrivialzeros=600
trivialzeros=200
pi_0_plot = plot(pi_0(x, ceil(log(maxdistance)-log(2))+1, nontrivialzeros, trivialzeros), (x, mindistance, maxdistance), rgbcolor = 'blue', legend_label=r'$\pi_0(x)$, 600 non-trivial "zeros" ')
primes = plot(prime_pi(x), (x, 0, maxdistance), rgbcolor = 'red', legend_label = r'$\pi(x)$')
p = plot(primes + pi_0_plot)
p.set_legend_options(loc=(0.05, 0.8))
p.show(legend_font_size=18, fontsize = 12)
p.save('newrat3.pdf',legend_font_size=18, fontsize=12)