Speedup for Project Euler 44 - pentagon numbers

I'm trying to solve the problems of Project Euler with Haskell within less than 10s with optimized code (ghc -O2). Unfortunately, I'm struggling to find a correct algorithm for Problem 44, which proves the correctness of the solution without setting a predefined upper limit (setting a limit and omitting the proof yields a fraction of a second... unfortunately, you'll almost only find this solution while searching the web, which might be due to a slightly different wording of the question in the past).


Pentagonal numbers are generated by the formula, \$P_n=n(3n−1)/2\$. The first ten pentagonal numbers are:


\$1, 5, 12, 22, 35, 51, 70, 92, 117, 145, \ldots\$


It can be seen that \$P_4 + P_7 = 22 + 70 = 92 = P_8\$. However, their difference, \$70 − 22 = 48\$, is not pentagonal.


Find the pair of pentagonal numbers, \$P_j\$ and \$P_k\$, for which their sum and difference are pentagonal and \$D = |P_k − P_j|\$ is minimised; what is the value of \$D\$?

The way I want to go is straightforward. Check every possible difference D, till a solution is found. Unfortunately, this means, that a total of around \$1.2\cdot10^{9}\$ checks for pentagonal numbers have to be made.

My solution is the following:

``
-- calculate the n-th pentagonal number
pn :: Int -> Int
pn n = n(3n-1) quot 2

-- check if number is pentagonal
ispn :: Int -> Bool
ispn x = x == pn n
where n = round $ (\x -> 1/6 (1 + sqrt(1+24x))) $ fromIntegral x

-- solution finder
-- l > difference index
-- m > prior limit of list with relevant pentagonal numbers
-- pl > list of pentagonal numbers
finder :: Int
finder = finder' 1 0 []
finder' :: Int -> Int -> [Int] -> Int
finder' l m pl = if length fpl /= 0
then l
else finder' (l+1) nm npl
where pnl = pn l
nm = (pnl - 1) div` 3 + 1 -- +1 as backup
npl = (map pn [(m+1)..nm]) ++ pl
-- pnl -> differenc

You're right that the problem condition that \$\left|P_k - P_j\right|\$ be minimized means that to properly solve it, it's not sufficient to find a pair of pentagonal numbers whose sum and difference are both pentagonal: you have to also establish that there is no pair of pentagonal numbers with the required property and a smaller difference.

Once you've found one pair of pentagonal numbers with the required property, this gives you a bound on \$j\$ and \$k\$; unfortunately the bound is in the millions, so that any \$\Omega(n^2)\$ algorithm (in particular, any algorithm that compares pairs of pentagonal numbers below the bound) will take far too long. It's no good just speeding up what you've got: you have to take a completely new approach.

So here's a rough sketch of an alternative algorithm.

The \$n\$th pentagonal number is $$P_n = {n(3n-1)\over 2}.$$ By completing the square, $$P_n = {(6n - 1)^2 - 1 \over 24}.$$ If \$P_k - P_j = P_x\$ and \$P_k + P_j = P_y\$ are both pentagonal numbers, then $$ \eqalign{ (6k - 1)^2 - (6j - 1)^2 &= (6x - 1)^2 - 1 \\ (6k - 1)^2 + (6j - 1)^2 &= (6y - 1)^2 + 1 } $$ Write \$K = 6k - 1\$, \$X = 6x - 1\$, and \$Y = 6y - 1\$. Add the two equations above, getting $$ 2K^2 = X^2 + Y^2. $$

So the plan is:

• Iterate over all values for \$k\$ below the bound that you already established.

• For each \$k\$, let \$K = 6k - 1\$, and enumerate all the ways to write \$2K^2\$ as the sum of two squares \$X^2 + Y^2\$.

• For each such decomposition, if \$X ≡ -1 \pmod 6\$ and \$Y ≡ -1 \pmod 6\$ and \$J = \sqrt{K^2 - X^2 + 1} ≡ -1 \pmod 6\$, then this is a solution.

In step 2, the enumeration of \$2K^2\$ as sums of two squares can be computed efficiently by factorizing \$K\$ and then using the Brahmagupta–Fibonacci identity. Since you are going to be needing a lot of these factorizations, it will be most efficient to sieve for them.

This ought to be doable within your 10 second goal.

(Credit: I adapted this strategy from Daniel Fischer's post in the Project Euler forum.)

StackExchange Code Review Q#93232, answer score: 7