Calculate the density of prime numbers, using a couple of different methods, and compare them to "measured" values

prime_density.py

import pylab as plt
import numpy as np

def primes(x):
    p = []
    mask = range(2,x-1)
    
    while len(mask)>0:
        p.append(mask[0])
        del mask[::p[-1]]
        
    return p

def prime_density(x):
    p = primes(x)
    #print p
    Z = [(1 - 1./pi) for pi in p if pi < x/2.]
    z = reduce( lambda x,y: x*y, Z )
    
    rho = (len(primes(x+2000)) - len(primes(x-2000)) )/4e3
    
    print z,rho," ",x
    return (z,rho)

pd = []
y = np.arange(3.4,8,0.2)
x = map(int,10**y)
for xi in x:
    try: pd.append(prime_density(xi))
    except KeyboardInterrupt: break


z = np.array([pdi[0] for pdi in pd])
rho = np.array([pdi[1] for pdi in pd])
x = x[:len(z)]
f = 1./(0.5772 + np.log(x))

plt.figure()
plt.xlabel("N")
plt.ylabel("Density of primes")

plt.semilogx(x,rho,'bo',label="Measured")
plt.semilogx(x,f,'g',label="1/($\gamma$ + log x)")
plt.semilogx(x,z,'r',label="product formula")
plt.semilogx(x,(z+f)/2.,'b',label="average formula")
plt.legend(loc="upper right")
plt.show()