# 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. Jun 10, 2015 at 16:37
• @Nobody This is correct. Otherwise it would be much easier. Jun 10, 2015 at 16:39
• i have one in java with 4 sec Jun 10, 2015 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. Jun 11, 2015 at 18:32
• That's right, a check on $J$ is necessary too. Fixed. Jun 11, 2015 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.