Project Euler Problem 27: Quadratic Primes

The source code for this problem can be found here.

Problem Statement

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 - 79 \times n + 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:

\[\begin{split}&n^2 + an + b, \mbox{ where } |a| \lt 1000 \mbox{ and } |b| \le 1000 \\ & \\ &\mbox{where } |n| \mbox{ is the modulus/absolute value of } n \\ &\mbox{e.g. } |11| = 11 \mbox{ and } |-4| = 4\end{split}\]

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\).

Solution Discussion

If \(n^2 + an + b\) is prime for \(0 \le n \le n^{\prime}\) where \(n^{\prime} \gt 40\), then we can assume that the solution maintains \(n^2 + an + b\) as prime for the two cases \(n = 0, 1\). First, we’ll consider these cases.

Case (\(n=0\)):

\[\begin{split}&n^2 + an + b = b \\ \Rightarrow &b \mbox{ must be a prime}\end{split}\]

Case (\(n=1\)):

\[\begin{split}&n^2 + an + b = 1 + a + b \\ \Rightarrow &1 + a + b \mbox{ must be a prime} \\ \Rightarrow &a = p - b - 1 \mbox{ for some primes } b, p\end{split}\]

Now, this has set up a search space on \(a,b\). For each such \(a,b\) we must determine the \(n^{\prime}\) s.t. \(n^2 + an + b\) is prime for \(0 \le n \le n^{\prime}\).

The answer is simply the maximal value of \(n^{\prime}\) where \(|a| \lt 1000 \mbox{ and } |b| \le 1000\).

Note

the work in this algorithm is dominated by repeatably checking whether \(n^2 + an + b\) is prime for various values of \(a,b,n\). Many tuples in this search will result in repeated values for \(n^2 + an + b\) and so this algorithm benefits heavily from applying memoization to cache these results, avoiding redundant calculations.

Solution Implementation

from itertools import product
from math import ceil, log
from typing import Callable

from lib.numbertheory import is_probably_prime
from lib.sequence import Primes
from lib.util import memoize


def find_n_prime(is_prime: Callable[[int], bool], a: int, b: int) -> int:
    """ Find :math:`n^{\\prime}` s.t. :math:`n^2 + an + b` is prime for :math:`0 \\le n \\le n^{\\prime}`

    :param is_prime: a primality testing function
    :param a: the parameter :math:`a`
    :param b: the parameter :math:`b`
    :return: :math:`n^{\\prime}`
    """

    n = 1  # we have already tested n=0 since 0^2 + a*0 + b = b, and b is a prime
    while is_prime(n ** 2 + a * n + b):
        n += 1
    return n


def solve():
    """ Compute the answer to Project Euler's problem #27 """

    range_limit = 1000

    maxsize = 2 ** int(ceil(log(2 * range_limit, 2)))  # round 2 * range_limit up to a power of two
    is_prime = memoize(is_probably_prime, maxsize=maxsize)  # memoization to avoid redundant calculations

    primes = list(Primes(upper_bound=range_limit))  # build a list of primes in the search space

    # Perform the search for the maximal n^{\\prime} given by an a, b in the search space
    max_n_prime = 0
    for b, p in product(primes, repeat=2):
        a = p - b - 1  # a is determined by b, p
        new_n_prime = find_n_prime(is_prime, a, b)
        if new_n_prime > max_n_prime:
            max_n_prime = new_n_prime
            best = a, b

    # Calculate the answer (product of the coefficients)
    a, b = best
    answer = a * b

    return answer


expected_answer = -59231

solutions.problem27.find_n_prime(is_prime, a, b)

Find \(n^{\prime}\) s.t. \(n^2 + an + b\) is prime for \(0 \le n \le n^{\prime}\)

Parameters:

• is_prime (Callable[[int], bool]) – a primality testing function
• a (int) – the parameter \(a\)
• b (int) – the parameter \(b\)

Return type:

int

Returns:

\(n^{\prime}\)

solutions.problem27.solve()

Compute the answer to Project Euler’s problem #27