Project Euler Problem 27: Quadratic Primes¶
The source code for this problem can be found here.
Problem Statement¶
Euler discovered the remarkable quadratic formula:
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:
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\)):
Case (\(n=1\)):
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}\)
- is_prime (
-
solutions.problem27.
solve
()¶ Compute the answer to Project Euler’s problem #27