Quadratic primes

Euler discovered the remarkable quadratic formula:

$$ n^2+n+41 $$ It turns out that the formula will produce 40 primes for the consecutive integer values $ 0 \le n \le 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.

The incredible formula $ n^2−79n+1601 $ was discovered, which produces 80 primes for the consecutive values $ 0 \le n \le 79 $. The product of the coefficients, −79 and 1601, is −126479.

Considering quadratics of the form:

$ n^2+an+b $, where |a|<1000 and |b|≤1000

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.

---

Idea

A problem of enumeration.

Key to this problem is to reduce search space.

Some observations:

• when n = 0, then b must be prime, so $ 2 \le b \le 1000 $
• when n = 1, then 1 + a + b must be prime, so $ 1-b \le a \lt 1000 $

Some more constraints might be applied:

• 1 + a + b be prime, and assume b != 2 and 1 + a + b != 2, we could constraint a to be odd, so could reduce search sapce by half

---

def is_prime(n):
    if n < 0:
        return False
    d = 2
    while pow(d, 2) <= n:
        if n % d == 0:
            return False
        d += 1
    return True

[(i, is_prime(i)) for i in range(1, 21)]

[(1, True),
 (2, True),
 (3, True),
 (4, False),
 (5, True),
 (6, False),
 (7, True),
 (8, False),
 (9, False),
 (10, False),
 (11, True),
 (12, False),
 (13, True),
 (14, False),
 (15, False),
 (16, False),
 (17, True),
 (18, False),
 (19, True),
 (20, False)]

def get_expression(a, b):
    return lambda n: pow(n, 2) + a * n + b

def solve(a_abs_bound, b_abs_bound):
    max_cnt = (-1, (0, 0, 0))
    b_candidates = filter(is_prime, range(2, b_abs_bound+1))
    for b in b_candidates:
        for a in range(1-b, a_abs_bound):
            e = get_expression(a, b)
            cnt = 0
            while is_prime(e(cnt)):
                cnt += 1
            max_cnt = max(max_cnt, (cnt, (a, b, a * b)))
    return max_cnt

solve(1000, 1000)

(71, (-61, 971, -59231))