Prime less than N that is a sum of the most consecutive other primes

Question

Here's a Python implementation of the Sieve of Eratosthenes:

def return_primes(upto=100):
    primes = []
    sieve = set()
    for i in range(2, upto+1):
        if i not in sieve:
            primes.append(i)
            sieve.update(range(i, upto+1, i))
    return primes

I use the above function in the below algorithm to find the prime that is the sum of the most consecutive primes. I'm sure it has room for improvement. It starts with a small sliding window which iteratively increases in size, and each iteration goes until the found sum is greater than the largest prime, and then breaks to the outer iteration. Optimized for Python 2:

def find_prime_sum_of_most_consecutive_primes_under(x):
    '''
    given x, an integer, return the prime that is the sum
    of the greatest number of consecutive primes under x
    '''
    primes = return_primes(x)
    length = len(primes)
    primeset = set(primes)
    maxprime = primes[-1] # last prime
    result = None
    for window in xrange(1, length): # xrange should be range in Python 3
        iters = length - window
        for it in xrange(iters):
            s = sum(primes[it:window+it])
            if s in primeset:
                result = s
            if s > maxprime:
                break
    return result

Can this be improved?

Answer 1

Multiples of 2 cannot be prime and primes greater than 2 are odd. This means that the sum to be found must be the sum of an odd number of primes if 2 is not in the window. This lets you skip half the work.

You can update a sliding sum by subtracting the outgoing number and adding the incoming one. Perhaps this could be the basis of a window that moves forward as it expands and contracts, much like a caterpillar moves... The idea is to avoid redoing lots of sums without having to store lots of intermediary results. The basic sliding sum should already save quite a bit of work, though.

This can give you a quick exploratory run with sequences (sums) of maximum possible length, if you expand two sums from 2 and 5 respectively - skipping every other prime - and start sliding once the sums exceed the maximum prime. The disadvantage is that you do not gain information for bounding future searches, as in the bottom-up case where you never have to consider sequences shorter than the longest you've found so far. Perhaps it is better to caterpillar right from the start, to quickly collect growing lower bounds on sequence length.

Growing pair of snakes: slide two sums forward, offset by one prime, that are one double step longer than the longest sequence found so far. Expand as you find longer and longer sequences. That should give you some decent results in one stereo slither through the range, which can help bound future searches. (hindsight: one snake, single steps; it's only the growing that can be double-stepped)

UPDATE

I did a bit of snake charming and it didn't pan out - snakes grow very slowly and they have to move around a lot in order to do so. By and large there is little that can be done in the way of improving/refining bounds for the search, apart from using the length of the best sequence found so far (plus two) as a lower bound for future searches.

A lot of re-computation can be avoided by computing cumulative sums for the primes and storing them in a table; that way the sum from prime i to prime k is available as p_sums[k] - p_sums[i - 1] for i greater than zero.

Also, searches go faster if all possible sequence lengths are tried for a given starting point before moving to the next starting point, rather than fixing the sequence length for a given round and then trying all starting points. However, I haven't studied this point thoroughly; the success of this strategy may depend on the fact that a lot of winning sequences happen to begin at very low primes for the small range I tested (2³²).

'Brute force with superior firepower' is several orders of magnitude faster than the plain brute-force algorithm:

*** g++ 4.8.1 64 ***

v0:          10       0.000 ms          5 = {2} @ 2
v1:          10       0.000 ms          5 = {2} @ 2
v0:         100       0.001 ms         41 = {6} @ 2
v1:         100       0.000 ms         41 = {6} @ 2
v0:        1000       0.010 ms        953 = {21} @ 7
v1:        1000       0.001 ms        953 = {21} @ 7
v0:       10000       0.122 ms       9521 = {65} @ 3
v1:       10000       0.002 ms       9521 = {65} @ 3
v0:      100000       1.505 ms      92951 = {183} @ 3
v1:      100000       0.013 ms      92951 = {183} @ 3
v0:     1000000      19.793 ms     997651 = {543} @ 7
v1:     1000000       0.104 ms     997651 = {543} @ 7
v0:    10000000     254.493 ms    9951191 = {1587} @ 5
v1:    10000000       0.782 ms    9951191 = {1587} @ 5
v0:   100000000    4485.791 ms   99819619 = {4685} @ 7
v1:   100000000       6.226 ms   99819619 = {4685} @ 7
v0:  1000000000   68065.223 ms  999715711 = {13935} @ 11
v1:  1000000000      52.905 ms  999715711 = {13935} @ 11
v0: 10000000000 1359343.949 ms 9999419621 = {41708} @ 2
v1: 10000000000     485.861 ms 9999419621 = {41708} @ 2

The last round used more than 7 GB of core for the two tables, which means this is pretty much the end of the line for this simple algorithm even though the actual search took only half a second. The time for sieving the primes - about one second per 2^31 numbers for my simple code - is not included.

I'm showing the routine as C++ since my Python-fu is not sufficient for doing either python or the algorithm justice. At least that way you have an accurate description for the algorithm, and can ponder about devious ways of having Python do all the hard work:

typedef uint32_t num_t;
typedef uint64_t sum_t;

std::vector<num_t> primes;
std::vector<sum_t> p_sums;    // sum(primes[0] .. primes[i]) for i < primes.size()

void prepare_prime_tables_for_limit (num_t n);
bool is_prime (num_t n);      // uses std::lower_bound() on primes[], i.e. binary search

num_t brute_force_v1 (unsigned n = 1000000000)
{
   assert( n >= 10 );

   prepare_prime_tables_for_limit(n);

   num_t prime_count = num_t(primes.size());
   num_t max_prime = primes.back();
   num_t best_len = 1;
   num_t best_idx = 0;
   sum_t best_sum = 0;

   // separate pass for the prime that is so odd because it's even

   for (num_t i = 1; ; i += 2)
   {           
      sum_t sum = p_sums[i];

      if (sum <= max_prime)
      {
         if (is_prime(num_t(sum)))
         {
            best_len = i + 1;
            best_sum = sum;
         }
      }
      else break;
   }

   --best_len;  // adjust for the evenness of 2

   for (num_t i = 1; i + best_len + 1 < prime_count; ++i)
   {
      for (num_t k = i + best_len + 1; ; k += 2)
      {
         sum_t sum = p_sums[k] - p_sums[i - 1];

         if (sum <= max_prime)  
         {
            if (is_prime(num_t(sum)))
            {
               best_len = k - i + 1;
               best_idx = i;
               best_sum = sum;
            }
         }
         else break;
      }
   }

   best_len += best_idx == 0;  // restore the stolen bit if 2 is still best

   // ... print stuff ...

   return num_t(best_sum);
}

I left index math etc. unoptimised in order to make the code easier to understand. For the same reasons I didn't harden it by adding the complications needed for handling boundary cases (n < 8 etc.). My intent was to make the basic algorithm as clear as possible.

Sums need a bigger type than indexes and numbers, hence the occasional cast to num_t when it is safe (and necessary) to do so.

Answer 2

Your code is pretty good. For example, I didn't know that searching membership on a set() is faster than doing the same on a list().

Minor updates from me:

1. Function name is too long.
2. Variables names can be improved. Use of iters and it confuses with pythonic concept of iterators.
3. The most important update is that we can return from the function when, for a given window, the first offset has a sum > maxprime. We will not find a smaller sum in higher windows.
4. Minor optimization with the use of elif statement.

The following code is written in Python3:

def max_cumulative_prime(x):
    '''
    given x, an integer, return the prime that is the sum
    of the greatest number of consecutive primes under x
    '''
    primes = return_primes(x)
    length = len(primes)
    primeset = set(primes)
    maxprime = primes[-1] # last prime
    result = None
    for window in range(1, length):
        for offset in range(length - window):
            s = sum(primes[offset:offset+window])
            if s in primeset:
                result = s
            elif s > maxprime:
                if offset==0: return result
                else: break
    return result

Answer 3

Here's a modification that makes it run in seconds as opposed to minutes. I add a minimum length parameter for the window so that I don't waste time checking fruitless short runs:

def return_primes(upto=100):
    primes = []
    sieve = set()
    for i in xrange(2, upto+1):
        if i not in sieve:
            primes.append(i)
            sieve.update(range(i, upto+1, i))
    return primes

def find_prime_sum_of_most_consecutive_primes_under(x):
    '''
    given x, an integer, return the prime that is the sum
    of the greatest number of consecutive primes under x
    '''
    primes = return_primes(x)
    length = len(primes)
    primeset = set(primes)
    maxprime = primes[-1] # last prime
    result = None
    min_length = 1
    for start in xrange(0, length):
        for window_length in xrange(start + min_length, length - start + 1):
            check_primes = primes[start:window_length]
            s = sum(check_primes)
            if s in primeset:
                min_length = len(check_primes)
                result = s
            elif s > maxprime:
                break
    return result