Pentagon Numbers

It's not difficult to test if a number is pentagonal.

from math import sqrt

def is_pentagonal(x):
    n = (sqrt(24*x+1) + 1) / 6
    return n == int(n)

With this function, and given a value $j$, we can attempt to find a value $k$ such that $P_j + P_k$ and $P_j - P_k$ are both pentagonal by just iterating over all $k < j$.

def pentagonal_pair(j):
    a = polygonal_number(5, j)
    for k in reversed(range(1, j)):
        b = polygonal_number(5, k)
        if is_pentagonal(a + b) and is_pentagonal(a - b):
            return k
        
    return None

A naive approach that happens to work is to just iterate over values of $j$ until this function returns a solution.

from itertools import count

for j in count(2):
    k = pentagonal_pair(j)
    if k is not None:
        break

j, k

(2167, 1020)

We can then recover the pentagonal numbers from these indices.

a, b = polygonal_number(5, j), polygonal_number(5, k)
a, b

(7042750, 1560090)

And calculate their difference.

a - b

5482660

This turns out to be the answer! The thing is... we sort of just got lucky that the first $j,k$ pair we found happened to minimize $|P_j - P_k|$. If you care to find out how we can verify that this is the minimal solution, read on.

Diophantine approach

Let's restate the problem. We are looking for positive integers $j, k, s, d$ such that $$P_j + P_k = P_s$$ $$P_j - P_k = P_d$$

Suppose $d$ is some constant. Then we can solve for $j$ and $k$ in the second equation using Diophantine methods. SageMath can do this for us with solve_diophantine, but we can do it more efficiently ourselves.

$P_j - P_k = P_d$ expands to $$3j^2 - j - 3k^2 + k = 3d^2 - d$$ We can factor the left hand side into $$(j-k)(3j + 3k - 1) = 3d^2 - d$$ This means that any solution $j,k$ corresponds to a factorization of $3d^2 - d$. Importantly, this means there's a finite number of solutions for a fixed $d$.

So if we factor $3d^2 - d = ab$, and set $a=j-k$ and $b=3j+3k-1$, we can solve for $j$ and $k$: $$k = \frac{b+1-3a}{6}$$ $$j = a + k$$ If we iterate over all possible factorizations $ab$ of $3d^2 - d$, this will let us find all $j,k$ that solve $P_j - P_k = P_d$.

Once we have found candidate pairs for $j$ and $k$, we can just check if $P_j + P_k$ is also a pentagonal number using the function from above - if so, $j$ and $k$ also satisfy $P_j + P_k = P_s$, and we've found a solution for our fixed value of $d$.

Here's a function that applies all of this to find solutions (if they exist) for a given $d$.

def pentagonal_sums(d):
    psums = set()
    N = 3*d^2 - d
    for a in divisors(N):
        b = N / a
        
        k = (b + 1 - 3*a) / 6
        if k != int(k) or k <= 0:
            continue
            
        j = a + k
            
        p = polygonal_number(5, j) + polygonal_number(5, k)
        if is_pentagonal(p):
            psums.add((j, k))
    
    return psums

Since we want to minimize $P_d$, we can just iterate over increasing values of $d$ and stop once this function finds a solution. This way, we know our solution is minimal since we didn't find solutions for any smaller $d$.

for d in count(1):
    sols = pentagonal_sums(d)
    if sols != set():
        break
            
sols

{(2167, 1020)}

This confirms that the solution we found earlier is in fact the solution that minimizes $P_d$.

Relevant sequences

• Pentagonal numbers: A000326
• Pentagonal numbers which are the sum of two other positive pentagonal numbers: A136117

Copyright (C) 2025 filifa

This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International license and the BSD Zero Clause license.