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*(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?

• Now I remember my thought: The condition that the difference has to be a pentagonal number as well might be satisfied much earlier for non consecutive pentagonal numbers than for consecutive. At least to me this is not quite obvious. – Nobody Jun 10 '15 at 16:37
• @Nobody This is correct. Otherwise it would be much easier. – Stefan Jun 10 '15 at 16:39
• i have one in java with 4 sec – RE60K Jun 10 '15 at 18:17

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.)

• Thank you! Implementing it straightforward (without factorization) and limiting x to the prior calculated limit, I already reach a time of 70ms. Though, the check if $P_j$ is pentagonal has to be added. – Stefan Jun 11 '15 at 18:32
• That's right, a check on $J$ is necessary too. Fixed. – Gareth Rees Jun 11 '15 at 19:05

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.