Problem 131: Prime Cube Partnership

Project Euler Problem 131 Solution

EulerSolve provides an optimized solution for Project Euler Problem 131, Prime Cube Partnership, with C++, Python, Java, and a step-by-step mathematical explanation.

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

Detailed mathematical approach

Problem Summary

We look for primes $p$ for which some positive integer $n$ makes $n^3+n^2p$ a perfect cube, that is, equal to $m^3$ for a positive integer $m$. For $p=19$ the statement gives $8^3+8^2\cdot19=12^3$. It adds that for every prime with this property the value of $n$ is unique, and that there are only four such primes below $100$. We must count the primes below one million that have the property. ^[1]

The whole solution in one sentence: for a prime $p$, the number $n^3+n^2p$ is a cube exactly when $n=k^3$ and $p=(k+1)^3-k^3=3k^2+3k+1$ for some $k\ge1$, so we only have to test the $576$ numbers $3k^2+3k+1$ below one million for primality.

Mathematical Approach

Trying $n=1,2,3,\dots$ finds a solution when one exists, but it can never show that a prime has none. Instead, a small search suggests a pattern, and prime factorisations prove that it describes every solution.

1. A search for small primes

For each of the $25$ primes below $100$ we tried every $n\le10^5$. Only four primes have a solution in this range, and each has only one:

$p=7$: $1^3+1^2\cdot7=8=2^3$;
$p=19$: $8^3+8^2\cdot19=512+1216$${}=1728=12^3$;
$p=37$: $27^3+27^2\cdot37=19683+26973$${}=46656=36^3$;
$p=61$: $64^3+64^2\cdot61=262144+249856$${}=512000=80^3$.

Three patterns stand out. The values $n=1,8,27,64$ are the cubes $k^3$ for $k=1,2,3,4$. The ratios $m/n=2,\frac32,\frac43,\frac54$ are $(k+1)/k$, so $m=k^2(k+1)$. And each prime is the difference of two consecutive cubes: $7=8-1$, $19=27-8$, $37=64-27$ and $61=125-64$. The next sections prove that these patterns hold for every solution.

2. Cubes and prime factorisations

By the fundamental theorem of arithmetic, every integer greater than $1$ is a product of primes, unique apart from the order of the factors. ^[2] Cubing a number triples every exponent, and because the factorisation is unique, a positive integer is a cube exactly when every exponent in its factorisation is a multiple of $3$; it is then the cube of the number with a third of each exponent. For example $1728=2^6\cdot3^3$ is the cube of $2^2\cdot3=12$, while $686=2\cdot7^3$ is not a cube, because the prime $2$ has the exponent $1$.

Lemma. If two positive integers $a$ and $b$ have no common prime factor and $ab$ is a cube, then $a$ and $b$ are cubes.

Proof. Every prime of $ab$ occurs in exactly one of the two factors, so its exponent in that factor equals its exponent in $ab$, a multiple of $3$. Hence all exponents of $a$ and of $b$ are multiples of $3$.

The condition on common factors matters: $2\cdot4=8$ is a cube, but neither $2$ nor $4$ is.

3. The common divisors of $n$ and $n+p$

The expression factors as $n^3+n^2p=n^2(n+p)$, so we need $n^2(n+p)=m^3$. A number that divides $n$ and $n+p$ also divides their difference $p$, and a number that divides $n$ and $p$ also divides $n+p$. So $n$ and $n+p$ have the same common divisors as $n$ and $p$; this is the observation on which the Euclidean algorithm rests. ^[3] Since $p$ is prime, there are two cases: either $p$ does not divide $n$, and then $n$ and $n+p$ have no common prime factor, or $p$ divides $n$.

4. When $p$ does not divide $n$

The primes of $n^2$ are those of $n$, so $n^2$ and $n+p$ have no common prime factor either, and the lemma makes both of them cubes. Every exponent of $n^2$ is twice an exponent $e$ of $n$, and $2e$ is a multiple of $3$ only when $e$ is, so $n$ itself is a cube: $n=k^3$ with $k\ge1$. Write $n+p=j^3$; then $j\gt k$, because $p\gt0$. The difference of two cubes factors as ^[4]

$$p=j^3-k^3=(j-k)\,(j^2+jk+k^2).$$

Both factors are positive integers, and the second is at least $2^2+2\cdot1+1^2=7$, since $j\ge2$ and $k\ge1$. A prime is not a product of two integers greater than $1$, so the first factor is $1$. Hence $j=k+1$ and

$$p=(k+1)^3-k^3=3k^2+3k+1,\qquad n=k^3.$$

Finally $m^3=n^2(n+p)=k^6(k+1)^3$, so $m=k^2(k+1)$: these are the three patterns of Section 1.

5. When $p$ divides $n$

Write $n=pt$ with $t\ge1$. Then

$$n^2(n+p)=p^2t^2\cdot p\,(t+1)=p^3\,t^2(t+1).$$

Dividing by $p^3$ lowers the exponent of $p$ by $3$ and keeps all other exponents, so $t^2(t+1)$ would have to be a cube. A common divisor of $t$ and $t+1$ divides their difference $1$, so the argument of Section 4 applies again: $t=s^3$ and $t+1=r^3$ for positive integers $s$ and $r$. But $s^3\lt s^3+1\lt(s+1)^3$ for $s\ge1$, so $s^3+1$ lies strictly between two consecutive cubes and is not a cube. This case never occurs. For instance, $p=n=7$ gives $7^2\cdot14=686=2\cdot7^3$, the non-cube of Section 2.

6. The characterisation and the uniqueness of $n$

Conversely, if $p=3k^2+3k+1$ is prime for some $k\ge1$, then $n=k^3$ works, because $k^3+3k^2+3k+1=(k+1)^3$:

$$\begin{aligned}n^3+n^2p&=k^9+k^6\,(3k^2+3k+1)\\&=k^6\,(k+1)^3=\bigl(k^2(k+1)\bigr)^3.\end{aligned}$$

So a prime has the property exactly when it equals $3k^2+3k+1$ for some $k\ge1$. Sections 4 and 5 also show that every solution is $n=k^3$ with this $k$, and since $3k^2+3k+1$ grows strictly with $k$, the prime determines $k$ and therefore $n$. This proves the uniqueness that the statement calls surprising. Primes that are differences of two consecutive cubes are called cuban primes of the first kind; the first of them are $7,19,37,61,127,271,331$. ^[5] ^[6]

7. Worked example: the primes below $100$

For $k=1,2,3,4,5$ the candidates are $7,19,37,61,91$, and $k=6$ already gives $127\gt100$. The first four are prime, while $91=7\cdot13$. So exactly four primes below $100$ have the property, as the statement says. For $p=19$ we have $k=2$, so $n=2^3=8$ and $m=2^2\cdot3=12$, the example of the statement. In the same way there are $11$ such primes below $10^3$, $28$ below $10^4$ and $64$ below $10^5$.

8. The candidates below one million

The candidates $3k^2+3k+1$ increase with $k$, so they lie below $10^6$ exactly for $k=1,2,\dots,576$:

$$\begin{aligned}3\cdot576^2+3\cdot576+1&=997057\lt10^6,\\3\cdot577^2+3\cdot577+1&=1000519\gt10^6.\end{aligned}$$

Each candidate is tested by trial division: a number $c\gt1$ without a divisor $d$ with $1\lt d\le\sqrt c$ is prime, because a factorisation $c=de$ with $1\lt d\le e$ would give $d^2\le c$. ^[7] Since $\sqrt{997057}\lt999$, divisors up to $998$ suffice for every candidate. Every candidate is odd, because $k$ or $k+1$ is even and so $3k(k+1)$ is even.

How the Code Works

All three programs apply Section 6 directly: they never compute $n$ or $m$ and only test the candidates of Section 8. The C++ program ^[8] keeps the default limit $=10^6$ and the switch run_checkpoints in the structure Options. main first calls parse_arguments, which accepts --skip-checkpoints and --limit= followed by decimal digits; parse_int_after_prefix collects these digits in an int without checking for overflow. Any other argument, or a limit below $2$, ends the program with exit code $1$. Unless the checkpoints are skipped, run_checkpoints requires solve(100) to return $4$, the count from the statement, and otherwise ends the program with exit code $2$. Finally main prints solve(limit).

In solve the loop variable a plays the role of $k$. For $a=1,2,3,\dots$ the loop computes p = 3 * a * a + 3 * a + 1, leaves at the first value that is not below the limit, which is correct because the candidates increase, and counts the values for which is_prime returns true. is_prime rejects numbers below $2$, accepts $2$ as the only even prime and then divides by $3,5,7,\dots$ as long as the square of the divisor does not exceed the number; it forms this square in std::int64_t, so the square cannot overflow. The even case never occurs for a candidate (Section 8). For the default limit the loop computes $577$ values and tests $576$ of them. With the usual 32-bit int, every limit up to $2147409811$, the candidate for $a=26754$, is safe; a larger limit overflows in 3 * a * a at $a=26755$ or, above $2^{31}-1$, already in the parser. That is undefined behaviour in C++, but it does not affect the default limit or the checkpoint.

The Python version is one function solve with the limit $10^6$ built in, without options or checkpoint. Its loop is the same, but its primality test tries every $d=2,3,4,\dots$ while d*d <= p, even divisors included, and its guard p > 1 is always true for a candidate; Python integers have unlimited precision. The Java version has no solve method: it contains the test as isPrime and the same loop directly in main, with the fixed limit and without the checkpoint. Like C++, it keeps the candidates in 32-bit int variables, which is safe at this limit, and forms the divisor square as (long) p * p. All three programs print the count.

Complexity and Verification

Below a limit $L$ there are about $K=\sqrt{L/3}$ candidates. The $k$-th is about $3k^2$, so its trial division tests at most about $\sqrt3\,k$ divisors, half as many in C++ and Java, which skip even divisors. Since $1+2+\dots+K\approx K^2/2$, the total is at most about $\sqrt3\,K^2/2=L/(2\sqrt3)$, roughly $0.29\,L$ divisions: the time is $O(L)$, and the memory is a fixed number of integer variables. For $L=10^6$ the full ranges would take $143586$ divisions in C++ and Java and $287462$ in Python, but most composite candidates stop much earlier, at their smallest prime factor, and the main pass actually performs $51767$ and $103622$ divisions, respectively. The C++ checkpoint adds $9$, and C++ and Java also make one parity test per candidate. On a desktop computer each program finishes in well under a tenth of a second. A sieve of Eratosthenes up to $L$ would need $O(L\log\log L)$ operations and $O(L)$ memory, ^[9] more than the few hundred tests need here.

The C++ program checks the statement's count of four primes below $100$ before the main computation. For this article we also ran checks that are not part of the published solutions. A brute-force search over all $1229$ primes below $10^4$ and all $n\le10^6$, with an exact integer cube test for each of the $1.229\cdot10^9$ pairs, finds exactly $28$ solutions: one for each of the $28$ primes $3k^2+3k+1$ below $10^4$, always with $n=k^3$ and $m=k^2(k+1)$, and none for any other prime; the largest of these $n$, $55^3=166375$ for $p=9241$, lies well inside the searched range. A sieve of Eratosthenes up to $10^6$ gives the same count as the three programs, which agree with each other. The C++ program with --limit=1000 and --limit=100000 prints $11$ and $64$, the counts of Section 7, and the first $20$ primes given by the formula agree with OEIS A002407. ^[6]

Footnotes and References

The references below supply the background facts used above. The reduction to two cases, the proof that only consecutive cubes can occur and the operation counts are derived explicitly in this article.

Project Euler 131 — Prime Cube Partnership. The official statement asks for the primes p for which some positive integer n makes n³ + n²p a perfect cube, gives the example 8³ + 8² × 19 = 12³, notes that n is unique for each such prime and that only four such primes lie below one hundred, and asks how many primes below one million have the property.
Wikipedia — Fundamental theorem of arithmetic. Every integer greater than 1 is a product of primes in only one way, apart from the order of the factors, and multiplying two numbers adds the exponents of each prime.
Wikipedia — Euclidean algorithm. The algorithm rests on the principle that the greatest common divisor of two numbers does not change when the larger one is replaced by its difference with the smaller, because both pairs have the same common divisors.
Wikipedia — Sum of two cubes. The difference of two cubes factors as a³ − b³ = (a − b)(a² + ab + b²), which the article checks by multiplying out.
Wikipedia — Cuban prime. The cuban primes of the first kind are the primes that are differences of two successive cubes, of the general form 3y² + 3y + 1: 7, 19, 37, 61, 127, 271, 331, …
OEIS A002407 — Cuban primes. The primes which are the difference of two consecutive cubes, equivalently of the form 1 + 3k(k + 1); a comment describes them as the primes p for which n²(p + n) is a cube for some n > 0.
Wikipedia — Trial division. Trial divisors up to the square root of n suffice, because a composite n has a factor no larger than its square root.
C++ — Euler131.cpp. The linked immutable C++ revision contains parse_arguments, is_prime, solve and run_checkpoints with the checkpoint for the limit 100. It is the implementation record for this article.
Wikipedia — Sieve of Eratosthenes. The sieve marks the multiples of each prime, starting with 2, as composite; its time complexity is O(n log log n), and the basic algorithm needs O(n) memory.

Problem 131 source code

C++

#include <cstdint>
#include <iostream>
#include <string>

namespace {

struct Options {
    int limit = 1000000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.limit >= 2;
}

bool is_prime(int n) {
    if (n < 2) {
        return false;
    }
    if ((n % 2) == 0) {
        return n == 2;
    }
    for (int p = 3; static_cast<std::int64_t>(p) * p <= n; p += 2) {
        if ((n % p) == 0) {
            return false;
        }
    }
    return true;
}

int solve(const int limit) {
    int count = 0;
    for (int a = 1;; ++a) {
        const int p = 3 * a * a + 3 * a + 1;
        if (p >= limit) {
            break;
        }
        if (is_prime(p)) {
            ++count;
        }
    }
    return count;
}

bool run_checkpoints() {
    if (solve(100) != 4) {
        std::cerr << "Checkpoint failed for limit=100" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve(options.limit) << '\n';
    return 0;
}

Python

# Problem 131: Prime cube partnership
# Primes p < 10^6 where there exists n such that n^3 + n^2*p is a perfect cube.
# This simplifies to p = 3a^2 + 3a + 1 (difference of consecutive cubes).

def solve():
    limit = 1000000
    count = 0
    a = 1
    while True:
        p = 3*a*a + 3*a + 1
        if p >= limit: break
        # check primality
        if p > 1:
            is_p = True
            d = 2
            while d*d <= p:
                if p % d == 0: is_p = False; break
                d += 1
            if is_p: count += 1
        a += 1
    print(count)

solve()

Java

public class Euler131 {
    static boolean isPrime(int n) {
        if (n < 2)
            return false;
        if (n % 2 == 0)
            return n == 2;
        for (int p = 3; (long) p * p <= n; p += 2)
            if (n % p == 0)
                return false;
        return true;
    }

    public static void main(String[] args) {
        int limit = 1000000, count = 0;
        for (int a = 1;; a++) {
            int p = 3 * a * a + 3 * a + 1;
            if (p >= limit)
                break;
            if (isPrime(p))
                count++;
        }
        System.out.println(count);
    }
}