Large Number Limit Extravaganza

Question

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 = 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.

• Phrancis: Welcome to Code Review! I have rolled back the last edit. Please see what you may and may not do after receiving answers .
• Graipher: In your code you seem to calculate n*(1 - Pi(1 - 1/p)), instead of n*(1 - Pi(1 - p)).

Answer 1

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

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_than

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 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:

1. Always put the expression after an if in a new line
2. Use spaces around operators
3. Use lower_case names for variables and functions
4. old is a bad name for a variable, amta even 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 += 1

Just use this simple syntax:

for x in l:
 print x

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

• Linus Rastegar: Thanks! What is the "[True]" syntax you were using in line 2?
• Graipher: @LinusRastegar [True] is a list with one element and that element is the value True. In addition, [True] * 3 == [True, True, True]. This also works e.g. with strings, "a" * 3 == "aaa".
• Linus Rastegar: Thanks for this clarification! I was able to run the code and of course everything runs much faster now! I am now trying to test with limit = 1 billion. Is there some multithreading trick I could use to make the program run faster?
• Graipher: Probably not with a sieve, because you need to mark off all multiples of prime numbers off as being composite, this won't work in parallel. One thing you could do is write the primes to a file. For some bound this might be faster than calculating them all, because increasing the limit gives not linear increase in runtime.
• Graipher: @LinusRastegar For all primes up to 1,000,000,000, you probably want to comment out the print, though.