Detailed mathematical approach

Problem Summary

We seek primes \(p \lt 10^6\) such that there exists a positive integer \(n\) with \(n^3+n^2p\) equal to a perfect cube. Writing that cube as \(m^3\), the equation becomes

$m^3-n^3=pn^2.$

A naive search over pairs \((n,m)\) would be pointless, because the interesting part of the problem is structural: the prime factor \(p\) forces the difference of the two cubes into a very narrow form. Once that form is identified, the problem collapses to testing a short explicit sequence for primality.

Mathematical Approach

The implementations rely on an exact characterization of all possible primes \(p\), not on a broad search. The whole argument comes from factoring the cube equation carefully and using the fact that \(p\) is prime.

Reduce the equation to a coprime core

Assume

$m^3=n^3+pn^2,$

with \(m,n \gt 0\). Since the right-hand side is larger than \(n^3\), we must have \(m \gt n\).

Let

$d=\gcd(m,n),\qquad m=du,\qquad n=dv,\qquad \gcd(u,v)=1.$

Substituting into the equation gives

$d^3u^3-d^3v^3=pd^2v^2,$

so after dividing by \(d^2\),

$d(u^3-v^3)=pv^2.$

Now \(\gcd(u,v)=1\) implies \(\gcd(u^3-v^3,v)=1\), because any common divisor of \(u^3-v^3\) and \(v\) would also divide \(u\). Therefore the factor \(v^2\) on the right must come entirely from \(d\). In particular,

$v^2 \mid d.$

Why the prime forces a single cube difference

Write

$d=v^2c.$

Then the previous equation becomes

$c(u^3-v^3)=p.$

Both factors on the left are positive integers. Because \(p\) is prime, there are only two possibilities: either \(c=1\) and \(u^3-v^3=p\), or \(c=p\) and \(u^3-v^3=1\).

The second case cannot happen for positive \(v\). If \(v \ge 1\), then the smallest possible positive cube difference is

$ (v+1)^3-v^3 = 3v^2+3v+1 \ge 7.$

So we must have

$c=1,\qquad p=u^3-v^3,\qquad d=v^2.$

This already tells us something strong: every valid prime is literally a difference of two cubes, and the associated \(n\) is

$n=dv=v^3.$

Only consecutive cubes can produce a prime

Factor the difference of cubes:

$p=u^3-v^3=(u-v)(u^2+uv+v^2).$

Because \(m \gt n\), we have \(u \gt v \ge 1\). The second factor is therefore greater than 1. Since \(p\) is prime, the first factor must be 1:

$u-v=1.$

Set \(v=k\), so \(u=k+1\). Then

$p=(k+1)^3-k^3=3k^2+3k+1.$

Conversely, whenever \(3k^2+3k+1\) is prime, choosing

$n=k^3,\qquad m=k^2(k+1)$

gives

$n^3+n^2p=k^9+k^6(3k^2+3k+1)=k^6(k+1)^3=(k^2(k+1))^3=m^3.$

So the condition is exact: the admissible primes are precisely the prime numbers of the form

$p=3k^2+3k+1,\qquad k \ge 1.$

These are often called Cuban primes of the first kind.

Worked Example

Take \(k=2\). Then

$p=3(2)^2+3(2)+1=19,$

and the construction gives

$n=2^3=8,\qquad m=2^2\cdot 3=12.$

Checking directly,

$8^3+8^2\cdot 19=512+64\cdot 19=1728=12^3.$

The small sample below 100 follows immediately from the same formula:

$k=1,2,3,4 \Longrightarrow p=7,19,37,61,$

while \(k=5\) gives \(p=91\), which is composite. That is why the sample count is 4.

Bound the search range

To solve the actual problem, we only need values with

$3k^2+3k+1 \lt 10^6.$

The largest such integer is \(k=576\), because

$3\cdot 576^2+3\cdot 576+1=997057 \lt 10^6,$

whereas

$3\cdot 577^2+3\cdot 577+1=1000519 \gt 10^6.$

So there are only 576 candidates to test.

How the Code Works

The C++, Python, and Java implementations do not search over arbitrary values of \(n\) or \(m\). They use the exact parameterization above and examine only the numbers that can possibly work.

Generate the only possible candidates

The implementation starts with \(k=1\), computes

$p=3k^2+3k+1,$

and stops as soon as \(p\) reaches the upper limit. Because the derivation is an equivalence, this loop misses nothing and produces no duplicates: each loop value corresponds to one possible prime and to the unique witness \(n=k^3\).

Test primality by trial division

Each candidate \(p\) is then checked for primality. The C++ and Java implementations handle the even case separately and test only odd divisors up to \(\sqrt{p}\). The Python implementation uses straightforward trial division from small divisors upward until the square-root bound is reached. Since there are only 576 candidates, this simple approach is entirely sufficient. The C++ version also contains a small checkpoint confirming that the count below 100 is 4.

Complexity Analysis

There are exactly 576 candidate values below \(10^6\). With trial division, each primality test costs \(O(\sqrt{p})\), so the total running time is

$O\!\left(\sum_{k=1}^{576}\sqrt{3k^2+3k+1}\right),$

which is comfortably tiny for this input size. In simpler terms, the whole job is just a few hundred square-root-bounded divisibility checks. Memory usage is \(O(1)\), because the implementations keep only a handful of integers and a running count.

Footnotes and References

1. Project Euler problem page: https://projecteuler.net/problem=131
2. Difference of two cubes: Wikipedia - Difference of two cubes
3. Greatest common divisor: Wikipedia - Greatest common divisor
4. Prime number: Wikipedia - Prime number
5. Trial division: Wikipedia - Trial division
6. Cuban prime: Wikipedia - Cuban prime