### Problem

Euler published the remarkable quadratic formula:

\[n^2+n+41\]

It turns out that the formula will produce 40 primes for the consecutive values $n=0$ to $39$. However, when $n=40$, $40^2+40+41=40(40+1)+41$ is divisible by 41, and certainly when $n=41$, $41^2+41+41$ is clearly divisible by 41.

Using computers, the incredible formula $n^2-79n+1601$ was discovered, which produces 80 primes for the consecutive values $n=0$ to $79$. The produce of the coefficients, $-79$ and $1601$ is $-126479$.

Considering quadratics of the form:

$n^2+an+b$ where $|a|

where $|n|$ is the modulus/absolute value of $n$, e.g., $|11|=11$ and $|-4|=4$.

Find the product of the coefficients $a$ and $b$, for the quadratic expression that produces the maximum number of primes for consecutive values of $n$, starting with $n=0$.

### Sage Solution

Some observations:

- Clearly, $b$ must be a prime or a negative prime, since $n=0$ must result in a prime.
- We can refine a lower bound on $b$ as follows: If a given combination $(a,b)$ results in $m$ consecutive primes, then clearly $b>m$, since otherwise we would obtain a factor of $m$ in the given polynomial. Also, when $a$ is negative, we must have $b>-(n^2+an)$, since the prime values $n^2+an+b$ must be positive. Since we know that $n^2+n+41$ returns 40 primes, we know that any interesting values for $a$ and $b$ must then satisfy $b>-(40^2+40a)$.

Given these observations, I’ll write a routine in Sage.

import time
start = time.time()
P = prime_range(1000)
a_max, b_max, c_max = 0, 0, 0
for a in range(-1000,1001):
for b in P:
if b < -1600-40*a or b < c_max: continue
c, n = 0, 0
while is_prime(n**2 + a*n + b):
c += 1
n += 1
if c > c_max:
a_max, b_max, c_max = a, b, c
prod = a_max * b_max
elapsed = time.time() - start
print "a=%s, b=%s, longest sequence = %s, prod=%s\nfound in %s seconds"\
% (a_max, b_max, c_max, prod, elapsed)

When executed, we obtain the following.

a=-61, b=971, longest sequence = 71, prod=-59231
found in 8.53578400612 seconds

It’s entirely possible that we’re being a bit redundant with the integer values that we’re throwing at Sage’s “is_prime” function. We’ll cache a list of primes within the range that we’re considering. (I came up with the number 751000 below in the Python section.)

import time
start = time.time()
P = prime_range(751000)
L = [False] * 751000
for p in P: L[p] = True
a_max, b_max, c_max = 0, 0, 0
for a in range(-1000,1001):
for b in range(1,1001):
if L[b] is False: continue
if b < -1600-40*a or b < c_max: continue
c, n = 0, 0
while L[n**2 + a*n + b] is True:
c += 1
n += 1
if c > c_max:
a_max, b_max, c_max = a, b, c
prod = a_max * b_max
elapsed = time.time() - start
print "a=%s, b=%s, longest sequence = %s, prod=%s\nfound in %s seconds"\
% (a_max, b_max, c_max, prod, elapsed)

This indeed executes more quickly, as is seen in the output.

a=-61, b=971, longest sequence = 71, prod=-59231
found in 1.39028286934 seconds

### Python Solution

Of course, Sage’s prime/factorization functionality isn’t available in Python, and so we’re going to need to create a prime sieve and store primes in some form. One of the things I need to know is how large the values $n^2+an+b$ can be. If we were to assume that the maximum value of the length of a prime sequence is 500 (which is an overestimate as we now know), then we’d have the maximum value of the quadratic given as:

500**2 + 1000 * 500 + 1000

Let’s create a prime sieve that will return a list of all primes below a given number.

def primes_xrange(stop):
primes = [True] * stop
primes[0], primes[1] = [False] * 2
L = []
for ind, val in enumerate(primes):
if val is True:
primes[ind*2::ind] = [False] * (((stop - 1)//ind) - 1)
L.append(ind)
return L

Let’s verify that it works.

[2, 3, 5, 7, 11, 13, 17, 19]

How long does this function take to form all of the primes under 751000?

time P = primes_xrange(751000)

Time: CPU 0.39 s, Wall: 0.39 s

Now, we’ll take the Sage routine from above and write it completely in Python. Here, I’m using only the sieve list portion of the primes_xrange function I just created. As above, I’m just checking to see if values in that list are True or False, corresponding to primes and composite numbers. This should take a bit longer to run than the Sage version above.

import time
def primes_xrange(stop):
primes = [True] * stop
primes[0], primes[1] = [False] * 2
for ind, val in enumerate(primes):
if val is True:
primes[ind*2::ind] = [False] * (((stop - 1)//ind) - 1)
return primes
start = time.time()
P = primes_xrange(751000)
a_max, b_max, c_max = 0, 0, 0
for a in range(-1000,1001):
for b in range(1,1001):
if P[b] is False: continue
if b < -1600-40*a or b < c_max: continue
c, n = 0, 0
while P[n**2 + a*n + b] is True:
c += 1
n += 1
if c > c_max:
a_max, b_max, c_max = a, b, c
prod = a_max * b_max
elapsed = time.time() - start
print "a=%s, b=%s, longest sequence = %s, prod=%s\nfound in %s seconds"\
% (a_max, b_max, c_max, prod, elapsed)

And it does run slightly slow compared to the Sage version, as expected.

a=-61, b=971, longest sequence = 71, prod=-59231
found in 1.73107814789 seconds

### Cython Solution

I’m going to take roughly the same approach as with the Python solution, but I’m going to malloc a C array instead of using the Python list. That should make things a bit more efficient. Before we do that, let’s see how much faster our routine runs if we just slap a “%cython” on the first line.

%cython
import time
def primes_xrange(stop):
primes = [True] * stop
primes[0], primes[1] = [False] * 2
for ind, val in enumerate(primes):
if val is True:
primes[ind*2::ind] = [False] * (((stop - 1)//ind) - 1)
return primes
start = time.time()
P = primes_xrange(751000)
a_max, b_max, c_max = 0, 0, 0
for a in range(-1000,1001):
for b in range(1,1001):
if P[b] is False: continue
if b < -1600-40*a or b < c_max: continue
c, n = 0, 0
while P[n**2 + a*n + b] is True:
c += 1
n += 1
if c > c_max:
a_max, b_max, c_max = a, b, c
prod = a_max * b_max
elapsed = time.time() - start
print "a=%s, b=%s, longest sequence = %s, prod=%s\nfound in %s seconds"\
% (a_max, b_max, c_max, prod, elapsed)

We find that this does improve the timing, and our Python code (which hasn’t been modified at all) now runs more quickly than the Sage code.

a=-61, b=971, longest sequence = 71, prod=-59231
found in 1.02033305168 seconds

Rewriting to take advantage of Cython’s C datatypes, how much faster can we obtain the result?

%cython
import time
from libc.stdlib cimport malloc, free
start = time.time()
cdef long i, j
cdef bint *primes = malloc(751000 * sizeof(bint))
primes[0], primes[1] = 0, 0
i = 2
while i < 751000:
primes[i] = 1
i += 1
i = 2
while i < 751000:
if primes[i] == 1:
j = 2 * i
while j < 751000:
primes[j] = 0
j += i
i += 1
cdef long long a, b, c, t
cdef long long pol, n
cdef long a_max, b_max, c_max, prod
a_max, b_max, c_max = 0, 0, 0
a = -1000
while a < 1001:
b = 2
while b < 1001:
if primes[b] == 0:
b += 1
continue
t = -1600-40*a
if b < t or b < c_max:
b += 1
continue
c, n = 0, 0
pol = n**2 + a*n + b
while primes[pol] == 1:
c += 1
n += 1
pol = n**2 + a*n + b
if c > c_max:
a_max, b_max, c_max = a, b, c
b += 1
a += 1
prod = a_max * b_max
free(primes)
elapsed = time.time() - start
print "a=%s, b=%s, longest sequence = %s, prod=%s\nfound in %s seconds"\
% (a_max, b_max, c_max, prod, elapsed)

Now, this should run much more quickly.

a=-61, b=971, longest sequence = 71, prod=-59231
found in 0.0381400585175 seconds

Thus, the Sage version runs 1.245 times as fast as the Python version. The Cython version runs 45.387 times as fast as the Python version and 36.452 times as fast as the Sage version.