Speedup for Project Euler 44 - pentagon numbers

Question

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*(3*n-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+24*x))) $ 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     -> difference
                       -- x       -> smaller number
                       -- x+pnl   -> larger number
                       -- 2*x+pnl -> sum
                       fpl = filter (\x -> ispn (x+pnl) && ispn (2*x+pnl)) npl

main = print finder

As you can see, the code goes through every pentagonal difference \$D=P_l\$ and checks all relevant numbers by calculating the upper limit for the smaller of both numbers (index nm; above, the distance between two pentagonal numbers is larger than \$D\$) and testing the resulting larger number and sum for their pentagonality.

Profiling with

ghc -O2 -fforce-recomp -rtsopts -prof -fprof-auto main.hs
./main +RTS -sstderr -p

yields an increase of runtime by 200% (from 34s to 103s) and reports the following costs:

COST CENTRE   MODULE  %time %alloc

finder'.fpl.\ Main     30.7    0.0
ispn          Main     23.7    0.0
ispn.n        Main     22.5    0.0
ispn.n.\      Main     12.6    0.0
finder'.fpl   Main     10.0    0.0
finder'.npl   Main      0.1   99.9

I already tried to convert the check for pentagonality to a set lookup, but this didn't really help because this will grow up to 2 million elements. Also, converting the to-be-filtered list to an array or set didn't work as expected.

For a given difference \$D\$ (given by the pentagonal number \$P_l\$ which is pnl = pn l in the code), we know that starting with a given \$k\$ (nm in the code) the difference between \$P_k\$ and \$P_j\$ (\$j > k\$) will always be larger than \$D\$. This starting index can be calculated as \$k = \lfloor\frac{D-1}{3}\rfloor + 1\$.

Is there any way to squeeze some additional time out of this code, preferably while sticking to base libraries?

Answer 1

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:

1. Iterate over all values for \$k\$ below the bound that you already established.
2. 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\$.
3. 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.)

Answer 2

I found an additional approach by browsing the solver forum of Project Euler. This approach is inspired by the solution of the user "observ" (post #68 for those, who have access), which attempts a proof of the solution, but stops before doing this.

The approach is to iterate over the larger number's index \$k\$. During this calculation, the smallest occurring difference is stored (initialized with \$\infty\$). Upto the point, where the first hit is found, all indices \$j < k\$ for the smaller number have to be checked.

After the first hit, the range of \$j\$ is limited by the distance \$D > P_k - P_j\$. The calculation finally ends when \$D < P_k - P_{k-1}\$.

My implementation for this takes around 400ms in Haskell and 200ms in Fortran. Accordingly, the algorithm of Gareth Rees is preferred.