Plot_Pi

#coding: utf-8
from operator import mul
import console
console.clear()
print 'Generating plot... (this may take a little while)'

import numpy as np
import matplotlib.pyplot as plt

GAMMA = 0.57721566490153286061
li2 = 1.045163780117492784844588889194 # Li(x) = li(x) - li(2)

def ambi_sieve(n):
#http://tommih.blogspot.com/2009/04/fast-prime-number-generator.html
	s = np.arange(3, n, 2)
	for m in xrange(3, int(n ** 0.5)+1, 2): 
		if s[(m-3)/2]: 
			s[(m*m-3)/2::m]=0
	return np.r_[2, s[s>0]]
	
def factors(n):
	return (p for p in ambi_sieve(n + 1) if n % p == 0) #generator not list
	
def Mobius(n):
#Möbius function that takes the value 1 when n is 1, 0 when the factorization of n contains a duplicated factor, and (−1)**k, where k is the number of factors in the factorization of n, otherwise.
	if n == 1: 
		return 1
	else: 
		m = lambda f: 0 if (n / f) % f == 0 else -1 #zero if duplicated factor, -1 otherwise
		return reduce(mul, [m(k) for k in factors(n)], 1) #multiply all. if either is duplicated = 0

def li(x):
#Gauss gives an estimate π(x) ∼ Li(x), the logarithmic integral. Can be calculated by expanding the infinite series
#programmingpraxis.com/2011/07/29/approximating-pi
	return GAMMA + np.log(np.log(x)) + sum(pow(np.log(x), k) 
	                                       / (reduce(mul, xrange(1, k + 1)) * k) 
	                                       for k in xrange(1, 101))

def Li(x):
	return li(x) - li2

def R(x):
#In 1859, Bernhard Riemann described a different, and much more accurate, approximation to the prime counting function:
#mathworld.wolfram.com/RiemannPrimeCountingFunction.html
	return sum(Mobius(k) / k * li(pow(x, 1.0 / k)) for k in xrange(1, 101) if Mobius(k) != 0)
	
def Ri2(x):
#mathworld.wolfram.com/RiemannPrimeCountingFunction.html
#with tuple ignoring k where Mobius == 0
	mu = tuple(Mobius(k) for k in xrange(1,101))
	return sum(mu[k-1] / k * li(pow(x, 1.0 / k)) for k in xrange(1, 101) if mu[k-1] != 0)
	
moby = tuple(Mobius(k) for k in xrange(1,101)) #tuple of Mobius(k) computed once outside R(x) calls

def Ri4(x):
#mathworld.wolfram.com/RiemannPrimeCountingFunction.html
#with tuple ignoring k where Mobius == 0
	return sum(moby[k-1] / k * li(pow(x, 1.0 / k)) for k in xrange(1, 101) if moby[k-1] != 0)

plt.grid(False)
plt.title('n/log(n) ~ Li(n) ~ Ri(n)')
plt.yscale('log')
n = np.logspace(1, 3)
plt.ylim(10, 200)
p1 = plt.plot(n, n/np.log(n) , lw=2, c='r')
p2 = plt.plot(n, Li(n), lw=2, c='b')
p3 = plt.plot(n, Ri4(n), lw=2, c='g')
plt.legend([p1[0], p2[0], p3[0]], ['n/log(n)', 'Li(n)', 'Ri(n)'], loc=4)
plt.show()