|
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
|
|
|
|
Loading history...
This check looks for lines that are too long. You can specify the maximum line length.