1
|
|
|
""" |
2
|
|
|
Project Euler Problem 27: Quadratic Primes |
3
|
|
|
========================================== |
4
|
|
|
|
5
|
|
|
.. module:: solutions.problem27 |
6
|
|
|
:synopsis: My solution to problem #27. |
7
|
|
|
|
8
|
|
|
The source code for this problem can be |
9
|
|
|
`found here <https://bitbucket.org/nekedome/project-euler/src/master/solutions/problem27.py>`_. |
10
|
|
|
|
11
|
|
|
Problem Statement |
12
|
|
|
################# |
13
|
|
|
|
14
|
|
|
Euler discovered the remarkable quadratic formula: |
15
|
|
|
|
16
|
|
|
.. math:: |
17
|
|
|
|
18
|
|
|
n^2 + n + 41 |
19
|
|
|
|
20
|
|
|
It turns out that the formula will produce :math:`40` primes for the consecutive integer values |
21
|
|
|
:math:`0 \\le n \\le 39`. However, when :math:`n = 40`, :math:`40^2 + 40 + 41 = 40(40 + 1) + 41` is divisible by |
|
|
|
|
22
|
|
|
:math:`41`, and certainly when :math:`n = 41`, :math:`41^2 + 41 + 41` is clearly divisible by :math:`41`. |
|
|
|
|
23
|
|
|
|
24
|
|
|
The incredible formula :math:`n^2 - 79 \\times n + 1601` was discovered, which produces :math:`80` primes for the |
|
|
|
|
25
|
|
|
consecutive values :math:`0 \\le n \\le 79`. The product of the coefficients, :math:`-79` and :math:`1601`, is |
|
|
|
|
26
|
|
|
:math:`-126479`. |
27
|
|
|
|
28
|
|
|
Considering quadratics of the form: |
29
|
|
|
|
30
|
|
|
.. math:: |
31
|
|
|
|
32
|
|
|
&n^2 + an + b, \\mbox{ where } |a| \\lt 1000 \\mbox{ and } |b| \\le 1000 \\\\ |
33
|
|
|
& \\\\ |
34
|
|
|
&\\mbox{where } |n| \\mbox{ is the modulus/absolute value of } n \\\\ |
35
|
|
|
&\\mbox{e.g. } |11| = 11 \\mbox{ and } |-4| = 4 |
36
|
|
|
|
37
|
|
|
Find the product of the coefficients, :math:`a` and :math:`b`, for the quadratic expression that produces the maximum |
|
|
|
|
38
|
|
|
number of primes for consecutive values of :math:`n`, starting with :math:`n = 0`. |
39
|
|
|
|
40
|
|
|
Solution Discussion |
41
|
|
|
################### |
42
|
|
|
|
43
|
|
|
If :math:`n^2 + an + b` is prime for :math:`0 \\le n \\le n^{\\prime}` where :math:`n^{\\prime} \\gt 40`, then we can |
|
|
|
|
44
|
|
|
assume that the solution maintains :math:`n^2 + an + b` as prime for the two cases :math:`n = 0, 1`. First, we'll |
|
|
|
|
45
|
|
|
consider these cases. |
46
|
|
|
|
47
|
|
|
Case (:math:`n=0`): |
48
|
|
|
|
49
|
|
|
.. math:: |
50
|
|
|
|
51
|
|
|
&n^2 + an + b = b \\\\ |
52
|
|
|
\\Rightarrow &b \\mbox{ must be a prime} |
53
|
|
|
|
54
|
|
|
Case (:math:`n=1`): |
55
|
|
|
|
56
|
|
|
.. math:: |
57
|
|
|
|
58
|
|
|
&n^2 + an + b = 1 + a + b \\\\ |
59
|
|
|
\\Rightarrow &1 + a + b \\mbox{ must be a prime} \\\\ |
60
|
|
|
\\Rightarrow &a = p - b - 1 \\mbox{ for some primes } b, p |
61
|
|
|
|
62
|
|
|
Now, this has set up a search space on :math:`a,b`. For each such :math:`a,b` we must determine the :math:`n^{\\prime}` |
|
|
|
|
63
|
|
|
s.t. :math:`n^2 + an + b` is prime for :math:`0 \\le n \\le n^{\\prime}`. |
64
|
|
|
|
65
|
|
|
The answer is simply the maximal value of :math:`n^{\\prime}` where :math:`|a| \\lt 1000 \\mbox{ and } |b| \\le 1000`. |
|
|
|
|
66
|
|
|
|
67
|
|
|
.. note:: the work in this algorithm is dominated by repeatably checking whether :math:`n^2 + an + b` is prime for |
|
|
|
|
68
|
|
|
various values of :math:`a,b,n`. Many tuples in this search will result in repeated values for |
|
|
|
|
69
|
|
|
:math:`n^2 + an + b` and so this algorithm benefits heavily from applying memoization to cache these results, |
|
|
|
|
70
|
|
|
avoiding redundant calculations. |
71
|
|
|
|
72
|
|
|
Solution Implementation |
73
|
|
|
####################### |
74
|
|
|
|
75
|
|
|
.. literalinclude:: ../../solutions/problem27.py |
76
|
|
|
:language: python |
77
|
|
|
:lines: 80- |
78
|
|
|
""" |
79
|
|
|
|
80
|
|
|
from itertools import product |
81
|
|
|
from math import ceil, log |
82
|
|
|
from typing import Callable |
83
|
|
|
|
84
|
|
|
from lib.numbertheory import is_probably_prime |
85
|
|
|
from lib.sequence import Primes |
86
|
|
|
from lib.util import memoize |
87
|
|
|
|
88
|
|
|
|
89
|
|
|
def find_n_prime(is_prime: Callable[[int], bool], a: int, b: int) -> int: |
|
|
|
|
90
|
|
|
""" Find :math:`n^{\\prime}` s.t. :math:`n^2 + an + b` is prime for :math:`0 \\le n \\le n^{\\prime}` |
|
|
|
|
91
|
|
|
|
92
|
|
|
:param is_prime: a primality testing function |
93
|
|
|
:param a: the parameter :math:`a` |
94
|
|
|
:param b: the parameter :math:`b` |
95
|
|
|
:return: :math:`n^{\\prime}` |
96
|
|
|
""" |
97
|
|
|
|
98
|
|
|
n = 1 # we have already tested n=0 since 0^2 + a*0 + b = b, and b is a prime |
|
|
|
|
99
|
|
|
while is_prime(n ** 2 + a * n + b): |
100
|
|
|
n += 1 |
|
|
|
|
101
|
|
|
return n |
102
|
|
|
|
103
|
|
|
|
104
|
|
|
def solve(): |
105
|
|
|
""" Compute the answer to Project Euler's problem #27 """ |
106
|
|
|
|
107
|
|
|
range_limit = 1000 |
108
|
|
|
|
109
|
|
|
maxsize = 2 ** int(ceil(log(2 * range_limit, 2))) # round 2 * range_limit up to a power of two |
110
|
|
|
is_prime = memoize(is_probably_prime, maxsize=maxsize) # memoization to avoid redundant calculations |
|
|
|
|
111
|
|
|
|
112
|
|
|
primes = list(Primes(upper_bound=range_limit)) # build a list of primes in the search space |
113
|
|
|
|
114
|
|
|
# Perform the search for the maximal n^{\\prime} given by an a, b in the search space |
115
|
|
|
max_n_prime = 0 |
116
|
|
|
for b, p in product(primes, repeat=2): |
|
|
|
|
117
|
|
|
a = p - b - 1 # a is determined by b, p |
|
|
|
|
118
|
|
|
new_n_prime = find_n_prime(is_prime, a, b) |
119
|
|
|
if new_n_prime > max_n_prime: |
120
|
|
|
max_n_prime = new_n_prime |
121
|
|
|
best = a, b |
122
|
|
|
|
123
|
|
|
# Calculate the answer (product of the coefficients) |
124
|
|
|
a, b = best |
|
|
|
|
125
|
|
|
answer = a * b |
126
|
|
|
|
127
|
|
|
return answer |
128
|
|
|
|
129
|
|
|
|
130
|
|
|
expected_answer = -59231 |
|
|
|
|
131
|
|
|
|
This check looks for lines that are too long. You can specify the maximum line length.