I have been solving some random Project Euler problems to practice my haskell. After solving the problem, I usually look up the solution on the haskell wiki.
For Problem 27, I solved it the regular way, i.e, using a combination of isPrime
and map
s. But then, I saw this solution on the wiki. I have no idea how this solution works. The only other mention I find of it is from this closed thread on math stackexchange.
problem_27 :: Integer
problem_27 = -(2 * a - 1) * (a ^ 2 - a + 41)
where
n = 1000
m = head $ filter (\x -> x ^ 2 - x + 41 > n) [1 ..]
a = m - 1
As I understand, this code does the following
But I do not understand why this program gives the correct answer. Or in other words, What is the reasoning behind this algorithm?
Notice that
(n-x)^2 + a(n-x) + b = n^2 + (a-2x)n + x^2 - ax + b
which means that a change of variables can turn any (integer-coefficient) quadratic into one of the forms n^2+d or n^2+n+d with a judicious choice of x. Now the first form is no good: every other term is even, so we can't have a long run of primes. For quadratics of the second form, Rabinowitz proved that the quadratic that produces the maximum sequence of primes must have 4d-1 a Heegner number, of which there are only 6. It's quick from there to check that n^2+n+41 is the best of those six, and the next best, n^2+n+17, produces fewer consecutive primes total (i.e. both positive and negative inputs) than n^2+n+41 does starting at 0.
Now, (-1-n)^2+(-1-n)+41 = n^2+n+41, so if n produces a prime, then -1-n does, too. This means about half the consecutive primes produced by this quadratic are produced by a negative input, so the more of those we can shift into range without blowing our coefficient budget, the better; and, what's more, we can't do better than trying to shift in more values, because shifting any other polynomial around gets fewer consecutive primes than before we even shift this one.
What is our budget? Well, looking back at our very first equation, now with a=1 and b=41, we must have:
|1-2x| < 1000
|x^2-x+41| < 1000
because these are the final two coefficients. The second is going to be our limiter, since x^2 generally grows much faster than 2x, so we are simply looking for the largest x that validates that equation. This is what the
m = head $ filter (\x -> x ^ 2 - x + 41 > n) [1 ..]
a = m-1
snippet does. Once we have that, 1-2x and x^2-x+41 are our two coefficients, and we want to compute their product, which is what the
-(2 * a - 1) * (a ^ 2 - a + 41)
snippet does.