Project Euler 50 in Python

Question

Project Euler, problem #50:

The prime 41, can be written as the sum of six consecutive primes:

41 = 2 + 3 + 5 + 7 + 11 + 13 This is the longest sum of consecutive primes that adds to a prime below one-hundred.

The longest sum of consecutive primes below one-thousand that adds to a prime, contains 21 terms, and is equal to 953.

Which prime, below one-million, can be written as the sum of the most consecutive primes?

I came up with this code, but it only works decently for primes below ten thousand.

Any ideas on how to optimize it?

from pyprimes import *

def sequence_exists(l,ls,limit = 100):
    for x in range(0,len(ls)-l):
        if x+l > len(ls): return False
        if any (ls[i] > limit/6 for i in range(x,x+l,1)) :
            return False
        test_sum = sum(ls[x:x+l:1])
        if (test_sum <limit) and is_prime(test_sum) :
            return True
    return False

def main():
    n = prime_count(10000)
    prime_list = list(nprimes(n))
    l = 6
    for x in range(6,len(prime_list)):
        if sequence_exists(x,prime_list,10000):
            l=x
    print l

if __name__ == '__main__':
    main()

Answer 1

Don't import 'wildcards'

Especially if the module is not part of the Standard Library people will be confused and not understand where the function is coming from. You can cut down typing as follows:

 import pyprimes as pr

Use longer names

def sequence_exists(l,ls,limit = 100):

In the above l and ls tell nothing to the reader.

Answer 2

Any ideas on how to optimize it?

You are calculating sums too many times. Sub-sequences of primes are also recalculated multiple times without reusing previous calculations.

I think using prefix sums can help here. With prefix sums you can get the sum of any sub-sequence in \$O(1)\$ time instead of \$O(N)\$ which should help.

Redundant code

This won't make a difference in the speed of your solution, but the if condition here is redundant and you can safely delete it, the range already includes that condition:

for x in range(0,len(list_of_primes)-length):
    if x+length > len(list_of_primes): return False

By the way, when a range starts from 0, you can omit the 0.

Redundant range steps

In multiple places, you use the step parameter for ranges = 1 when it's unnecessary, as the default is 1 anyway, for example here:

    if any (list_of_primes[i] > limit/6 for i in range(x,x+length,1)) :

    test_sum = sum(list_of_primes[x:x+length:1])

You could write these slightly simpler as:

    if any (list_of_primes[i] > limit/6 for i in range(x,x+length)) :

    test_sum = sum(list_of_primes[x:x+length])

Naming

@Caridorc already pointed this out, but I have to re-state it too because it's so important. The logic of this program is really not trivial, and having descriptive variable names is crucial for understanding the code.

Follow PEP8

PEP8 is the recommended coding style guide for Python, and your formatting violates it at multiple places. Give it a good read and follow it in the future.

Answer 3

Any ideas on how to optimize it?

I'm not going to include a lot of detail, because this is a project Euler problem and I don't want to give too much away. However, here are some ideas which might help you along.

• You're looking for the largest number which specifies a certain condition. Your program runs an iterative loop over a set of numbers to test for the condition. Where does your loop start? How could you change it to find the answer faster? What assumptions can you make about the data and the problem statement to focus your search?
• You use an is_prime function to test a computed summation for primality. You don't show us the implementation, but this is guaranteed to be nonoptimal. At a bare minimum, you could inline it and eliminate the overhead of a function call (and you're making this call a lot). But, overhead aside, there's almost certainly a much better way to test for primality. Rather than performing e.g. Miller Rabin on every candidate, is there some simple, in-place check you could make instead, to see if your candidate is part of a known group of primes? How could you structure such a check efficiently?
• There are a lot of unnecessary recomputations. The biggest offender is covered in janos' answer, but there are some small ones as well, such as the length of the primes list and limit/6.
• Prime generation: you're counting the number of primes below 10000 and then generating that many primes (hopefully in sequence). The earlier Euler problems should have led you to develop an efficient way to sift out primes below a given number. You should be doing that. If you haven't done the earlier problems, I encourage you to go try them; the later problems build upon knowledge gained by completing earlier ones.

For reference, my python 3 answer to this problem finishes in under half a second on a poor quality laptop, and doesn't use precomputed sums. The substantial differences between our algorithms are the ones hinted at above.

Style

• Naming. Like the other two answers say, this is critical. Your code is confusing and hard to follow. The local variable names are often not meaningful.
• Spacing and line breaks aren't consistent. Pick a convention and stick to it. Preferably, add spaces between comparison operators and always include line breaks after conditionals. It helps readability a lot. I personally also really like spaces between function arguments. Example of inconsistency:
  • if x+l > len(ls): return False
  • if (test_sum <limit) and is_prime(test_sum) : return True
• When making a range or slice, if you're starting at 0 or using a step size of 1, you don't need to declare it. You even do this a little already (another example of inconsistency):
  • for i in range(x,x+l,1)
  • for x in range(6,len(prime_list)):
• Magic numbers. What is the meaning of 10000? Why is it that value and not something else? Declare a variable at the start with a descriptive name and set it instead.

Answer 4

1. Review

1. Using from module import * is not normally a good idea, because you don't know all the names that are being imported. This might accidentally import a name that you want to use for something else, or worse, it might work fine for now but break after an update when the module adds some new exported names. It's more reliable just to import the names you need:

   from pyprimes import is_prime, nprimes, prime_count
   

2. Instead of:

   n = prime_count(10000)
   prime_list = list(nprimes(n))
   

   Use the end keyword argument to primes, like this:

   prime_list = list(primes(end=10000))
   

3. The number 10000 should be a parameter, to avoid repetition and to make it convenient to test the program.

4. The function sequence_exists is only called from one place in the code, so little is gained by making it a separate function. In fact, doing so makes it harder to spot performance optimizations. Better to inline it at the single point of use.

5. range(0, end) can be simplified to range(end) since ranges start at 0 by default. Similarly, range(start, end, 1) can be simplified to range(start, end) since ranges have step 1 by default.

6. Project Euler problem 50 asks you "Which prime, below one-million, can be written as the sum of the most consecutive primes?" (my emphasis), so you need to return the sum, not the number of consecutive primes.

7. The second line here is unnecessary:

   for x in range(0,len(ls)-l):
       if x+l > len(ls): return False
   

   because the loop ensures that x is always smaller than len(ls)-l, and so the condition will always be false.

8. The test:

   if any (ls[i] > limit/6 for i in range(x,x+l,1)) :
       return False
   

   is wrong! For suppose that we have x = 0, l = 6 and limit = 42. Then ls[4] = 11 which is greater than limit/6 = 7 and so we'll execute return False. But in fact sum(ls[0:6]) is 2 + 3 + 5 + 7 + 11 + 13 = 41, which is prime and less than limit.

   The test that you need here is:

   if ls[x] > limit / l:
       return False
   

9. Starting the search with sequences of length 6 means that the code can't be tested on very small instances. For example, if we reduce the limit from 10000 to 10, main still prints 6, even though the longest consecutive prime sum in this case is 2 + 3 = 5.

   I understand why you've done this: it's because you know from the example in the Project Euler problem statement (2 + 3 + 5 + 7 + 11 + 13 = 41) that when the limit is greater than 41, the longest consecutive prime sum must have length at least 6, so you can save time by not trying sums of length 1 to 5.

   But by implementing this optimization, you've made your code give incorrect results for small limits, and that makes it hard to check for correctness. It would be better to make it correct in all circumstances, and then figure out a general way to avoid doing unnecessary work (see §2.4 below).

10. Python ranges are exclusive at the top, so this loop will never consider summing all the primes in the list:

    for x in range(6,len(prime_list)):
    

    You might say that we never need to sum all the primes in the list, because they can't sum to another prime in the list, but that's false if the list only has one prime! Correctness on small examples is important for testing, so the loop needs to be:

    for x in range(1, len(prime_list) + 1):
    

11. Similarly, this loop never considers summing the last l primes in the list:

    for x in range(0,len(ls)-l):
    

    It needs to be:

    for x in range(0, len(ls) - l + 1):
    

2. Performance

Timing the original code, using limit=40000 as our test case:

>>> from timeit import timeit
>>> timeit(lambda:main(40000), number=1)
38.903770233970135

1. Applying all the changes in §1 above results in this code:

   from pyprimes import is_prime, primes
   
   def longest_prime_sum(limit):
       """Return the prime below limit that can be written as as the sum of
       the most consecutive primes.
   
       >>> longest_prime_sum(100)
       41
   
       """
       prime_list = list(primes(end=limit))
       result = None
       for l in range(1, len(prime_list) + 1):
           for x in range(len(prime_list) - l + 1):
               if prime_list[x] > limit / l:
                   break
               test_sum = sum(prime_list[x:x+l])
               if test_sum < limit and is_prime(test_sum):
                   result = test_sum
       return result
   

   This is already more than twenty times faster (mainly due to §1.8):

   >>> timeit(lambda:longest_prime_sum(40000), number=1)
   1.983643549028784
   

2. If you look at the implementation of is_prime, you'll see that it uses trial division, followed by a Miller–Rabin primality test if trial division fails. But that's wasteful in this case, because we already know all the primes below limit as we just made a list of them. So instead of calling is_prime(test_sum), we can make a set of primes:

   prime_set = set(prime_list)
   

   (so that membership testing is efficient) and then:

   if test_sum < limit and test_sum in prime_set:
   

   This is now about 25 times faster than the original:

   >>> timeit(lambda:longest_prime_sum(40000), number=1)
   1.5118930450407788
   

3. Instead of exiting the loop when prime_list[x] > limit / l, it would be better to exit when test_sum >= limit. This is a tighter condition, so exits earlier. The code is now more than 100 times faster than the original:

   >>> timeit(lambda:longest_prime_sum(40000), number=1)
   0.35101655102334917
   

4. Since the goal is to produce the longest sum, why do we start with small sums and work our way up? If we started with big sums and worked our way down, then as soon as we found a sum we could stop, knowing immediately that it was the longest.

   What's the longest sum we might have to consider? Well, we can work it out by adding up primes 2 + 3 + 5 + … until they reach limit. No sum can be longer than that and still meet the conditions of the problem.

   Here's the revised code:

   def longest_prime_sum(limit):
       """Return the prime below limit that can be written as as the sum of
       the most consecutive primes.
   
       >>> longest_prime_sum(100)
       41
   
       """
       prime_list = list(primes(end=limit))
       prime_set = set(prime_list)
   
       # Find max number of consecutive primes whose sum is below limit.
       total = 0
       for max_length, p in enumerate(prime_list):
           total += p
           if total >= limit:
               break
   
       for length in range(max_length + 1, 0, -1):
           for x in range(len(prime_list) - length + 1):
               total = sum(prime_list[x:x+length])
               if total >= limit:
                   break
               elif total in prime_set:
                   return total
   

   This is now about 1500 times as fast as the original:

   >>> timeit(lambda:longest_prime_sum(40000), number=1)
   0.021531195961870253
   

   And solves the Project Euler problem in less than a second:

   >>> timeit(lambda:longest_prime_sum(1000000), number=1)
   0.54454751603771
   

5. The runtime is now dominated by the call to pyprimes.primes: this takes more than 90% of the time, so that's the place to put effort. Using sieve5 from this answer instead of pyprimes speeds things up by another four times, and using sieve8 (which requires NumPy) gives a 5× speedup on top of that, solving the Project Euler problem in 27 milliseconds:

   >>> timeit(lambda:longest_prime_sum(1000000), number=1)
   0.026997262961231172