Project Euler Solution 44: Pentagon numbers

Today in the Project Euler series we are going to look at Problem 44: Pentagon numbers where one has to find a pair of pentagonal numbers such that both their sum and difference is also pentagonal.

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

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$?

This is one of the problems where one has to search a two dimensional lattice, spanned by $j$ and $k$. This lattice is infinite, so we cannot search all of it sensibly. If it had only one dimension, we could stop as soon as we had the solution. But with two dimensions, we have to think of a sensible search path through this lattice.

We want two pentagonal numbers such that their sum and difference are also pentagonal. Once we have a solution, we will have four pentagonal numbers. Without loss of generality I will have $k > j$ and therefore we can write the sum as $P_k + P_j = P_s$ and $P_k - P_j = p_d$. We need to find $k$ and $j$ such that $s$ and $d$ exist and are integers. But as the problem asks us to find $k$ and $j$ such that $d$ is minimized, we could also just start to loop over $d$ and then pick either $k$, $j$ or $s$ for the second loop. We generate two pentagonal numbers and then have to check whether the other two numbers are also pentagonal.

Let's start writing a check which verifies whether something is a pentagonal number. This can be done by inverting the condition, just like we did for Solution 42: Coded triangle numbers. A number $y$ is pentagonal if it can solve the equation $$ \frac{n(3n-1)}{2} = y $$ which $n$ being an integer. We can apply the abc-formula and solve this equation as $$ n = \frac{1 \pm \sqrt{1 + 24 y}}{6} \,. $$ We don't care about the negative branch because $n$ must be a positive number. So a number $y$ is a pentagonal number if $1 + 24 y$ is a perfect square and also if $1 + \sqrt{1 + 24 y}$ is divisible by 6.

We can implement that:

def is_pentagonal(number: int) -> bool:
    s = int(round(math.sqrt(1 + 24 * number)))
    return s**2 == 1 + 24 * number and (1 + s) % 6 == 0

Then it is handy to have an iterator for pentagonal numbers, also easy to write as a generator. Here I use that this is a quadratic formula such that the distance between the numbers grows linearly by 3. This gives the following code:

def iter_pentagonal_numbers() -> Iterator[int]:
    p = 1
    d = 4
    while True:
        yield p
        p += d
        d += 3

Combining these we can write a little test to make sure that the check for pentagonal numbers work:

def test_is_pentagonal() -> None:
    assert is_pentagonal(1)
    assert is_pentagonal(5)
    assert is_pentagonal(12)
    assert is_pentagonal(22)
    assert not is_pentagonal(2)
    assert not is_pentagonal(3)
    assert not is_pentagonal(4)

    for p in itertools.islice(iter_pentagonal_numbers(), 1_000_000):
        assert is_pentagonal(p)

So it does what we expect for a few hand-picked cases and also gets the first million pentagonal numbers classified correctly without rounding issues. Good!

My attempt then was to parametrize $P_d$ first and then iterate over the lower numbers $P_j$. Then I can construct the upper number $K = P_d + P_j$. I need to make sure whether $K$ is a pentagonal number and could be written as $P_k$. Then I take the sum $S = P_j + K$ and check whether this can be written as $P_s$.

I do need some criterion to stop with the $P_j$ and move to the next $P_d$. This is reached when $P_d < P_{j+1} - P_j$ because then $P_j + P_d$ cannot be a pentagonal number. In the code it looks like this:

def solution_diff() -> int:
    for p_diff in iter_pentagonal_numbers():
        for p_lower, p_lower_next in itertools.pairwise(iter_pentagonal_numbers()):
            upper = p_lower + p_diff
            if upper < p_lower_next:
                break
            sum_ = p_lower + upper
            if is_pentagonal(upper) and is_pentagonal(sum_):
                return p_diff

I am pretty sure that it would work, it just doesn't produce a solution in a sensible time. I couldn't really think of a different sensible way to slice the problem. Because if I parametrize $k$ and $j$, I would not know when to stop. And if I parametrized $d$ and $s$, then I wouldn't know when to move on to the next $d$. Parametrizing with $s$ and $d$ would give an easy stopping for $d$ as $P_d < P_s$ must hold. But it wouldn't guarantee that the first element that I find is the one with minimal $d$. I was stuck.

Different way

Then I looked at somebody else's solution and found that it is pretty much that. It goes over the sum $P_s$ and then over $P_j$. This means that the $P_j$ start off small and the $K$ are pretty large. This makes the difference $D$ pretty large as well. Still, the first $P_d$ which is found is actually the solution of the problem.

def solution_sum_lower() -> int:
    for p_sum in iter_pentagonal_numbers():
        for p_lower in iter_pentagonal_numbers():
            upper = p_sum - p_lower
            if upper < p_lower:
                break
            diff = upper - p_lower
            if is_pentagonal(upper) and is_pentagonal(diff):
                return diff

This produces the correct answer of 5,482,660 within 901 ms. I am a bit flabbergasted that this has worked, although I didn't think that this would come up with the correct solution.