Experiments_with_Mobius

from operator import mul
import itertools
import numpy as np
import datetime

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

def factors(n):
#with list
	return [p for p in ambi_sieve(n + 1) if n % p == 0]
	
def factors2(n):
#with generator
	return (p for p in ambi_sieve(n + 1) if n % p == 0) #generator

def factors3(n):
#with yield
	for p in ambi_sieve(n + 1):
		if n % p == 0:
			yield p	#yield

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 all_factors(prime_dict):
	series = [[p**e for e in range(maxe+1)] for p, maxe in prime_dict.items()]
	for multipliers in itertools.product(*series):
		yield reduce(mul, multipliers)

def divisors(factors):
#Generates all divisors, unordered, from the prime factorization.
	ps = sorted(set(factors))
	omega = len(ps)
	def rec_gen(n = 0) :
		if n == omega :
			yield 1
		else :
			pows = [1]
			for j in xrange(factors.count(ps[n])) :
				pows += [pows[-1] * ps[n]]
			for q in rec_gen(n + 1) :
				for p in pows :
					yield p * q

	for p in rec_gen() :
		yield p

def mu(n):
# Works. But difficult to understand code for n == 1
	m = lambda p: -1 if (n / p) % p else 0
	return reduce(mul, [m(p) for p in factors3(n)], 1)
	
def mobius(n):
#My definition of Mobius.
	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 m(n, f):
#lambda from above in explicit way
	if (n/f)%f == 0:
		return 0
	else:
		return -1
		
def li(x):
#programmingpraxis.com/2011/07/29/approximating-pi
#mathworld.wolfram.com/LogarithmicIntegral.html
	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):
#mathworld.wolfram.com/RiemannPrimeCountingFunction.html
#suboptimal. don't bother with k where mobius(k)==0
	return sum(mobius(k) / k * li(pow(x, 1.0 / k)) for k in xrange(1, 101))
	
def Ri(x):
#mathworld.wolfram.com/RiemannPrimeCountingFunction.html
#with generator computing Mobius(k) twice but ignoring k where Mobius == 0
	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
	mob = tuple(mobius(k) for k in xrange(1,101))
	return sum(mob[k-1] / k * li(pow(x, 1.0 / k)) for k in xrange(1, 101) if mob[k-1] != 0)
	
def Ri3(x):
#mathworld.wolfram.com/RiemannPrimeCountingFunction.html
#with generator ignoring k where Mobius == 0
	mob = ((k, mobius(k)) for k in xrange(1,101))
	return sum(mu / k * li(pow(x, 1.0 / k)) for (k, mu) in mob if mu != 0)

moby = tuple(mobius(k) for k in xrange(1,101)) #tuple of Mobius(k) computed 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)


N = 1000000
#print ambi_sieve(n), factors(n), divisors(factors(n)), mobius3(n)

#print divisors([2, 3, 5])

#print sorted(all_factors({2:3, 3:2, 5:1}))
		
#print [m(N, k) for k in factors(N)], mobius3(N)

#print [mobius3(n) for n in xrange(100)]

#print [mu(n) for n in xrange(N)]

#mob = [mobius3(k) for k in xrange(1, N)]
#print mob, len(mob)

tc=datetime.datetime.now()
print Ri4(N), datetime.datetime.now() - tc
tc=datetime.datetime.now()
print Ri2(N), datetime.datetime.now() - tc
tc=datetime.datetime.now()
print Ri(N), datetime.datetime.now() - tc
tc=datetime.datetime.now()
print R(N), datetime.datetime.now() - tc

#print reduce(mul, [m(n, k) for k in factors(n)], 1)

#print [mobius2(i) for i in xrange(1000)]

#print [factors(i) for i in xrange(1000)]

#for N in range(10, 10000, 1000):
#	print R(N), li(N), Li(N), N/np.log(N)

#print factors(N), factors2(N), factors3(N)