12
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ Welcome to Code Review! I have rolled back the last edit. Please see what you may and may not do after receiving answers. \$\endgroup\$ – Phrancis Dec 21 '16 at 2:41
  • \$\begingroup\$ In your code you seem to calculate n*(1 - Pi(1 - 1/p)), instead of n*(1 - Pi(1 - p)). \$\endgroup\$ – Graipher Dec 21 '16 at 8:01
16
\$\begingroup\$

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)
\$\endgroup\$
  • \$\begingroup\$ Thanks! What is the "[True]" syntax you were using in line 2? \$\endgroup\$ – Linus Rastegar Dec 21 '16 at 2:31
  • \$\begingroup\$ @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". \$\endgroup\$ – Graipher Dec 21 '16 at 2:33
  • \$\begingroup\$ 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? \$\endgroup\$ – Linus Rastegar Dec 21 '16 at 2:38
  • \$\begingroup\$ 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. \$\endgroup\$ – Graipher Dec 21 '16 at 2:40
  • 2
    \$\begingroup\$ @LinusRastegar For all primes up to 1,000,000,000, you probably want to comment out the print, though. \$\endgroup\$ – Graipher Dec 21 '16 at 9:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.