patternpythonMajor
Large Number Limit Extravaganza
Viewed 0 times
numberlargelimitextravaganza
Problem
I am writing a program that computes
$$n-n \cdot \left(\prod_{i=1}^n (1-\frac{1}{p_i})\right)$$
which is rewritten in my code as:
$$\left(1-\left(\prod_{i=1}^n(1-\frac{1}{p_i})\right)\right) \cdot n$$
There are lots of primes involved here especially as \$n\$ gets larger and larger. I do not want to use Atkin's method for generating them since my teacher and I will probably lose oversight as to what is going on, but would still like to optimize my algorithm for speed and accuracy (not sure if the latter can be better).
I am open to multithreading, and to translating this to a different language like C++, but do not know what optimizations I could make there either.
$$n-n \cdot \left(\prod_{i=1}^n (1-\frac{1}{p_i})\right)$$
which is rewritten in my code as:
$$\left(1-\left(\prod_{i=1}^n(1-\frac{1}{p_i})\right)\right) \cdot n$$
There are lots of primes involved here especially as \$n\$ gets larger and larger. I do not want to use Atkin's method for generating them since my teacher and I will probably lose oversight as to what is going on, but would still like to optimize my algorithm for speed and accuracy (not sure if the latter can be better).
i = 0
j = 3
old = 1.0
primes = [2]
while True:
if (j>240000): break
cont = True
enumerator = 0
for e in primes:
if j%e==0: cont = False
if enumerator>=len(primes)/2 + 1: break
enumerator += 1
if cont: primes.append(j)
j+=2
#print primes
while True:
if (i>len(primes)-1): break
old = old * (1 - 1.0/primes[i])
print old
i+=1;
primeslessthan = j-1
amta = 1-old
print "The convergence: " + str(primeslessthan)
print "All the way up to the number: " + str(amta)
print "The amount of active inhibitors are: " + str((1-amta)*primeslessthan)I am open to multithreading, and to translating this to a different language like C++, but do not know what optimizations I could make there either.
Solution
Your method of generating primes is horribly slow. Running your code takes on the order of 5 minutes on my machine.
Consider as an alternative a prime sieve, even a simple one like the Sieve of Eratosthenes will do (you don't need to go full Atkinson on it):
This runs in less than a second.
Accuracy
In your question you said you want to calculate $$n \cdot \left(1-\left(\prod_{i=1}^n(1-p_i)\right)\right)$$ But in your code you implemented $$n \cdot \left(1-\left(\prod_{i=1}^n(1-1/p_i)\right)\right)$$.
You do
Review:
Python has an official style-guide, which programmers are encouraged to follow, PEP8.
Your code violates it in the following points:
In addition, the following are discouraged for iterating in python:
Just use this simple syntax:
Note that the
Finally, the
Consider as an alternative a prime sieve, even a simple one like the Sieve of Eratosthenes will do (you don't need to go full Atkinson on it):
def prime_sieve(limit):
a = [True] * limit
a[0] = a[1] = False
for i, isprime in enumerate(a):
if isprime:
yield i
for n in xrange(i * i, limit, i):
a[n] = False
if __name__ == "__main__":
prod = 1
for p in prime_sieve(240000):
prod *= (1 - 1.0 / p)
print prod
primes_less_than = p - 1
amta = 1 - prod
print "The convergence:", primes_less_than
print "All the way up to the number:", amta
print "The amount of active inhibitors are:", prod * primes_less_thanThis runs in less than a second.
Accuracy
In your question you said you want to calculate $$n \cdot \left(1-\left(\prod_{i=1}^n(1-p_i)\right)\right)$$ But in your code you implemented $$n \cdot \left(1-\left(\prod_{i=1}^n(1-1/p_i)\right)\right)$$.
You do
amta = 1 -old and then later (1 - amta) * primeslessthan. Here you are needlessly loosing precision, because 1 - (1 - x) != x due to floating point precision. Better use old directly in the output.Review:
Python has an official style-guide, which programmers are encouraged to follow, PEP8.
Your code violates it in the following points:
- Always put the expression after an
ifin a new line
- Use spaces around operators
- Use
lower_casenames for variables and functions
oldis a bad name for a variable,amtaeven worse
In addition, the following are discouraged for iterating in python:
l = [1, 2, 3, ...]
for i in range(len(l)):
print l[i]
i = 0
while True:
if i > len(l) - 1:
break
print l[i]
i += 1Just use this simple syntax:
for x in l:
print xNote that the
if does not need any parenthesis.Finally, the
print statement can take multiple expressions, separated by commas. Alternatively, you can use str.format:print "The convergence: {}".format(primeslessthan)
print "All the way up to the number: {}".format(amta)
print "The amount of active inhibitors are: {}".format((1 - amta) * primeslessthan)Code Snippets
def prime_sieve(limit):
a = [True] * limit
a[0] = a[1] = False
for i, isprime in enumerate(a):
if isprime:
yield i
for n in xrange(i * i, limit, i):
a[n] = False
if __name__ == "__main__":
prod = 1
for p in prime_sieve(240000):
prod *= (1 - 1.0 / p)
print prod
primes_less_than = p - 1
amta = 1 - prod
print "The convergence:", primes_less_than
print "All the way up to the number:", amta
print "The amount of active inhibitors are:", prod * primes_less_thanl = [1, 2, 3, ...]
for i in range(len(l)):
print l[i]
i = 0
while True:
if i > len(l) - 1:
break
print l[i]
i += 1for x in l:
print xprint "The convergence: {}".format(primeslessthan)
print "All the way up to the number: {}".format(amta)
print "The amount of active inhibitors are: {}".format((1 - amta) * primeslessthan)Context
StackExchange Code Review Q#150442, answer score: 23
Revisions (0)
No revisions yet.