Problem 283: Integer Sided Triangles with Integral Area/perimeter Ratio

Project Euler Problem 283 Solution

EulerSolve provides an optimized solution for Project Euler Problem 283, Integer Sided Triangles with Integral Area/perimeter Ratio, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let a triangle with integer sides $a,b,c$ have area $A$ and perimeter $P=a+b+c$. We want all triangles for which $$\frac{A}{P}=k$$ is an integer. For each $k$, let $p(k)$ be the sum of the corresponding perimeters. The program computes $$\sum_{k=1}^{1000} p(k),$$ but the final numeric answer is intentionally omitted here. Mathematical Approach 1) Replace side lengths by semiperimeter gaps. Let $$s=\frac{a+b+c}{2},\qquad x=s-a,\qquad y=s-b,\qquad z=s-c.$$ Then $x,y,z$ are positive integers, and conversely $$a=y+z,\qquad b=z+x,\qquad c=x+y,\qquad s=x+y+z.$$ So every integer triangle corresponds to one positive integer triple $(x,y,z)$. If we sort the sides as $a \ge b \ge c$, then $x \le y \le z$, which is exactly the ordering used in the code to avoid duplicates. 2) Heron's formula gives a Diophantine equation. Heron's formula is $$A^2=s(s-a)(s-b)(s-c)=sxyz.$$ The condition $A/P=k$ means $$A=k(a+b+c)=2ks.$$ Squaring and cancelling one factor of $s$ gives $$xyz=4k^2s=4k^2(x+y+z).$$ If we define $$n=(2k)^2=4k^2,$$ the triangle problem becomes the pure integer equation $$xyz=n(x+y+z),\qquad x \le y \le z.$$ The perimeter of the recovered triangle is $$P=2s=2(x+y+z),$$ which is why the program adds $2(x+y+z)$ for each valid triple. 3) Why this is a bijection. Starting from a triangle, we get one positive triple $(x,y,z)$....

Detailed mathematical approach

Problem Summary

A triangle with integer side lengths $a$, $b$, $c$ has the perimeter $P=a+b+c$ and an area $A$, which Heron's formula computes from the sides. ^[2] We must find the sum of the perimeters of all such triangles for which $A/P$ is a positive integer $k\le1000$; a triangle is determined by its three side lengths, so each one is counted once. The statement gives two examples: the triangle with sides $6$, $8$, $10$ has perimeter and area both equal to $24$, so its ratio is $1$, and the triangle with sides $13$, $14$, $15$ has perimeter $42$ and area $84$, so its ratio is $2$. ^[1] If $p(k)$ denotes the sum of the perimeters of the triangles with ratio exactly $k$, the requested number is $p(1)+p(2)+\dots+p(1000)$.

The whole solution in one sentence: the ratio $A/P$ is half the radius of the inscribed circle, so we describe each triangle by the three whole-number distances from its corners to the points where that circle touches the sides, which must satisfy one simple equation, and once the smallest of these distances is chosen, the other two are read off from the factor pairs of a single known number.

Mathematical Approach

We turn the ratio into a statement about the inscribed circle, then into an equation in whole numbers, and finally into a short list of divisors, trying every step on the smallest case $k=1$.

1. The ratio is half the inradius

Every triangle has an inscribed circle, the incircle, which touches all three sides; its radius $r$ is the inradius and its centre $I$ the incentre. The segments from $I$ to the three vertices cut the triangle into three smaller triangles. Each of them has one side of the big triangle as its base, and its height is $r$, because the radius to the touching point is perpendicular to that side. Adding the three areas gives, with the semiperimeter $s=P/2$, ^[3]

$$A=\frac{ra}{2}+\frac{rb}{2}+\frac{rc}{2}=rs,$$

and therefore $A/P=rs/(2s)=r/2$. So $A/P=k$ means exactly $r=2k$: we need all integer triangles whose inradius is one of the even numbers $2,4,\dots,2000$. For the triangle $6,8,10$ we have $s=12$ and $A=24$, so $r=2$; for $13,14,15$ we have $s=21$ and $A=84$, so $r=4$. Triangles whose area equals their perimeter, the case $k=1$, are called equable, and every equable triangle has inradius $2$. ^[4]

2. The perimeter is even

Heron's formula can be written without fractions as ^[2] ^[5]

$$16A^2=P\,(P-2a)(P-2b)(P-2c).$$

If $A/P=k$ is a whole number, then $A=kP$ is a whole number, so our triangles are Heronian triangles, with integer sides and integer area, and the left side is even. If $P$ were odd, the three numbers $P-2a$, $P-2b$, $P-2c$ would be odd too, and the right side would be a product of four odd numbers, which is odd. Hence $P$ is even, and the semiperimeter $s=P/2$ is a whole number. For example, the triangle $2,3,4$ has the odd perimeter $9$, and $16A^2=9\cdot5\cdot3\cdot1=135$ is odd, so its area is not a whole number and its ratio is not one either.

3. Tangent lengths

The two tangent segments from a point to a circle have the same length. So every vertex has one tangent length, its distance to the two touching points next to it, and every side is the sum of the tangent lengths at its two ends. Call them $x$, $y$, $z$ at the vertices opposite the sides $a$, $b$, $c$. Then

$$a=y+z,\quad b=z+x,\quad c=x+y,$$

so $s=x+y+z$, and $x=s-a$, $y=s-b$, $z=s-c$. ^[3] For our triangles these are integers by Section 2, and they are positive because each side is shorter than the sum of the other two. Conversely, any positive integers $x,y,z$ give the sides of a triangle, since for instance $(z+x)+(x+y)$ exceeds $y+z$ by $2x$. Sorting the sides as $a\ge b\ge c$ is the same as asking for $x\le y\le z$, because $a-b=y-x$ and $b-c=z-y$. So the integer triangles with even perimeter correspond one to one to the triples $x\le y\le z$ of positive integers, and all the triangles we want are among them. The triangle $10,8,6$ has $s=12$ and $(x,y,z)=(2,4,6)$; the triangle $15,14,13$ has $s=21$ and $(6,7,8)$.

4. The equation

In tangent lengths Heron's formula reads $A^2=s\,xyz$, and Section 1 gives $A=rs$. Hence $r^2s^2=s\,xyz$, and dividing by $s$ gives, with $n=r^2=4k^2$,

$$xyz=n\,(x+y+z).$$

This is the formula $r^2=(s-a)(s-b)(s-c)/s$ for the inradius. ^[3] Conversely, if positive integers satisfy the equation, then $A^2=s\,xyz=n\,s^2$, so $A=2ks$ and $A/P=2ks/(2s)=k$. The task is therefore: for each $k\le1000$, find all integer triples $x\le y\le z$ with $xyz=4k^2(x+y+z)$, and add the perimeter $2(x+y+z)$ of each. Both examples pass: $2\cdot4\cdot6=48=4\cdot12$ for $k=1$, and $6\cdot7\cdot8=336=16\cdot21$ for $k=2$.

5. The smallest tangent length is small

Solved for $z$, the equation reads $z\,(xy-n)=n\,(x+y)$. The right side is positive, so $xy\gt n$, and once $x$ and $y$ are chosen, $z=n(x+y)/(xy-n)$ is fixed. The smallest number is bounded too: because $x\le y\le z$, each of $x^2$, $xy$ and $xz$ is at most $yz$, so $x(x+y+z)\le3yz$. Multiplying by $x$ and using the equation gives

$$x^2\,(x+y+z)\le3xyz=3n\,(x+y+z),$$

so $x^2\le3n=12k^2$. Equality would need $x=y=z$ and $x^2=12k^2$, which is impossible because $12k^2=3(2k)^2$ is not a square. Hence $x\le\lfloor2\sqrt3\,k\rfloor$: there are at most $3$ values of $x$ for $k=1$, $6$ for $k=2$ and $3464$ for $k=1000$. A bound for $y$ follows in the same way: $z\ge y$ means $n(x+y)\ge y(xy-n)$, that is $xy^2-2ny-nx\le0$, and solving this quadratic inequality gives

$$y\le\frac{n+\sqrt{n^2+nx^2}}{x}.$$

These two bounds already give a correct but slow method: try every $x$ and every $y$ in range and compute $z$. The C++ program keeps it as a check.

6. One divisor instead of two unknowns

The slow method spends almost all of its time on values of $y$ that fail. For a fixed $x$ the equation can be rewritten so that $y$ and $z$ appear in a product. Expanding the left side and cancelling shows

$$\begin{aligned}(xy-n)(xz-n)&-n\,(x^2+n)\\&=x\bigl(xyz-n(x+y+z)\bigr).\end{aligned}$$

Since $x\gt0$, the equation of Section 4 is therefore equivalent to

$$(xy-n)(xz-n)=n\,(x^2+n).$$

Put $m=x^2+n$, $M=nm$, $u=xy-n$ and $v=xz-n$. Claim: for fixed $k$ and $x$, the solutions $x\le y\le z$ correspond one to one to the positive divisors $u$ of $M$ with $u^2\le M$ for which $x$ divides both $u+n$ and $v+n$, where $v=M/u$, and $(u+n)/x\ge x$. Indeed, for a solution we have $u\gt0$ by Section 5 and $v\ge u$ because $z\ge y$, so $u$ is a divisor of $M$ with $u^2\le uv=M$, and $u+n=xy$, $v+n=xz$ are multiples of $x$. Conversely, for such a divisor the numbers $y=(u+n)/x$ and $z=(v+n)/x$ are positive integers that satisfy the product equation, hence the equation of Section 4; $u\le v$ gives $y\le z$, and the last condition is $y\ge x$. Different divisors give different $y$. Every triangle is therefore found exactly once, at its own smallest tangent length $x$, and nothing else is found.

7. Worked example: the five equable triangles

For $k=1$ we have $n=4$, and Section 5 allows $x=1,2,3$.

$x=1$: $m=5$ and $M=20$, with the divisors $1,2,4,5,10,20$. Those with $u^2\le20$ are $u=1,2,4$, with $v=20,10,5$. Every number is a multiple of $1$, so all three give solutions: $(y,z)=(5,24)$, $(6,14)$ and $(8,9)$.
$x=2$: $m=8$ and $M=32$, with the divisors $1,2,4,8,16,32$. $u=1$ fails, because $u+n=5$ is odd; $u=2$, $v=16$ gives $(y,z)=(3,10)$, and $u=4$, $v=8$ gives $(4,6)$.
$x=3$: $m=13$ and $M=52$, with the divisors $1,2,4,13,26,52$. $u=1$ and $u=4$ fail, because $5$ and $8$ are not multiples of $3$; $u=2$ gives $y=2\lt x$, the triple $(2,3,10)$ again, which belongs to $x=2$.

The five triples $(1,5,24)$, $(1,6,14)$, $(1,8,9)$, $(2,3,10)$ and $(2,4,6)$ give the triangles $29,25,6$, $20,15,7$, $17,10,9$, $13,12,5$ and $10,8,6$, with the perimeters $60$, $42$, $36$, $30$ and $24$. Heron's formula confirms that each of them has an area equal to its perimeter. These are exactly the five equable triangles with integer sides, a list proved complete by Whitworth and Biddle in 1904, ^[4] and $p(1)=192$.

8. Why both divisibility tests are needed

One might expect $x\mid u+n$ to imply $x\mid v+n$. Writing $u+n=qx$ and substituting $u=qx-n$ into $uv=nx^2+n^2$ gives $qxv-nv=nx^2+n^2$, that is

$$n\,(v+n)=x\,(qv-nx).$$

So $x$ always divides $n(v+n)$, and if $x$ and $n$ have no common prime factor, then $x\mid v+n$ by the general form of Euclid's lemma. ^[6] But whenever $x$ shares a prime factor with $n=4k^2$, for instance whenever $x$ is even, the second test is not automatic. For $k=3$, $n=36$ and $x=2$ we get $M=36\cdot40=1440$; the divisor $u=32$ has $u+n=68=2\cdot34$, but $v=45$ has $v+n=81$, which is odd, so $z$ would be $81/2$. Both conditions must therefore be tested.

9. A second example and the first ten values

For $k=2$ we have $n=16$ and $x\le6$, and the same procedure gives $18$ triangles: $5$ with $x=1$, $5$ with $x=2$, $4$ with $x=3$, $3$ with $x=4$, none with $x=5$ and one with $x=6$, the triple $(6,7,8)$ of the triangle $15,14,13$ from the statement. For it $u=6\cdot7-16=26$ and $v=6\cdot8-16=32$, and indeed $26\cdot32=832$, which is $16\cdot52=n\,(x^2+n)$. The eighteen perimeters add up to $p(2)=2604$. For $k=1,2,\dots,10$:

the numbers of triangles are $5$, $18$, $45$, $45$, $52$, $139$, $80$, $89$, $184$, $145$;
the sums of their perimeters are $192$, $2604$, $14848$, $40212$, $82242$, $262840$, $309306$, $494586$, $1234250$, $1340706$.

They add up to $140098$ for $k\le5$ and to $3781786$ for $k\le10$, the two checkpoints stored in the C++ program. The counts jump around, $139$ for $k=6$ but $80$ for $k=7$: they depend on the divisors of the numbers $M$, and $n=4k^2$ has $15$ divisors for $k=6$ but only $9$ for $k=7$.

10. How large the triangles get

For a fixed $x$, $x+y+z=x+(u+v+2n)/x$. With the product $uv=M$ fixed, $(u-1)(v-1)\ge0$ gives $u+v\le1+M$, and with $M=nx^2+n^2$ this gives

$$x+y+z\le(n+1)\,x+\frac{(n+1)^2}{x}.$$

The right side is at most $(n+1)(n+2)$, its value at $x=1$: the difference is $(n+1)(x-1)(n+1-x)/x$, which is not negative because $1\le x\le\sqrt{3n}\le n+1$. So every triangle with ratio $k$ has a perimeter of at most $2(n+1)(n+2)$, and the triple $(1,n+1,n^2+2n)$, from $x=1$ and $u=1$, attains it. For $k=1$ this is the triangle $29,25,6$ with perimeter $60$; for $k=1000$ the bound is about $3.2\cdot10^{13}$. This bound fixes the size of the numbers in the programs and limits the brute-force check below.

11. Divisors from a table of smallest prime factors

To list the divisors of $M=nm$ we need its prime factorization, and both factors are small: $n=4k^2\le4\cdot10^6$ and, by Section 5, $m=x^2+n\lt4n\le16\cdot10^6$. So a table of the smallest prime factor of every number up to $16\,000\,000$ is built once, with a variant of the sieve of Eratosthenes: the numbers are visited in increasing order, an entry that is still empty belongs to a prime $p$, and $p$ is written into its own entry and into every still-empty entry among its multiples $p^2,p^2+p,\dots$. ^[7] A composite number is a multiple of its smallest prime factor $p$ and at least $p^2$, so it receives $p$, and no smaller prime divides it. A number is factorized by dividing it repeatedly by its table entry. The factorizations of $n$ and $m$ are merged, adding the exponents of common primes, and the divisors are generated prime by prime: starting from the list $1$, each prime power $p^e$ appends the list as it was before this prime, multiplied by $p$, by $p^2$, and so on up to $p^e$. Every divisor has one exponent of $p$ and one remaining factor in that old list, so a number $p_1^{e_1}\cdots p_t^{e_t}$ receives all its $(e_1+1)\cdots(e_t+1)$ divisors, ^[8] each exactly once. For $k=1$ and $x=1$, the factors $n=2^2$ and $m=5$ give first $1,2,4$ and then $5,10,20$, the divisors used in Section 7.

How the Code Works

The C++ program ^[9] starts in main, which reads the optional arguments --k-limit=, the largest ratio $K$ (default $1000$), and --skip-checkpoints. It sets n_max to $4(2K)^2$, which is $16\,000\,000$ for $K=1000$, builds the table of Section 11 with build_spf, runs run_checkpoints unless told to skip it, calls solve_fast and prints the total with to_string_u128. The checkpoints always run up to $k=10$, so for a limit $K\lt10$ the table is too small for them, and --skip-checkpoints must be given as well. The helpers factorize_int, merge_factorizations and generate_divisors are the three steps of Section 11; a factorization is a list of pairs (prime, exponent) in increasing order of the primes, and generate_divisors keeps the length of the old list in base.

The main work happens in solve_range, which adds the contributions of a range of ratios. For each $k$ it sets r = 2k, the inradius of Section 1, and n = r*r, and factorizes $n$ once. The bound x_max of Section 5 is the floor of a long double square root of $3n$; a check over all $k\le1000$ confirms that rounding never changes it. For each $x$ it forms $m$, factorizes it, merges the result with the factorization of $n$, generates the divisors of $M$ and runs through them. A divisor is skipped if u * u > M, which would mean $y\gt z$; if $x$ does not divide both $u+n$ and $v+n$ (Section 8); or if $y\lt x$. Every remaining divisor adds the perimeter $2(x+y+z)$ to a total of type unsigned __int128, a 128-bit integer type of GCC and Clang. ^[10] solve_fast splits $1\le k\le K$ into consecutive blocks, evaluates each with solve_range and adds the block totals.

solve_bruteforce_small is the slow method of Section 5: for every $x$ and every $y$ from $x$ up to the bound of Section 5 it computes $d=xy-n$, requires $d\gt0$ and $d\mid n(x+y)$, sets $z=n(x+y)/d$, requires $z\ge y$, checks the equation once more and adds the perimeter. run_checkpoints computes the sums for $K=5$ and $K=10$ with both methods and compares them with each other and with the stored values $140098$ and $3781786$ of Section 9; on a mismatch it prints both values, and the program stops with exit status $2$.

By Section 5, $M=nm\lt4n^2$, at most about $6.4\cdot10^{13}\lt2^{46}$ for $k\le1000$, and by Section 10 every perimeter is at most about $3.2\cdot10^{13}$, so all these numbers fit easily into the 64-bit type u64. The product u * u does not always fit: for a divisor $u\ge2^{32}$ it wraps around modulo $2^{64}$, and the test can then keep a divisor with $u^2\gt M$. A replay of the loops with exact 128-bit products shows that this happens for $930$ of the $573\,323\,938$ divisors generated for $k\le1000$, and that each of them then fails a divisibility test, so the printed total is not affected. The Java version writes the test as u > M / u, which is equivalent for positive integers and cannot overflow; Python's integers have no size limit.

The Python and Java versions follow the same steps without the checkpoints. Python has a single function solve with $K=1000$ fixed and the nested helpers factorize, merge_factors and gen_divisors; its sieve starts from list(range(n_max + 1)) and overwrites entries for the primes up to $\sqrt{16\cdot10^6}=4000$. Java has buildSpf, factorizeInt, mergeFactorizations, generateDivisors and solveRange, keeps the divisors in a List<Long> and the total in a long, and its solve adds the totals of consecutive blocks of $k$ like solve_fast.

Complexity and Verification

The table has $16\cdot10^6+1$ entries, $64$ MB of 32-bit integers in C++ and Java, and the sieve takes $O(L\log\log L)$ time for $L=16K^2$. ^[7] Since $x\le\lfloor2\sqrt3\,k\rfloor$, the main loops visit about $\sqrt3\,K^2$ pairs $(k,x)$, exactly $1\,733\,282$ for $K=1000$. Each pair costs a factorization of $m$, a few dozen divisions at most, and one step per divisor of $M$: $573\,323\,938$ divisors in total, about $331$ per pair on average and at most $9720$, for $k=924$ and $x=1986$. The number of divisors grows more slowly than any power of the number, ^[8] so the work is about $K^2$ times a slowly growing factor. The slow method of Section 5 would try $21\,210\,932\,262$ values of $y$, about $37$ times as many. On a desktop computer the C++ and Java programs take about a second and the Python version well under a minute. Python's list(range(n_max + 1)) first creates sixteen million integer objects, so its peak memory use is several hundred megabytes, far more than the $64$ MB table of C++ and Java.

The C++ program checks itself with run_checkpoints before the main computation. For this article, further checks that are not part of the published solutions were run:

A brute force that uses nothing but Heron's formula in the form of Section 2 went through all $60\,760\,114\,290$ integer triangles with perimeter at most $20604=2\cdot101\cdot102$, the bound of Section 10 for $k=5$. It found $5$, $18$, $45$, $45$ and $52$ triangles for $k=1,\dots,5$, with the sums of Section 9 and the total $140098$; none of them has an odd perimeter, and for each $k$ the largest perimeter is exactly $2(n+1)(n+2)$.
For every $k\le10$ the divisor method and the slow method of Section 5 produce the same triples.
The replay with 128-bit products for $k\le1000$ checked every triple against $xyz=n(x+y+z)$ and $y\le z$, found $3\,269\,538$ triangles, gave the same total as the published C++ program and confirmed that the total is below $2^{63}$, so the long of the Java version cannot overflow.
The C++, Java and Python programs print the same total.

Footnotes and References

The references below supply the background facts used above. The divisor parametrization of Section 6, the bounds of Sections 5 and 10 and the analysis of the two divisibility tests are derived explicitly in this article.

Project Euler 283 — Integer Sided Triangles with Integral Area/perimeter Ratio. The official statement: the triangle with sides 6, 8, 10 has perimeter and area 24, so its area/perimeter ratio is 1, and the triangle with sides 13, 14, 15 has perimeter 42 and area 84, so its ratio is 2; find the sum of the perimeters of all integer sided triangles whose area/perimeter ratio is a positive integer not exceeding 1000.
Wikipedia — Heron's formula. The area of a triangle with sides a, b, c and semiperimeter s is √(s(s − a)(s − b)(s − c)), or equivalently 4A = √((a + b + c)(−a + b + c)(a − b + c)(a + b − c)), the form without fractions used in Section 2.
Wikipedia — Incircle and excircles. The area is Δ = rs; the inradius is r = √((s − a)(s − b)(s − c)/s); and each vertex is at distance s − a, s − b or s − c from its two nearest touch points of the incircle.
Wikipedia — Equable shape. There are exactly five equable triangles with integer sides, (5, 12, 13), (6, 8, 10), (6, 25, 29), (7, 15, 20) and (9, 10, 17), as Whitworth and Biddle proved in 1904, and every equable triangle has inradius 2.
Wikipedia — Heronian triangle. Heronian triangles have integer sides and integer area; by Heron's formula they are the positive integer solutions of 16T² = (a + b + c)(a + b − c)(b + c − a)(c + a − b). Heron's own example was the triangle 13, 14, 15 with area 84.
Wikipedia — Euclid's lemma. In its general form: if an integer divides a product ab and is coprime with a, then it divides b. Section 8 applies it to x, which divides the product n(v + n).
Wikipedia — Sieve of Eratosthenes. The sieve marks the multiples of every prime p, starting from p² because the smaller multiples have already been marked, and runs in O(n log log n) time.
Wikipedia — Divisor function. A number p₁^a₁ ⋯ pᵣ^aᵣ has (a₁ + 1) ⋯ (aᵣ + 1) divisors; the number of divisors satisfies d(n) = o(n^ε) for every ε > 0, and the numbers up to x have about log x divisors on average.
C++ — Euler283.cpp. The linked immutable C++ revision contains build_spf, solve_range, solve_bruteforce_small and run_checkpoints with the stored values 140098 and 3781786. It is the implementation record for this article.
GCC manual — 128-bit Integers. GCC provides the types __int128 and unsigned __int128 on targets that have an integer mode wide enough to hold 128 bits.

Problem 283 source code

C++

#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>
#include <thread>
#include <utility>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;
using i64 = std::int64_t;

constexpr int kDefaultLimit = 1000;

struct Options {
    int k_limit = kDefaultLimit;
    unsigned threads = std::thread::hardware_concurrency();
    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_unsigned_after_prefix(const std::string& arg, const std::string& prefix, unsigned& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    unsigned parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10U + static_cast<unsigned>(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, "--k-limit=", options.k_limit)) {
            continue;
        }
        if (parse_unsigned_after_prefix(arg, "--threads=", options.threads)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    if (options.k_limit <= 0) {
        return false;
    }
    if (options.threads == 0) {
        options.threads = 1;
    }
    return true;
}

std::vector<int> build_spf(const int limit) {
    std::vector<int> spf(static_cast<std::size_t>(limit + 1), 0);
    if (limit >= 1) {
        spf[1] = 1;
    }
    for (int i = 2; i <= limit; ++i) {
        if (spf[static_cast<std::size_t>(i)] != 0) {
            continue;
        }
        spf[static_cast<std::size_t>(i)] = i;
        if (static_cast<i64>(i) * i > limit) {
            continue;
        }
        for (int j = i * i; j <= limit; j += i) {
            if (spf[static_cast<std::size_t>(j)] == 0) {
                spf[static_cast<std::size_t>(j)] = i;
            }
        }
    }
    return spf;
}

void factorize_int(int n, const std::vector<int>& spf, std::vector<std::pair<int, int>>& out) {
    out.clear();
    while (n > 1) {
        const int p = spf[static_cast<std::size_t>(n)];
        int cnt = 0;
        while (n % p == 0) {
            n /= p;
            ++cnt;
        }
        out.push_back({p, cnt});
    }
}

void merge_factorizations(const std::vector<std::pair<int, int>>& a,
                         const std::vector<std::pair<int, int>>& b,
                         std::vector<std::pair<int, int>>& out) {
    out.clear();
    std::size_t i = 0;
    std::size_t j = 0;

    while (i < a.size() || j < b.size()) {
        if (j == b.size() || (i < a.size() && a[i].first < b[j].first)) {
            out.push_back(a[i]);
            ++i;
        } else if (i == a.size() || b[j].first < a[i].first) {
            out.push_back(b[j]);
            ++j;
        } else {
            out.push_back({a[i].first, a[i].second + b[j].second});
            ++i;
            ++j;
        }
    }
}

void generate_divisors(const std::vector<std::pair<int, int>>& factors, std::vector<u64>& divisors) {
    divisors.clear();
    divisors.push_back(1);

    for (const auto& [p, e] : factors) {
        const std::size_t base = divisors.size();
        u64 mul = 1;
        for (int i = 1; i <= e; ++i) {
            mul *= static_cast<u64>(p);
            for (std::size_t j = 0; j < base; ++j) {
                divisors.push_back(divisors[j] * mul);
            }
        }
    }
}

std::string to_string_u128(u128 value) {
    if (value == 0) {
        return "0";
    }
    std::string s;
    while (value > 0) {
        const int digit = static_cast<int>(value % 10);
        s.push_back(static_cast<char>('0' + digit));
        value /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

u128 solve_range(const int k_start, const int k_end, const std::vector<int>& spf) {
    std::vector<std::pair<int, int>> factors_n;
    std::vector<std::pair<int, int>> factors_m;
    std::vector<std::pair<int, int>> factors_all;
    std::vector<u64> divisors;

    u128 total = 0;

    for (int k = k_start; k <= k_end; ++k) {
        const u64 r = static_cast<u64>(2 * k);
        const u64 n = r * r;
        factorize_int(static_cast<int>(n), spf, factors_n);

        const u64 x_max = static_cast<u64>(std::sqrt(static_cast<long double>(3 * n)));
        for (u64 x = 1; x <= x_max; ++x) {
            const u64 m = x * x + n;
            factorize_int(static_cast<int>(m), spf, factors_m);
            merge_factorizations(factors_n, factors_m, factors_all);
            generate_divisors(factors_all, divisors);

            const u64 M = n * m;
            for (u64 u : divisors) {
                if (u * u > M) {
                    continue;
                }

                const u64 v = M / u;
                if ((u + n) % x != 0 || (v + n) % x != 0) {
                    continue;
                }

                const u64 y = (u + n) / x;
                if (y < x) {
                    continue;
                }
                const u64 z = (v + n) / x;

                total += static_cast<u128>(2) * (x + y + z);
            }
        }
    }

    return total;
}

u128 solve_fast(const int k_limit, const std::vector<int>& spf, unsigned threads) {
    if (threads <= 1 || k_limit < 100) {
        return solve_range(1, k_limit, spf);
    }

    const unsigned use_threads = std::min<unsigned>(threads, static_cast<unsigned>(k_limit));
    const int chunk = (k_limit + static_cast<int>(use_threads) - 1) / static_cast<int>(use_threads);

    std::vector<std::thread> pool;
    std::vector<u128> partial(use_threads, 0);
    pool.reserve(use_threads);

    for (unsigned t = 0; t < use_threads; ++t) {
        const int start = static_cast<int>(t) * chunk + 1;
        const int end = std::min(k_limit, start + chunk - 1);
        if (start > end) {
            continue;
        }

        pool.emplace_back([&, start, end, t]() {
            partial[t] = solve_range(start, end, spf);
        });
    }

    for (auto& th : pool) {
        th.join();
    }

    u128 total = 0;
    for (u128 part : partial) {
        total += part;
    }
    return total;
}

u128 solve_bruteforce_small(const int k_limit) {
    u128 total = 0;

    for (int k = 1; k <= k_limit; ++k) {
        const u64 r = static_cast<u64>(2 * k);
        const u64 n = r * r;

        const u64 x_max = static_cast<u64>(std::sqrt(static_cast<long double>(3 * n)));
        for (u64 x = 1; x <= x_max; ++x) {
            const long double disc = std::sqrt(static_cast<long double>(n) * (static_cast<long double>(n) + x * x));
            const u64 y_max = static_cast<u64>((n + disc) / x + 1e-12L);

            for (u64 y = x; y <= y_max; ++y) {
                const i64 d = static_cast<i64>(x * y) - static_cast<i64>(n);
                if (d <= 0) {
                    continue;
                }

                const u64 num = n * (x + y);
                if (num % static_cast<u64>(d) != 0) {
                    continue;
                }

                const u64 z = num / static_cast<u64>(d);
                if (z < y) {
                    continue;
                }

                if (x * y * z != n * (x + y + z)) {
                    continue;
                }

                total += static_cast<u128>(2) * (x + y + z);
            }
        }
    }

    return total;
}

bool run_checkpoints(const std::vector<int>& spf) {
    const u128 fast5 = solve_fast(5, spf, 1);
    const u128 brute5 = solve_bruteforce_small(5);
    if (fast5 != brute5 || fast5 != static_cast<u128>(140098ULL)) {
        std::cerr << "Checkpoint failed for k<=5: fast=" << to_string_u128(fast5)
                  << ", brute=" << to_string_u128(brute5) << '\n';
        return false;
    }

    const u128 fast10 = solve_fast(10, spf, 1);
    const u128 brute10 = solve_bruteforce_small(10);
    if (fast10 != brute10 || fast10 != static_cast<u128>(3781786ULL)) {
        std::cerr << "Checkpoint failed for k<=10: fast=" << to_string_u128(fast10)
                  << ", brute=" << to_string_u128(brute10) << '\n';
        return false;
    }

    return true;
}

}  // namespace

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

    // max(x^2 + n) with n=(2k)^2 and x<=sqrt(3n) is at most 4*n
    const int n_max = 4 * (2 * options.k_limit) * (2 * options.k_limit);
    const std::vector<int> spf = build_spf(n_max);

    if (options.run_checkpoints && !run_checkpoints(spf)) {
        return 2;
    }

    const u128 answer = solve_fast(options.k_limit, spf, options.threads);
    std::cout << to_string_u128(answer) << '\n';
    return 0;
}

Python

import math

def solve():
    k_limit = 1000
    n_max = 4 * (2 * k_limit) ** 2

    # Build smallest prime factor
    spf = list(range(n_max + 1))
    spf[0] = 0
    if n_max >= 1:
        spf[1] = 1
    for i in range(2, math.isqrt(n_max) + 1):
        if spf[i] == i:
            for j in range(i*i, n_max + 1, i):
                if spf[j] == j:
                    spf[j] = i

    def factorize(n):
        factors = []
        while n > 1:
            p = spf[n]
            e = 0
            while n % p == 0:
                n //= p
                e += 1
            factors.append((p, e))
        return factors

    def merge_factors(a, b):
        result = []
        i = j = 0
        while i < len(a) or j < len(b):
            if j == len(b) or (i < len(a) and a[i][0] < b[j][0]):
                result.append(a[i]); i += 1
            elif i == len(a) or b[j][0] < a[i][0]:
                result.append(b[j]); j += 1
            else:
                result.append((a[i][0], a[i][1] + b[j][1])); i += 1; j += 1
        return result

    def gen_divisors(factors):
        divs = [1]
        for p, e in factors:
            base = len(divs)
            mul = 1
            for _ in range(e):
                mul *= p
                for j in range(base):
                    divs.append(divs[j] * mul)
        return divs

    total = 0
    for k in range(1, k_limit + 1):
        r = 2 * k
        n = r * r
        factors_n = factorize(n)

        x_max = int(math.sqrt(3 * n))
        for x in range(1, x_max + 1):
            m = x * x + n
            factors_m = factorize(m)
            factors_all = merge_factors(factors_n, factors_m)
            divs = gen_divisors(factors_all)

            M = n * m
            for u in divs:
                if u * u > M:
                    continue
                v = M // u
                if (u + n) % x != 0 or (v + n) % x != 0:
                    continue
                y = (u + n) // x
                if y < x:
                    continue
                z = (v + n) // x
                total += 2 * (x + y + z)

    return str(total)

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

Java

import java.util.*;
import java.util.concurrent.*;

public class Euler283 {
    static class Pair {
        int first, second;

        Pair(int f, int s) {
            first = f;
            second = s;
        }
    }

    static int[] buildSpf(int limit) {
        int[] spf = new int[limit + 1];
        if (limit >= 1)
            spf[1] = 1;
        for (int i = 2; i <= limit; ++i) {
            if (spf[i] != 0)
                continue;
            spf[i] = i;
            if ((long) i * i > limit)
                continue;
            for (int j = i * i; j <= limit; j += i) {
                if (spf[j] == 0)
                    spf[j] = i;
            }
        }
        return spf;
    }

    static void factorizeInt(int n, int[] spf, List<Pair> out) {
        out.clear();
        while (n > 1) {
            int p = spf[n];
            int cnt = 0;
            while (n % p == 0) {
                n /= p;
                ++cnt;
            }
            out.add(new Pair(p, cnt));
        }
    }

    static void mergeFactorizations(List<Pair> a, List<Pair> b, List<Pair> out) {
        out.clear();
        int i = 0, j = 0;
        while (i < a.size() || j < b.size()) {
            if (j == b.size() || (i < a.size() && a.get(i).first < b.get(j).first)) {
                out.add(a.get(i));
                ++i;
            } else if (i == a.size() || b.get(j).first < a.get(i).first) {
                out.add(b.get(j));
                ++j;
            } else {
                out.add(new Pair(a.get(i).first, a.get(i).second + b.get(j).second));
                ++i;
                ++j;
            }
        }
    }

    static void generateDivisors(List<Pair> factors, List<Long> divisors) {
        divisors.clear();
        divisors.add(1L);
        for (Pair p : factors) {
            int base = divisors.size();
            long mul = 1;
            for (int i = 1; i <= p.second; ++i) {
                mul *= p.first;
                for (int j = 0; j < base; ++j) {
                    divisors.add(divisors.get(j) * mul);
                }
            }
        }
    }

    static long solveRange(int kStart, int kEnd, int[] spf) {
        List<Pair> factorsN = new ArrayList<>();
        List<Pair> factorsM = new ArrayList<>();
        List<Pair> factorsAll = new ArrayList<>();
        List<Long> divisors = new ArrayList<>();

        long total = 0;

        for (int k = kStart; k <= kEnd; ++k) {
            long r = 2L * k;
            long n = r * r;
            factorizeInt((int) n, spf, factorsN);

            long xMax = (long) Math.sqrt(3L * n);
            for (long x = 1; x <= xMax; ++x) {
                long m = x * x + n;
                factorizeInt((int) m, spf, factorsM);
                mergeFactorizations(factorsN, factorsM, factorsAll);
                generateDivisors(factorsAll, divisors);

                long M = n * m;
                for (long u : divisors) {
                    if (u > M / u)
                        continue;
                    long v = M / u;
                    if ((u + n) % x != 0 || (v + n) % x != 0)
                        continue;

                    long y = (u + n) / x;
                    if (y < x)
                        continue;
                    long z = (v + n) / x;

                    total += 2L * (x + y + z);
                }
            }
        }
        return total;
    }

    public static String solve() {
        int kLimit = 1000;
        int threads = Math.max(1, Runtime.getRuntime().availableProcessors());

        int nMax = 4 * (2 * kLimit) * (2 * kLimit);
        int[] spf = buildSpf(nMax);

        if (threads <= 1 || kLimit < 100) {
            long ans = solveRange(1, kLimit, spf);
            return String.valueOf(ans);
        }

        int useThreads = Math.min(threads, kLimit);
        int chunk = (kLimit + useThreads - 1) / useThreads;

        ExecutorService executor = Executors.newFixedThreadPool(useThreads);
        List<Future<Long>> futures = new ArrayList<>();

        for (int t = 0; t < useThreads; ++t) {
            final int start = t * chunk + 1;
            final int end = Math.min(kLimit, start + chunk - 1);
            if (start > end)
                continue;

            futures.add(executor.submit(() -> solveRange(start, end, spf)));
        }

        long total = 0;
        for (Future<Long> f : futures) {
            try {
                total += f.get();
            } catch (Exception e) {
                e.printStackTrace();
            }
        }
        executor.shutdown();

        return String.valueOf(total);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}