pi-minus-li.sagews

plot(lambda x : prime_pi(x) - Li(x), (2, 100), plot_points=1000)

plot(lambda x : prime_pi(x) - Li(x), (2, 1000), plot_points=1000)

plot(lambda x : prime_pi(x) - Li(x), (2, 10000), plot_points=1000)

︠8e4219a9-d2d3-44aa-b571-968704cc195d︠
B = 20000
v = stats.TimeSeries([Li(n) - prime_pi(n) for n in [2..B]])
v2 = stats.TimeSeries([0]+[v[:i].mean() for i in range(1,len(v))])
v.plot(plot_points=B) + v2.plot(color='red', plot_points=B) + plot(sqrt(x)*(13/sqrt(10000)), (2, B), color='purple')

%time
B = 40000
v = stats.TimeSeries([Li(n) - prime_pi(n) for n in [2..B]]).abs()
v2 = stats.TimeSeries([0]+[v[:i].mean() for i in range(1,len(v))])
v.plot(plot_points=B) + v2.plot(color='red', plot_points=B) + plot(sqrt(2/pi)*sqrt(x)/log(x), (2, B), color='purple')

CPU time: 19.64 s, Wall time: 17.60 s

%time
B = 100000
v = stats.TimeSeries([Li(n) - prime_pi(n) for n in [2..B]]).abs()
v2 = stats.TimeSeries([0]+[v[:i].mean() for i in range(1,len(v))])
︠b42efae5-1954-4787-a2be-5880253c8ff1︠

v.plot(plot_points=B) + v2.plot(color='red', plot_points=B) + plot(sqrt(x)/log(x), (2, B), color='purple')

%time
B = 100000
v = stats.TimeSeries([Li(n) - prime_pi(n) for n in [2..B]]).abs()
v2 = stats.TimeSeries([0]+[v[:i].mean() for i in range(1,len(v))])
︠7e8bd2f4-a308-4f0f-b2ec-b782653fa0f8︠
x=var('x')
v.plot(plot_points=B) + v2.plot(color='red', plot_points=B) + plot(sqrt(2/pi)*sqrt(x/log(x)), (2, B), color='purple')

%time v = stats.TimeSeries([Li(n) - prime_pi(n) for n in [2..B]]).abs()

CPU time: 36.28 s, Wall time: 36.13 s

def prime_pi_time_series(B):
    """
    Return the time series stats.TimeSeries([prime_pi(n) for n in [0..B-1]])
    but computed (vastly) more efficiently.
    """
    x = 0
    w = []
    pp = 0
    for p in prime_range(B):
        w.extend([x]*(p-pp))
        pp = p
        x += 1
    w.extend([x]*(B-pp))
    return stats.TimeSeries(w)

def running_average(v):
    """
    Return the running average of the time series... i.e., the Cesaro sum.
    i.e., stats.TimeSeries([0]+[v[:i].mean() for i in range(1,len(v))]),
    but computed (vastly) more efficiently.
    """
    s = v.sums()
    # now divide s[i] by i+1 for i>=1.
    for i in range(1,len(s)):
        s[i] /= i+1
    return s

running_average(stats.TimeSeries([1,2,3,4]))

[1.0000, 1.5000, 2.0000, 2.5000]

import math

import scipy.special

import mpmath

a = float(5)

%timeit float(mpmath.li(a,offset=True))

625 loops, best of 3: 75.5 µs per loop

%timeit Li(a)

625 loops, best of 3: 145 µs per loop

%time zz = [Li(float(n)) for n in range(B)]
︠fdd339fc-10c3-4d2e-9dc9-d06c7cbef114︠

%time zz = [Li(float(n)) for n in range(B)]

CPU time: 1.71 s, Wall time: 1.72 s

B = 100000
%time v = prime_pi_time_series(B)
%time li_minus_pi = stats.TimeSeries([0,0] + [mpmath.li(i,offset=True)-v[i] for i in range(2,B)])
%time v2 =running_average(li_minus_pi)

CPU time: 0.04 s, Wall time: 0.04 s CPU time: 10.63 s, Wall time: 10.62 s CPU time: 0.70 s, Wall time: 0.70 s

︠911d1322-5fa1-456f-9450-a9dd75692ae0︠
li_minus_pi.plot(plot_points=B) + v2.plot(color='red', plot_points=B) + plot(sqrt(2/pi)*sqrt(x/log(x)), (2, B), color='purple')

B = 250000
%time v = prime_pi_time_series(B)
%time li_minus_pi = stats.TimeSeries([0,0] + [mpmath.li(i,offset=True)-v[i] for i in range(2,B)])
%time v2 =running_average(li_minus_pi)

CPU time: 0.09 s, Wall time: 0.09 s CPU time: 27.17 s, Wall time: 27.16 s CPU time: 1.78 s, Wall time: 1.78 s

li_minus_pi.plot(plot_points=B) + v2.plot(color='red', plot_points=B) + plot(sqrt(2/pi)*sqrt(x/log(x)), (2, B), color='purple')

li_minus_pi=5101648384.7155276
x = 4*10^22
N(sqrt(2/pi)*sqrt(x/log(x))) - li_minus_pi

1.70185084142047e10

li = 783964159852157952242.7155
N(sqrt(2/pi)*sqrt(li)) - li_minus_pi

1.72386086637032e10

︠8a08f9f4-0e4c-40d0-a56c-9ccc259bc5a6︠


v3 = stats.TimeSeries([0]+[v2[i]/sqrt(i-1) for i in [2..len(v2)-1]])

%var x
v.plot(plot_points=B) + v2.plot(color='red', plot_points=B) + plot(sqrt(x)/log(x), (2, B), color='purple')

︠da657687-cf1a-400a-83b0-17d06f08f43b︠
v3.plot()

v4 = stats.TimeSeries([v3[i]*log(i-1) for i in [3..len(v3)-1]])

v4.plot()

v4[-1]/20000

6.024322782230395e-05

N(sqrt(10000))

100.000000000000

5101648384
10^11

5101648384 100000000000

12/sqrt(10000)
︠fdb5d3de-f211-4405-b671-48fcecd644ed︠
len(v)

9999

v.plot(plot_points=10000)

list(v[:20])

[-1.0, -0.8815751854503011, -0.0775786850784419, -0.41057547008484097, 0.1770586104216081, -0.2881120141198865, 0.20855451944143866, 0.6760739727100677, 1.1204357246698047, 0.545845435603721, 0.9553838412940756, 0.35138416807063333, 0.7356618154418326, 1.109661078651051, 1.4745526835935667, 0.8313008199008909, 1.1807105326411023, 0.5234620527702756, 0.8601361975153274, 1.1912335155313478]

v[:20].mean()

0.4641806806761636

B = 20000
v = stats.TimeSeries([Li(n) - prime_pi(n) for n in [2..B]])
v2 = stats.TimeSeries([0]+[v[:i].mean() for i in range(1,len(v))])
v.plot(plot_points=B) + v2.plot(color='red', plot_points=B) + plot(sqrt(x)*log(x), (2, B), color='purple')

x=10^8
z = N(Li(x) - prime_pi(x))

z

753.330284251233

N(sqrt(x)*log(x))

184206.807439524

x

100000000