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)
DensityOfPrimes