Problem 136: Singleton Difference

Project Euler Problem 136 Solution

EulerSolve provides an optimized solution for Project Euler Problem 136, Singleton Difference, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We seek the number of integers $n \lt 5 \times 10^7$ for which the equation $x^2 - y^2 - z^2 = n$ has exactly one solution in positive integers, under the extra condition that $x$, $y$, and $z$ are consecutive terms of an arithmetic progression. Because the progression condition ties the three variables together, the real task is not an unrestricted three-variable search but a structured counting problem over admissible factorizations of $n$. The C++, Python, and Java implementations do not test each $n$ separately. Instead they enumerate every valid progression once, convert it into the corresponding value of $n$, and keep only the information needed to distinguish “zero solutions”, “one solution”, and “more than one solution”. Mathematical Approach The key step is to rewrite the arithmetic progression in terms of its middle term and common difference, then convert the resulting quadratic expression into a product with simple congruence and inequality constraints....

Detailed mathematical approach

Problem Summary

Three positive integers $x$, $y$, $z$ are consecutive terms of an arithmetic progression, so that $y-x=z-y$. For a positive integer $n$ we count the triples of this kind that solve the equation

$$x^2-y^2-z^2=n.$$

For $n=20$ there is exactly one solution, $13^2-10^2-7^2=20$, and there are $25$ values of $n$ below $100$ with exactly one solution. We must find how many values of $n$ below fifty million, that is $n\lt5\cdot10^7$, have exactly one solution. ^[1]

The whole solution in one sentence: writing the progression as $y+d$, $y$, $y-d$ turns the equation into $n=y(4d-y)$, which makes every solution a factorization of $n$ with two simple side conditions, so one pass over all such factor pairs below $5\cdot10^7$, with a counter for each $n$ that stops at $2$, finds the values that are reached exactly once.

Mathematical Approach

We list small cases, turn the equation into a factorization of $n$ with two side conditions and prove exactly which $n$ have one solution. Then we count all $n$ at once by running through the factor pairs, as the programs do.

1. Small progressions

Three consecutive terms of an arithmetic progression have a middle term and a common difference, the constant step from one term to the next. ^[2] We write them as $x=y+d$, $y$ and $z=y-d$, so that $d=x-y=y-z$ is the step of the progression read backwards, from $z$ to $x$; at first $d$ could be positive, zero or negative. When the squares are expanded, the terms $d^2$ cancel:

$$\begin{aligned}n&=(y+d)^2-y^2-(y-d)^2\\&=4yd-y^2=y\,(4d-y).\end{aligned}$$

Section 2 shows that only $d\ge1$ and $d\lt y\lt4d$ can occur. For the three smallest differences this gives:

$d=1$: $y=2,3$ give $n=4,3$;
$d=2$: $y=3,\dots,7$ give $n=15,16,15,12,7$;
$d=3$: $y=4,\dots,11$ give $n=32,35,36,35,32,27,20,11$.

The pair $y=10$, $d=3$ is the progression $13,10,7$ of the statement. Some values repeat: $15$ comes from $5,3,1$ and from $7,5,3$, and once more from $19,15,11$ with $d=4$.

2. The progression decreases

In $n=y(4d-y)$ both $n$ and $y$ are positive, so $4d-y\gt0$, that is $y\lt4d$. In particular $d\gt0$, and the progression decreases: $x\gt y\gt z$. The third term $z=y-d$ is positive exactly when $y\gt d$. Hence the solutions are the pairs of integers $(y,d)$ with

$$d\ge1,\qquad d\lt y\lt 4d,$$

and conversely every such pair gives positive $x$, $y$, $z$ and a positive $n$.

3. Solutions are factor pairs

Put $u=y$ and $v=4d-y$. Then $n=uv$ and $u+v=4d$. Two integers $a$ and $b$ are congruent modulo $m$, written $a\equiv b\pmod m$, when $m$ divides $a-b$, ^[3] so $u+v\equiv0\pmod4$. The conditions of Section 2 become $v\ge1$ and, since $y\gt d$ means $4u\gt u+v$, also $v\lt3u$.

Claim. The number of solutions for $n$ equals the number of ordered pairs $(u,v)$ of positive integers with

$$\begin{aligned}uv&=n,\\u+v&\equiv0\pmod4,\\v&\lt3u.\end{aligned}$$

Proof. A solution is determined by its middle term $y$ and its difference $d=x-y$, and Section 2 turns it into such a pair. Conversely, for such a pair $d=(u+v)/4$ is a positive integer, and $y=u$ satisfies $d\lt y\lt4d$, because $v\lt3u$ gives $u+v\lt4u$ and $v\ge1$ gives $u\lt u+v$. The resulting progression is

$$\begin{aligned}x&=\frac{5u+v}{4},\\y&=u,\\z&=\frac{3u-v}{4}.\end{aligned}$$

The two maps $(y,d)\mapsto(y,4d-y)$ and $(u,v)\mapsto(u,(u+v)/4)$ undo each other, so they match solutions and pairs one to one.

The pair has a classical meaning: $n=(2d)^2-(y-2d)^2$, and $u$ and $v$ are the sum and the difference of $2d$ and $y-2d$, as in the identity $a^2-b^2=(a+b)(a-b)$. ^[4]

4. Example: the divisors of 20 and of 15

The ordered pairs with $uv=20$ are $(1,20)$, $(2,10)$, $(4,5)$, $(5,4)$, $(10,2)$ and $(20,1)$, with the sums $21,12,9,9,12,21$. Only $(2,10)$ and $(10,2)$ have a sum divisible by $4$, and $(2,10)$ fails $v\lt3u$, because $10\ge6$. The single remaining pair $(10,2)$ gives $d=3$ and the progression $13,10,7$, so $20$ has exactly one solution.

For $n=15$ the pairs $(1,15)$, $(3,5)$, $(5,3)$ and $(15,1)$ all have a sum divisible by $4$, and only $(1,15)$ fails $v\lt3u$. The other three give $d=2$, $2$ and $4$, and the progressions $5,3,1$, $7,5,3$ and $19,15,11$: three solutions. The same test for every $n\lt100$ finds $56$ values without a solution, $25$ with one, $8$ with two, $8$ with three and $2$ with four, namely $63$ and $99$.

5. Odd values, and even values not divisible by 4

If $n$ is odd, then $u$ and $v$ are odd, so each leaves the remainder $1$ or $3$ after division by $4$, and $u+v$ is divisible by $4$ exactly when the two remainders differ, since $1+3=4$ while $1+1$ and $3+3$ are not multiples of $4$. Congruences may be multiplied, ^[3] and $1\cdot3=3$ while $1\cdot1=1$ and $3\cdot3=9\equiv1\pmod4$, so the remainders differ exactly when $n\equiv3\pmod4$. Hence $n\equiv1\pmod4$ has no solution, and for $n\equiv3\pmod4$ every factorization satisfies the congruence, so only $v\lt3u$ has to be checked. If $n$ is even but not divisible by $4$, exactly one of $u$ and $v$ is even, $u+v$ is odd, and there is no solution either.

Now let $n\equiv3\pmod4$. The pair $(n,1)$ always works. If $n=p$ is prime, the only other pair is $(1,p)$, which fails because $p\ge3$. If $n$ is composite, let $p$ be its smallest prime factor and $k=n/p$. Every prime factor of $k$ is at least $p$, so $k\ge p$, and the pair $(k,p)$ satisfies $p\le k\lt3k$: a second solution. So an odd $n$ has exactly one solution if and only if it is a prime $p\equiv3\pmod4$, and that solution has $y=p$ and $d=(p+1)/4$; for $p=7$ it is $9,7,5$.

6. Multiples of 4

If $4$ divides $n$, then $u$ and $v$ are both even, since otherwise one factor would be odd, the other a multiple of $4$, and the sum odd. With $u=2a$ and $v=2b$ we get $n=4ab$, the congruence says that $a+b$ is even, and $v\lt3u$ becomes $b\lt3a$. Let $m=n/4=ab$.

If $m$ is odd, then $a$ and $b$ are odd, and the pairs with $ab=m$ and $b\lt3a$ behave as in Section 5: $(m,1)$ always works, it is the only pair for $m=1$ and for an odd prime $m=p$, and a composite $m$ has the second pair $(k,p)$. So $n=4$ and $n=4p$ have exactly one solution, for $n=4p$ the one with $y=2p$ and $d=(p+1)/2$, as in the statement's $20=4\cdot5$.
If $m$ is even but not divisible by $4$, the factors $a$ and $b$ cannot be both odd or both even, so there is no solution.
If $4$ divides $m$, then $a=2a'$ and $b=2b'$ are even, $n=16a'b'$, and only $b'\lt3a'$ remains, with no condition on parity. With $m'=n/16=a'b'$, again $(m',1)$ always works and is the only pair for $m'=1$, and a composite $m'$ has a second pair. For a prime $m'=p$, the pair $(1,p)$ also works exactly when $p\lt3$, that is $p=2$: so $n=32$ has the two solutions $7,4,1$ and $11,8,5$, while $n=16p$ for an odd prime $p$ has exactly one, with $y=4p$ and $d=p+1$.

7. Which values have exactly one solution

Sections 5 and 6 together prove the following.

Theorem. A positive integer $n$ has exactly one solution if and only if $n$ is a prime $p\equiv3\pmod4$, or $n=4p$ or $n=16p$ for an odd prime $p$, or $n=4$ or $n=16$.

Below $100$ these are the $13$ primes $3,7,11,19,23,31,43$, $47,59,67,71,79,83$, the $8$ numbers $4p$ from $12$ to $92$, the $2$ numbers $48$ and $80$, and the values $4$ and $16$: $13+8+2+2=25$, the number given in the statement. In general, for $N\gt16$ the count below $N$ is the number of primes $p\equiv3\pmod4$ below $N$, plus the number of odd primes below $N/4$, plus the number of odd primes below $N/16$, plus $2$. For $N=1000$ this is $87+52+17+2=158$, and for $N=10^6$ it is $39322+22043+6274+2=67641$. The theorem reduces the problem to counting primes, which a sieve of Eratosthenes does for all numbers below $N$ in $O(N\log\log N)$ operations. ^[5] The published programs count the solutions directly instead; the theorem serves as an independent check.

8. Counting all values at once

Instead of factoring every $n$ below the limit $N=5\cdot10^7$, we run through all pairs $(u,v)$ of Section 3 with $uv\lt N$ and add $1$ to a counter of $n=uv$ at each pair. By the claim of Section 3, the counter of $n$ then holds its number of solutions; the programs stop every counter at $2$, as Section 10 explains. For a fixed $u\ge1$, three conditions describe the admissible $v$:

$v\equiv-u\pmod4$: the smallest positive such $v$ is $4$, $3$, $2$ or $1$ when $u$ leaves the remainder $0$, $1$, $2$ or $3$ after division by $4$, and the others follow in steps of $4$;
$v\lt3u$, that is $v\le3u-1$;
$uv\le N-1$, that is $v\le\lfloor(N-1)/u\rfloor$, where $\lfloor t\rfloor$ is $t$ rounded down.

So $v$ runs in steps of $4$ from its smallest value up to the smaller of the two bounds. Since $(u,1)$ is a pair for every $u\equiv3\pmod4$, the values of $u$ run from $1$ to $N-1$.

9. Worked example: $N=25$

For $N=25$ the procedure visits these $12$ pairs, written as $u\cdot v=n$:

$2\cdot2=4$, $3\cdot1=3$, $3\cdot5=15$, $4\cdot4=16$, $5\cdot3=15$, $6\cdot2=12$,
$7\cdot1=7$, $10\cdot2=20$, $11\cdot1=11$, $15\cdot1=15$, $19\cdot1=19$, $23\cdot1=23$.

Every other $u$ gives no pair: $u=1$ would need $v=3$, but $3u-1=2$, and $u=8$ would need $v=4$, but $\lfloor24/8\rfloor=3$. The value $15$ is reached three times, and the nine values $3,4,7,11,12,16,19,20,23$ exactly once. This matches Section 7: five primes $p\equiv3\pmod4$, $12=4\cdot3$, $20=4\cdot5$, $4$ and $16$. For $N=100$ the procedure visits $73$ pairs, one for each of the $25\cdot1+8\cdot2+8\cdot3+2\cdot4=73$ solutions counted in Section 4, and leaves exactly $25$ counters at $1$.

10. Counters that stop at 2

The question only distinguishes no solution, exactly one, and more than one. So the counter of $n$ is increased only while it is below $2$; it ends at $\min(2,r)$, where $r$ is the number of solutions, and a final $1$ still means exactly one solution. One byte per value is therefore enough, whatever the size of $r$.

How the Code Works

The C++ program ^[6] keeps its settings in the structure Options: limit, by default $50\,000\,000$, and run_checkpoints, by default true. parse_arguments accepts --skip-checkpoints and --limit= followed by decimal digits, which parse_int_after_prefix accumulates in an int without a range check; any other argument, or a limit below $2$, ends the program with exit code $1$.

solve(limit) carries out Sections 8 and 10. It creates the table std::vector<std::uint8_t> count with one zero byte for each $n$ from $0$ to $N-1$, where $N$ is the limit. For each u from $1$ to $N-1$ it computes vmax, the bound $\lfloor(N-1)/u\rfloor$, and upper, the smaller of vmax and $3u-1$; for the default limit the test upper < 1 never succeeds, since both bounds are at least $1$. The first v comes from (-u) % 4. In C++ the quotient of two integers is rounded towards zero, so the remainder of a negative number is negative or zero; ^[7] for $u=5$ the expression gives $-1$, and the program adds $4$ to obtain $3$. A remainder of $0$ is replaced by $4$, because $v$ must be positive. Then v runs in steps of $4$ up to upper, and the counter of $n=uv$ is increased if it is below $2$. A last loop counts the $n$ from $1$ to $N-1$ whose counter equals $1$. run_checkpoints requires solve(100) to return $25$, the value from the statement, and otherwise ends the program with exit code $2$; finally main prints solve(options.limit).

For the default limit every intermediate value fits easily into a 32-bit int: $uv\le N-1$ and $3u-1\le149\,999\,996$. Larger values after --limit= lead to signed overflow, which is undefined behaviour in C++: ^[7] a value above $2^{31}-1$ overflows already while its digits are parsed, and from --limit=715827884 on, the product $3u$ overflows at $u=715\,827\,883$. This $u$ is a prime $\equiv3\pmod4$, and in a test with that limit the program missed its single pair $(u,1)$ and printed a count one too small.

The Python version has the limit $50\,000\,000$ built into solve(), without options or a checkpoint, and uses a bytearray as the table. Python's % gives a result with the sign of the divisor, ^[7] so (-u) % 4 already lies between $0$ and $3$ and only the replacement of $0$ by $4$ is needed; Python integers cannot overflow. solve returns the count as a string, which the main block prints. The Java version repeats the loops in main with a byte[] table and the same fixed limit; its remainder has the sign of the dividend, as in C++, hence ((-u) % 4 + 4) % 4. A Java byte is signed, which is harmless, since no counter exceeds $2$.

Complexity and Verification

The outer loop makes $N-1$ steps. The inner loop makes one step for each pair $(u,v)$, that is, for each solution of some $n\lt N$. For $u$ up to about $\sqrt{N/3}$, where $3u-1$ is the smaller bound, there are about $3u/4$ values of $v$, about $N/8$ in all. For larger $u$ there are about $N/(4u)$, which add up to $\frac N4(H_{N-1}-H_s)$ with $s\approx\sqrt{N/3}$, where $H_k=1+\frac12+\dots+\frac1k$ is the $k$-th harmonic number. Since $H_k=\ln k+\gamma+O(1/k)$, ^[8] this is about

$$\frac N4\Bigl(\ln N-\ln\sqrt{N/3}\Bigr)=\frac N8\ln(3N).$$

The rounding in each count of $v$ changes the total by $O(N)$, so the inner loop makes $\frac N8\ln N+O(N)$ steps, and the running time is $\Theta(N\log N)$. For $N=5\cdot10^7$ the inner loop makes exactly $118\,626\,557$ steps, about $2.4$ per value of $n$, besides the $49\,999\,999$ steps of the outer loop and the final scan of the table. The table takes one byte per value, $5\cdot10^7$ bytes or about $48$ MiB, in all three languages. In tests on a desktop computer the C++ program, compiled with -O2, took about a quarter of a second, the Java program under half a second on OpenJDK 17 and the Python program about $13$ seconds on CPython 3.11.

The C++ program checks the statement's value $25$ for $N=100$ before the main run, and with --limit= it reproduces the counts $9$, $158$ and $67641$ of Sections 7 and 9 for $N=25$, $1000$ and $10^6$. Separate computations for this article, not part of the published solutions, confirmed the mathematics: a direct enumeration of all $y\lt3000$ and all $d$ with $|d|\lt y$, including zero and negative differences, agrees with the pairs of Section 3 for every $n\lt3000$ and misses no solution, because $n=y(4d-y)\ge y$ when $n\gt0$; the counters of the published loop equal $\min(2,r)$, with $r$ found from the divisors of $n$, for every $n\lt20\,000$; and the theorem of Section 7 picks out exactly the values with one solution for every $n\lt10^6$. For $N=5\cdot10^7$, a separate program that enumerates the progressions directly with 64-bit arithmetic and a prime sieve that applies the theorem both give the same count as the three published programs; the direct enumeration finds $118\,626\,557$ progressions, the number of inner steps found above.

Footnotes and References

The references below supply the background facts used above. The correspondence between solutions and factor pairs, the theorem on the values with exactly one solution and the operation count are derived explicitly in this article.

Project Euler 136 — Singleton Difference. The official statement: the positive integers x, y and z are consecutive terms of an arithmetic progression, and the equation x² − y² − z² = n has exactly one solution for n = 20, namely 13² − 10² − 7² = 20; twenty-five values of n below one hundred have a unique solution, and the task asks how many values of n less than fifty million have exactly one solution.
Wikipedia — Arithmetic progression. An arithmetic progression is a sequence in which the difference between each term and the one before it stays the same throughout; this constant is the common difference.
Wikipedia — Modular arithmetic. Two integers are congruent modulo m when their difference is a multiple of m; congruence is compatible with addition, subtraction and multiplication.
Wikipedia — Difference of two squares. Every difference of two squares factors as the sum of the two numbers times their difference, a² − b² = (a + b)(a − b).
Wikipedia — Sieve of Eratosthenes. The sieve finds all primes up to a bound by marking the multiples of each prime as composite; it needs O(n log log n) operations and O(n) memory.
C++ — Euler136.cpp. The linked immutable C++ revision contains parse_arguments, solve and run_checkpoints with the checkpoint for the values below 100. It is the implementation record for this article.
Remainders and overflow: cppreference — Arithmetic operators (integer division truncates towards zero, and signed overflow is undefined behaviour); Java Language Specification, §15.17.3 (a remainder can be negative only if the dividend is negative); Python Language Reference, §6.7 (the result of the modulo operator has the sign of its second operand).
Wikipedia — Harmonic number. The harmonic numbers Hₙ = 1 + 1/2 + … + 1/n satisfy Hₙ = ln n + γ + 1/(2n) − …, with the Euler–Mascheroni constant γ ≈ 0.5772.

Problem 136 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

struct Options {
    int limit = 50000000;
    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;
}

int solve(const int limit) {
    std::vector<std::uint8_t> count(static_cast<std::size_t>(limit), 0);

    for (int u = 1; u < limit; ++u) {
        const int vmax = (limit - 1) / u;
        const int upper = std::min(vmax, 3 * u - 1);
        if (upper < 1) {
            continue;
        }

        int v = (-u) % 4;
        if (v < 0) {
            v += 4;
        }
        if (v == 0) {
            v = 4;
        }

        for (; v <= upper; v += 4) {
            const int n = u * v;
            if (count[static_cast<std::size_t>(n)] < 2) {
                ++count[static_cast<std::size_t>(n)];
            }
        }
    }

    int answer = 0;
    for (int n = 1; n < limit; ++n) {
        if (count[static_cast<std::size_t>(n)] == 1) {
            ++answer;
        }
    }
    return answer;
}

bool run_checkpoints() {
    if (solve(100) != 25) {
        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

def solve():
    limit = 50000000
    count = bytearray(limit)
    for u in range(1, limit):
        vmax = (limit - 1) // u
        upper = min(vmax, 3 * u - 1)
        if upper < 1:
            continue
        v = (-u) % 4
        if v == 0:
            v = 4
        while v <= upper:
            n = u * v
            if count[n] < 2:
                count[n] += 1
            v += 4
    answer = 0
    for n in range(1, limit):
        if count[n] == 1:
            answer += 1
    return str(answer)

if __name__ == '__main__':
    print(solve())

Java

public class Euler136 {
    public static void main(String[] args) {
        int limit = 50000000;
        byte[] count = new byte[limit];
        for (int u = 1; u < limit; u++) {
            int vm = Math.min((limit - 1) / u, 3 * u - 1);
            if (vm < 1)
                continue;
            int v = ((-u) % 4 + 4) % 4;
            if (v == 0)
                v = 4;
            for (; v <= vm; v += 4) {
                int n = u * v;
                if (count[n] < 2)
                    count[n]++;
            }
        }
        int ans = 0;
        for (int n = 1; n < limit; n++)
            if (count[n] == 1)
                ans++;
        System.out.println(ans);
    }
}