In [54]:
#--------------------------------------------------------------
#
# densityofprimes.py
#
# Paul Soper
#
# April 24, 2016
#
#--------------------------------------------------------------
In [55]:
import math
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [56]:
def isprime(q):
primeflag = 1;
if (q == 1):
return False
for x in range(2,int(math.sqrt(q)+1)):
if ((q % x) == 0):
return False
return True
In [57]:
vlog = np.vectorize(math.log)
In [58]:
p = []
upper_limit = 1000000
for x in range(2,upper_limit):
if isprime(x):
p.append(x)
In [59]:
num_bins = 50
bin_size = upper_limit/num_bins
n, bins, patches = plt.hist(p, num_bins, normed=0, facecolor='green', alpha=0.5)
t = np.linspace(2,upper_limit,num_bins)+(bin_size/2.0)
u = bin_size/vlog(t)
plt.plot(t,u,'r-',linewidth=3)