Problem 397: Triangle on Parabola

Project Euler Problem 397 Solution

EulerSolve provides an optimized solution for Project Euler Problem 397, Triangle on Parabola, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each integer $k \ge 1$, consider the parabola $y=\dfrac{x^2}{k}$. We choose three distinct integers $a \lt b \lt c$ with $|a|,|b|,|c|\le X$, form the points $A=\left(a,\dfrac{a^2}{k}\right)$, $B=\left(b,\dfrac{b^2}{k}\right)$, $C=\left(c,\dfrac{c^2}{k}\right)$, and ask whether triangle $ABC$ has at least one $45^\circ$ angle. Let $T_k(X)$ be that count for fixed $k$. The program computes $$F(K,X)=\sum_{k=1}^{K} T_k(X).$$ The key reduction is that chord slopes on this parabola depend only on pair sums such as $a+b$, so the geometry becomes a divisor problem for $2k^2$ plus interval counting for the remaining coordinate. Mathematical Approach Step 1: Replace Geometry by Pair Sums Define $$x=a+b,\qquad y=a+c,\qquad z=b+c,$$ so $x \lt y \lt z$ because $a \lt b \lt c$. The slope of a chord between two points on $y=\dfrac{x^2}{k}$ is $$m_{uv}=\frac{v^2-u^2}{k(v-u)}=\frac{u+v}{k}.$$ Therefore the three side slopes are exactly $$m_{AB}=\frac{x}{k},\qquad m_{AC}=\frac{y}{k},\qquad m_{BC}=\frac{z}{k}.$$ This is the crucial simplification: once $(x,y,z)$ are known, the angle conditions no longer involve quadratic expressions in $a,b,c$. Step 2: Convert the $45^\circ$ Condition into Diophantine Equations At vertex $A$, the two scaled direction vectors are proportional to $(k,x)$ and $(k,y)$....

Detailed mathematical approach

Problem Summary

On the parabola $y=x^2/k$, with $k\ge1$ an integer, three integers $a\lt b\lt c$ pick out the points $A(a,a^2/k)$, $B(b,b^2/k)$ and $C(c,c^2/k)$. Let $F(K,X)$ be the number of integer quadruplets $(k,a,b,c)$ with $1\le k\le K$ and $-X\le a\lt b\lt c\le X$ for which at least one angle of the triangle $ABC$ measures $45^\circ$. The statement gives the checkpoints $F(1,10)=41$ and $F(10,100)=12492$, and asks for $F(10^6,10^9)$. ^[1]

The whole solution in one sentence: a chord of this parabola has slope $(u+v)/k$, so a $45^\circ$ angle turns into a divisor of $2k^2$, and every divisor then costs no more than the length of one interval of integers.

Mathematical Approach

The geometry is first rewritten in the three pair sums $a+b$, $a+c$ and $b+c$, because the slopes of the three sides depend on nothing else. Each $45^\circ$ condition then becomes one equation with a product of two integers on the left and $2k^2$ on the right, so its solutions are read off the divisors of $2k^2$. One interval length counts the quadruplets of each solution, and triangles with two $45^\circ$ angles are subtracted at the end.

1. Chord slopes are pair sums

Take $k=1$ with $a=-1$, $b=1$, $c=2$, so that $A=(-1,1)$, $B=(1,1)$ and $C=(2,4)$. The side $AB$ is horizontal, and $AC$ climbs $3$ units while moving $3$ units to the right, so the angle at $A$ is $45^\circ$. The two slopes, $0$ and $1$, are exactly the sums $a+b$ and $a+c$. That is the general rule: for $u\ne v$ the chord of $y=x^2/k$ through the points at $u$ and $v$ has slope

$$\begin{aligned}\frac{v^2/k-u^2/k}{v-u}&=\frac{(v-u)(v+u)}{k(v-u)}\\&=\frac{u+v}{k}.\end{aligned}$$

Only the sum of the two abscissae matters. The tangent at $u$ has slope $2u/k$, ^[2] so a chord slope is the mean of the tangent slopes at its two ends. Write $x=a+b$, $y=a+c$ and $z=b+c$; then $x\lt y\lt z$ follows from $a\lt b\lt c$, and the sides $AB$, $AC$, $BC$ have the slopes $x/k$, $y/k$, $z/k$. The three points are never collinear, since a line $y=mx+n$ meets the parabola where $x^2-kmx-kn=0$, a quadratic with at most two roots; every triple therefore gives a genuine triangle.

2. The three angles, read off the three slopes

Give the three sides their direction angles

$$\begin{aligned}\alpha&=\arctan\frac{x}{k},\\\beta&=\arctan\frac{y}{k},\\\delta&=\arctan\frac{z}{k},\end{aligned}$$

all of them in $(-90^\circ,90^\circ)$ and, as $\arctan$ is increasing, with $\alpha\lt\beta\lt\delta$. The step from $A$ to $B$ is $(b-a)\bigl(1,x/k\bigr)$ with $b-a\gt0$, so the ray $AB$ points in the direction $\alpha$; likewise $AC$ points along $\beta$ and $BC$ along $\delta$, while the reversed rays $BA$, $CA$, $CB$ point along $\alpha+180^\circ$, $\beta+180^\circ$ and $\delta+180^\circ$. Each interior angle is the difference of the directions of its two rays:

$$\begin{aligned}\angle A&=\beta-\alpha,\\\angle C&=\delta-\beta,\\\angle B&=180^\circ-(\delta-\alpha),\end{aligned}$$

three positive numbers that add up to $180^\circ$. In the range $0^\circ\lt\theta\lt180^\circ$ the equation $\tan\theta=1$ has the single solution $\theta=45^\circ$, so a $45^\circ$ angle at $A$, at $C$ or at $B$ says exactly that

$$\begin{aligned}\tan(\beta-\alpha)&=1,\\\tan(\delta-\beta)&=1,\\\tan(\delta-\alpha)&=-1.\end{aligned}$$

3. Three equations in the pair sums

The tangent of a difference of angles ^[3] turns the first of them into a relation between the sums,

$$\begin{aligned}\tan(\beta-\alpha)&=\frac{\tan\beta-\tan\alpha}{1+\tan\alpha\tan\beta}\\&=\frac{k(y-x)}{k^2+xy},\end{aligned}$$

and the same substitution works for the other two pairs of directions. Clearing the denominators leaves

$$\begin{aligned}k(y-x)&=k^2+xy,\\k(z-y)&=k^2+yz,\\k(z-x)&=-(k^2+xz).\end{aligned}$$

Nothing is lost when a denominator vanishes: $1+\tan\alpha\tan\beta=0$ means the two directions differ by exactly $90^\circ$, neither $45^\circ$ nor $135^\circ$. All three equations are bilinear, so each factorises by completing the product:

$$\begin{aligned}(k-x)(k+y)&=2k^2,\\(k-y)(k+z)&=2k^2,\\(-k-x)(z-k)&=2k^2.\end{aligned}$$

Expanding the first product gives $k^2+ky-kx-xy$, and setting that equal to $2k^2$ returns $k(y-x)=k^2+xy$. Each equation describes a rectangular hyperbola, in the first case with the lines $x=k$ and $y=-k$ as asymptotes, ^[4] so counting its lattice points means counting factorisations of $2k^2$. The first two equations are one and the same, applied to the consecutive pair $(x,y)$ and to the consecutive pair $(y,z)$; that is why the program keeps one list of pairs for the angles at $A$ and at $C$, and a second list for the angle at $B$.

4. The divisor form, and why both factors are positive

Claim. Let $k\ge1$ and let $a\lt b\lt c$ be integers. The angle at the named vertex is $45^\circ$ if and only if the $d$ named with it is a positive divisor of $2k^2$ and the remaining sum has the stated value:

at $A$: $d=k-x$, and then $y=2k^2/d-k$;
at $C$: $d=k-y$, and then $z=2k^2/d-k$;
at $B$: $d=-k-x$, and then $z=2k^2/d+k$.

Conversely every positive divisor $d$ of $2k^2$ delivers one admissible pair in each family: $(k-d,\ 2k^2/d-k)$ for the first two conditions and $(-k-d,\ 2k^2/d+k)$ for the third.

Proof. Take the angle at $A$, which by Section 3 is $45^\circ$ exactly when $(k-x)(k+y)=2k^2$. That product is positive, so its two factors have the same sign. Were both negative, say $k-x=-f$ and $k+y=-g$ with $f,g\ge1$, then $x=k+f$ and $y=-k-g$, hence $x\gt y$, contradicting $x\lt y$. So both factors are positive, $d=k-x$ is a positive divisor of $2k^2$, and $k+y=2k^2/d$, as claimed. Conversely, for a given $d$ the pair $x=k-d$, $y=2k^2/d-k$ satisfies the equation, and $x\lt y$ because

$$\begin{aligned}y-x&=\frac{2k^2-2kd+d^2}{d}\\&=\frac{k^2+(k-d)^2}{d}\gt0.\end{aligned}$$

The angle at $C$ is the same statement with $(x,y)$ replaced by $(y,z)$. For the angle at $B$ the equation is $(-k-x)(z-k)=2k^2$. Were both factors negative, write $-k-x=-f$ and $z-k=-g$ with $f,g\ge1$ and $fg=2k^2$; then $x=f-k$ and $z=k-g$, so $x\lt z$ would force $f+g\lt2k$, while $f+g\ge2\sqrt{fg}=2k\sqrt2\gt2k$. Conversely $x=-k-d$ is negative and $z=2k^2/d+k$ is positive, so $x\lt z$ is automatic. $\square$

The divisors of $2k^2$ are therefore the whole search space for one value of $k$. If $k=2^em$ with $m$ odd and $m=\prod p^{e_p}$, then $2k^2=2^{2e+1}\prod p^{2e_p}$, and the product formula for the divisor function ^[5] counts

$$d(2k^2)=(2e+2)\prod(2e_p+1)$$

of them.

5. Which divisors can occur at all

Each of $x,y,z$ is a sum of two numbers in $[-X,X]$, so all three lie in $[-2X,2X]$. For the first family this bounds $d$ from both sides: $x=k-d\ge-2X$ gives $d\le k+2X$, and $y=2k^2/d-k\le2X$ gives $d\ge2k^2/(2X+k)$. The two remaining inequalities are free whenever $k\le2X$, which holds for the target parameters. So only the divisors with

$$\frac{2k^2}{2X+k}\le d\le k+2X$$

survive. For the $B$-family the same argument gives $2k^2/(2X-k)\le d\le2X-k$ when $k\lt2X$, and no pair at all when $k\ge2X$, since then $x=-k-d\lt-2X$. As long as $2k^2+k\le2X$, that is up to $k=31622$ for $X=10^9$, no divisor is lost in either family; beyond that the extreme ones drop out, a very small $d$ pushing the partner sum past $2X$ and a very large $d$ pushing $x$ below $-2X$. At $K=10^6$ and $X=10^9$, $155\,556\,242$ of the $202\,789\,372$ candidate pairs survive, just under $77\%$.

6. Counting the free coordinate

A surviving pair fixes two of the three sums, and the third point is then counted, not searched. For a $45^\circ$ angle at $A$ the pair is $(x,y)=(a+b,a+c)$, so $b=x-a$ and $c=y-a$, and every admissible integer $a$ gives exactly one quadruplet. The coordinate bounds $-X\le a,b,c\le X$ and the ordering $a\lt b\lt c$ amount to $a\ge-X$, $a\le x+X$, $a\ge y-X$ and $2a\lt x$: the condition $b\lt c$ is $x\lt y$ itself, and $a\le X$, $a\ge x-X$, $a\le y+X$ follow from the four because $x\lt y\le2X$. For integers $2a\lt x$ means $a\le\lfloor(x-1)/2\rfloor$. So $a$ runs through an interval $L\le a\le U$ with

$$\begin{aligned}L&=\max(-X,\,y-X),\\U&=\min\Bigl(X,\,x+X,\Bigl\lfloor\tfrac{x-1}{2}\Bigr\rfloor\Bigr),\end{aligned}$$

and the number of quadruplets is $\max(0,\,U-L+1)$, with both endpoints included. For a $45^\circ$ angle at $B$ the free coordinate is $b$, with $a=x-b$ and $c=z-b$, and $a\lt b\lt c$ becomes $x\lt2b\lt z$; for one at $C$ it is $c$, with $a=y-c$ and $b=z-c$, where $a\lt b$ is $y\lt z$ and $b\lt c$ becomes $z\lt2c$. The same bookkeeping gives

$$\begin{aligned}L&=\max\Bigl(-X,\,z-X,\Bigl\lfloor\tfrac{x}{2}\Bigr\rfloor+1\Bigr),\\U&=\min\Bigl(X,\,x+X,\Bigl\lfloor\tfrac{z-1}{2}\Bigr\rfloor\Bigr)\end{aligned}$$

for $b$, and

$$\begin{aligned}L&=\max\Bigl(-X,\,z-X,\Bigl\lfloor\tfrac{z}{2}\Bigr\rfloor+1\Bigr),\\U&=\min(X,\,y+X)\end{aligned}$$

for $c$.

7. Triangles counted twice

For one $k$, let $T_k(X)$ be the number of triples $a\lt b\lt c$ in $[-X,X]$ whose triangle has at least one $45^\circ$ angle, let $N_A$, $N_B$, $N_C$ count the triples with such an angle at the named vertex, and let $N_{AB}$, $N_{AC}$, $N_{BC}$ count those with $45^\circ$ at both named vertices. A triangle cannot have three $45^\circ$ angles, because $3\cdot45^\circ=135^\circ\ne180^\circ$, so the triple term of inclusion–exclusion ^[6] vanishes:

$$\begin{aligned}T_k(X)&=N_A+N_B+N_C\\&\quad-N_{AB}-N_{AC}-N_{BC}.\end{aligned}$$

Two $45^\circ$ angles are possible, and then the third angle is $90^\circ$. Such a triangle shows up in two families, which share one pair sum: the angles at $A$ and $B$ share $x$, those at $A$ and $C$ share $y$, those at $B$ and $C$ share $z$. Inside one family the shared sum already determines its partner, since $d$ can be recovered from it, so each value occurs at most once per list, and one linear merge of two sorted lists finds every coincidence.

A coincidence delivers all three sums $x\lt y\lt z$, from which the triple is recovered by

$$\begin{aligned}a&=\frac{x+y-z}{2},\\b&=\frac{x+z-y}{2},\\c&=\frac{y+z-x}{2},\end{aligned}$$

an integer triple only if $x+y+z$ is even. Then $a\lt b\lt c$ is automatic, because $b-a=z-y\gt0$ and $c-b=y-x\gt0$, so only $-X\le a$ and $c\le X$ remain to be tested. Each coincidence contributes at most one quadruplet, the three corrections are exact, and the problem asks for $F(K,X)=\sum_{k\le K}T_k(X)$.

8. A complete example: $F(1,10)=41$

Let $k=1$, so that $2k^2=2$ has the two divisors $d=1$ and $d=2$. The claim of Section 4 gives six pairs:

$d=1$: $(x,y)=(0,1)$ at $A$, $(y,z)=(0,1)$ at $C$, $(x,z)=(-2,3)$ at $B$;
$d=2$: $(x,y)=(-1,0)$ at $A$, $(y,z)=(-1,0)$ at $C$, $(x,z)=(-3,2)$ at $B$.

With $X=10$ the intervals of Section 6 contain $9$ and $10$ values for the two $A$-pairs, $10$ and $9$ for the two $C$-pairs and $2$ and $2$ for the two $B$-pairs, so $N_A=19$, $N_C=19$ and $N_B=4$. The nine triples of the first $A$-pair are $(a,-a,1-a)$ for $a=-1,-2,\dots,-9$, that is $(-1,1,2)$ — the triangle of Section 1 — up to $(-9,9,10)$; those of the first $C$-pair are $(-c,1-c,c)$ for $c=1,\dots,10$, and the two $B$-pairs give $(-2,0,3)$, $(-3,1,2)$, $(-2,-1,3)$ and $(-3,0,2)$.

The sum $19+4+19=42$ exceeds the checkpoint by one. Comparing the shared sums finds the reason: no $x$ of an $A$-pair equals an $x$ of a $B$-pair, and no $z$ of a $B$-pair equals a $z$ of a $C$-pair, but the $A$-pair of $d=2$ and the $C$-pair of $d=1$ share $y=0$. Together they give $(x,y,z)=(-1,0,1)$, whose sum is even, and the recovery formula returns $(a,b,c)=(-1,0,1)$: the triangle $(-1,1)$, $(0,0)$, $(1,1)$ on $y=x^2$, with $45^\circ$ at $A$ and at $C$ and the right angle at $B$. It stands in both lists above, and subtracting it once gives $F(1,10)=42-1=41$, the value in the statement.

9. A symmetry that checks the bookkeeping

Reflecting the picture in the $y$-axis maps the parabola to itself and sends $(a,b,c)$ to $(-c,-b,-a)$, which is again increasing; the sums become $(-z,-y,-x)$, and the vertices $A$ and $C$ trade places while $B$ stays in the middle. This bijection of the counted quadruplets forces $N_A=N_C$ for every $k$ and every $X$, and makes the $A$-$B$ coincidences correspond one to one to the $B$-$C$ coincidences. For $k=1$ above this reads $19=19$. Summed over $1\le k\le10$ with $X=100$ the parts are $N_A=N_C=5462$ and $N_B=1588$, with $5$, $10$ and $5$ coincidences of the three kinds, so that the total $5462+1588+5462$ less those $20$ gives $12492$, the second checkpoint.

How the Code Works

The C++ program is the reference implementation. ^[10] parse_arguments reads --k, --x and --skip-checkpoints, defaulting to the parameters of the problem, and main builds one table with build_spf: a sieve that stores, for every integer up to the largest $k$ needed, the smallest prime dividing it. ^[7] With that table factor_2k_squared factorises $2k^2$ in a few steps, emitting the prime $2$ with exponent $2e+1$ for the exponent $e$ of $2$ in $k$ and then peeling off the odd primes of $k$ with doubled exponents, exactly as in Section 4, and collect_divisors_recursive turns that factor list into all divisors of $2k^2$.

build_pairs_for_k implements Sections 4 and 5. For every divisor $d$ it forms both signs $p=\pm d$ and $q=-2k^2/p$, and from them the candidate $(p+k,\,q-k)$ for the $A$- and $C$-family and $(p-k,\,q+k)$ for the $B$-family, keeping a candidate only when its first entry is smaller than its second and both lie in $[-2X,2X]$. By the proof in Section 4 only the sign $p=-d$ ever passes the ordering test, so half of the candidates are discarded at once, and distinct divisors give distinct pairs, so the std::unique pass after the sort never removes anything; the sort itself is needed because the merges require sorted lists. Finally plus_in and minus_in hold the same pairs with their entries swapped and sorted again, as the merges on the shared sums $y$ and $z$ need.

count_angle_a_for_pair, count_angle_b_for_pair and count_angle_c_for_pair are the three interval formulas of Section 6, and count_interval returns $0$ for an empty interval. Because the sums are frequently negative, they call a hand-written floor_div: C++ and Java truncate integer division towards zero, so $-3/2$ would come out as $-1$ instead of the required $-2$. solve_range loops over $k$, accumulates the three counts over the pair lists and runs the three merges of Section 7, and valid_triangle_from_sums performs the parity test, the recovery of $(a,b,c)$ — its division by $2$ is exact because the parity is checked first — and the range test; every match lowers the subtotal by one. The subtotals are accumulated in unsigned __int128 and printed by to_string_u128, since the stream operators do not handle that type. ^[8] Unless --skip-checkpoints is given, main first recomputes the three entries of kCheckpoints and stops on any mismatch.

In the Java port, solve fixes the two target parameters, builds the same smallest-prime-factor table and adds up the long results of solveForK, which repeats the factorisation, the interval counts and the overlap merges; getDivisors generates the divisors, Factor stores a prime with its exponent, and Pair implements Comparable so that Collections.sort orders lexicographically before a scan removes duplicates. That port carries no checkpoint table of its own. The Python file is not a second implementation of the mathematics but a wrapper: it expects the C++ source in the repository's solutionsCpp directory, compiles it with clang++ or g++ when the cached binary is missing or older, passes --skip-checkpoints, and returns the value behind the last = of the output.

Complexity and Verification

The sieve costs $O(K\log\log K)$ time ^[7] and one 32-bit entry per integer, $4$ MB at $K=10^6$. For one $k$, factorising $2k^2$ takes $O(\log k)$ steps, generating the divisors and the pairs is linear in $D_k=d(2k^2)$, the sorts dominate with $O(D_k\log D_k)$, and the three merges are linear again, so the whole computation runs in

$$O\Bigl(K\log\log K+\sum_{k\le K}D_k\log D_k\Bigr)$$

time, with the sieve plus a few lists of at most $D_k$ pairs in memory. The divisor sums are small: $\sum_{k\le10^6}d(2k^2)=101\,394\,686$, an average of $101.39$ and a maximum of $4050$, reached at $k=720720$, which is $2^4\cdot3^2\cdot5\cdot7\cdot11\cdot13$, where $2k^2=2^9\cdot3^4\cdot5^2\cdot7^2\cdot11^2\cdot13^2$ has $10\cdot5\cdot3^4=4050$ divisors. The order of growth is $O(K\log^2K)$: with $k=2^em$ and $m$ odd the product formula gives $d(2k^2)/d(k^2)=(2e+2)/(2e+1)$, hence $d(k^2)\le d(2k^2)\le2d(k^2)$, and the partial sums of $d(k^2)$, which reach $73\,858\,790$ at $K=10^6$, grow like $(3/\pi^2)K\log^2K$ with slowly decaying further terms. ^[9] Compiled with -O2, one pass over $1\le k\le10^6$ at $X=10^9$ was measured at about three seconds on a current laptop.

The result fits in 64 bits, which is what the Java port relies on. Each family contributes at most one pair per divisor, and each pair at most $2X+1$ quadruplets, because the free coordinate stays inside $[-X,X]$. Hence

$$\begin{aligned}F(K,X)&\le3(2X+1)\sum_{k\le K}d(2k^2)\\&\lt6.1\cdot10^{17},\end{aligned}$$

about fifteen times below the largest long, which is $2^{63}-1$. ^[8] The C++ accumulators are wider than they need to be.

The packaged checks are the three checkpoints $F(1,10)=41$, $F(10,100)=12492$ and $F(1000,10^4)=329\,864\,366$, verified here for the C++ program and, with the two parameters changed, for the Java port. The method was also checked independently for this article. A separate brute force works on the original geometry: it scales the parabola by $k$, a similarity that leaves all angles unchanged, so that the points become the lattice points $(ku,u^2)$, and at each vertex it compares the cross product and the dot product of the two integer edge vectors, a $45^\circ$ angle being the case of a positive dot product equal to the absolute cross product. It reproduces both statement checkpoints, and agrees with the fast method for every $K\le6$ and $2\le X\le12$, as well as at $F(20,60)=12930$ and $F(3,200)=4783$. The divisor description of Section 4 was tested against the same geometry for all $k\le14$ and $|a|,|b|,|c|\le25$, and holds for exactly those triples whose angles the brute force finds. No triple with three $45^\circ$ angles turned up there, as proved in Section 7; the counts of the double cases match the merges, $5$, $10$ and $5$ at $K=10$ and $X=100$; and the symmetry $N_A=N_C$ of Section 9 was confirmed for $k\lt40$ at several $X$ up to $10^6$. These comparisons are validation for this article, not tests inside the published sources.

Footnotes and References

The references below supply the background facts used above. The reduction to pair sums, the three factorised equations, the divisor description with its sign argument, the interval formulas and the correction for double $45^\circ$ angles are derived in this article.

Project Euler 397 — Triangle on Parabola. The official statement puts the points A(a, a²/k), B(b, b²/k) and C(c, c²/k) on the parabola y = x²/k, defines F(K, X) as the number of integer quadruplets (k, a, b, c) with at least one 45-degree angle in the triangle ABC, 1 ≤ k ≤ K and −X ≤ a < b < c ≤ X, gives F(1, 10) = 41 and F(10, 100) = 12492, and asks for F(10⁶, 10⁹). The two examples test an implementation without disclosing the requested value.
Wikipedia — Parabola. For the parabola as the graph of a function the article gives the tangent of y = ax² at (x₀, ax₀²) as y = 2ax₀x − ax₀², that is with slope 2ax₀. With a = 1/k this is the tangent slope 2u/k of Section 1.
Wikipedia — List of trigonometric identities. Among the angle sum and difference identities the article states tan(α ± β) = (tan α ± tan β)/(1 ∓ tan α tan β); the lower signs give the difference formula used in Section 3.
Wikipedia — Hyperbola. The article states that the graph of y = A/x with A > 0 is a rectangular hyperbola in the first and third quadrants with the coordinate axes as asymptotes. The equations of Section 3 are this curve after shifting both coordinates.
Wikipedia — Divisor function. The article gives σ₀(n) = ∏(aᵢ + 1) for n = ∏ pᵢ^aᵢ and notes that σ₀ is multiplicative; that product formula is what Section 4 applies to 2k².
Wikipedia — Inclusion–exclusion principle. The article states the three-set case |A ∪ B ∪ C| = |A| + |B| + |C| − |A ∩ B| − |A ∩ C| − |B ∩ C| + |A ∩ B ∩ C|, which Section 7 uses with an empty triple intersection.
Wikipedia — Sieve of Eratosthenes. The article describes marking the multiples of each prime in turn and states the time complexity O(n log log n) used in the complexity section. The implementation keeps, instead of a single mark, the prime that first marks a composite, which makes every integer below the limit factorisable in O(log n) steps.
GCC manual — 128-bit integers documents __int128 and unsigned __int128 as an extension available wherever a 128-bit integer mode exists; the C++ accumulators use the unsigned form. The Java Language Specification, §4.2.1, fixes the range of long as −9223372036854775808 to 9223372036854775807, the bound used above.
OEIS A048691 — number of divisors of n². The entry defines a(n) = d(n²), gives the prime-power values a(pᵉ) = 2e + 1 and the Dirichlet generating function ζ(s)³/ζ(2s), and records an asymptotic for its partial sums whose leading term is (3/π²)n log²n, the growth that the complexity section carries over to the sum of d(2k²).
C++ — Euler397.cpp. The linked immutable revision contains build_spf, factor_2k_squared, build_pairs_for_k, the three interval counters, valid_triangle_from_sums and the kCheckpoints table. It is the implementation record for this article; the independent comparisons above are not part of it.

Problem 397 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <utility>
#include <vector>

namespace {

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

constexpr u64 kDefaultK = 1'000'000ULL;
constexpr i64 kDefaultX = 1'000'000'000LL;

struct Checkpoint {
    u64 k = 0ULL;
    i64 x = 0LL;
    u64 expected = 0ULL;
};

constexpr Checkpoint kCheckpoints[] = {
    {1ULL, 10LL, 41ULL},
    {10ULL, 100LL, 12'492ULL},
    {1'000ULL, 10'000LL, 329'864'366ULL},
};

struct Options {
    u64 k = kDefaultK;
    i64 x = kDefaultX;
    bool run_checkpoints = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_i64_after_prefix(const std::string& arg, const char* prefix, i64& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }

    if (parsed > static_cast<u64>(std::numeric_limits<i64>::max())) {
        return false;
    }

    value = static_cast<i64>(parsed);
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
        return false;
    }
    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(const 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 (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }

        u64 parsed_u64 = 0ULL;
        if (parse_u64_after_prefix(arg, "--k=", parsed_u64)) {
            options.k = parsed_u64;
            continue;
        }

        i64 parsed_i64 = 0LL;
        if (parse_i64_after_prefix(arg, "--x=", parsed_i64)) {
            options.x = parsed_i64;
            continue;
        }

        unsigned parsed_unsigned = 0U;
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
            options.requested_threads = parsed_unsigned;
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    if (options.k == 0ULL) {
        std::cerr << "--k must be positive.\n";
        return false;
    }
    if (options.x <= 0LL) {
        std::cerr << "--x must be positive.\n";
        return false;
    }

    return true;
}

unsigned choose_thread_count(const bool allow_multithreading,
                             const unsigned requested_threads,
                             const std::size_t workload_units) {
    constexpr std::size_t kMinUnitsForParallel = 1'000ULL;

    if (!allow_multithreading || workload_units < 2ULL || workload_units < kMinUnitsForParallel) {
        return 1U;
    }

    unsigned threads = requested_threads;
    if (threads == 0U) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0U) {
            threads = 1U;
        }
    }

    return std::max(1U, std::min<unsigned>(threads, static_cast<unsigned>(workload_units)));
}

i64 floor_div(const i64 a, const i64 b) {
    if (a >= 0LL) {
        return a / b;
    }
    return -(((-a) + b - 1LL) / b);
}

u64 count_interval(const i64 lo, const i64 hi) {
    if (lo > hi) {
        return 0ULL;
    }
    return static_cast<u64>(hi - lo + 1LL);
}

u64 count_angle_a_for_pair(const i64 u, const i64 v, const i64 x_bound) {
    const i64 lo = std::max(-x_bound, v - x_bound);
    const i64 hi = std::min({x_bound, u + x_bound, floor_div(u - 1LL, 2LL)});
    return count_interval(lo, hi);
}

u64 count_angle_b_for_pair(const i64 u, const i64 v, const i64 x_bound) {
    const i64 lo = std::max({-x_bound, v - x_bound, floor_div(u, 2LL) + 1LL});
    const i64 hi = std::min({x_bound, u + x_bound, floor_div(v - 1LL, 2LL)});
    return count_interval(lo, hi);
}

u64 count_angle_c_for_pair(const i64 u, const i64 v, const i64 x_bound) {
    const i64 lo = std::max({-x_bound, v - x_bound, floor_div(v, 2LL) + 1LL});
    const i64 hi = std::min(x_bound, u + x_bound);
    return count_interval(lo, hi);
}

bool valid_triangle_from_sums(const i64 s1, const i64 s2, const i64 s3, const i64 x_bound) {
    if (!(s1 < s2 && s2 < s3)) {
        return false;
    }

    const i64 parity_sum = s1 + s2 + s3;
    if ((parity_sum & 1LL) != 0LL) {
        return false;
    }

    const i64 a = (s1 + s2 - s3) / 2LL;
    const i64 b = (s1 + s3 - s2) / 2LL;
    const i64 c = (s2 + s3 - s1) / 2LL;

    if (a < -x_bound || a > x_bound) {
        return false;
    }
    if (b < -x_bound || b > x_bound) {
        return false;
    }
    if (c < -x_bound || c > x_bound) {
        return false;
    }

    return a < b && b < c;
}

std::string to_string_u128(u128 value) {
    if (value == 0U) {
        return "0";
    }

    std::string digits;
    while (value > 0U) {
        const unsigned digit = static_cast<unsigned>(value % 10U);
        digits.push_back(static_cast<char>('0' + digit));
        value /= 10U;
    }
    std::reverse(digits.begin(), digits.end());
    return digits;
}

std::vector<int> build_spf(const u64 limit) {
    std::vector<int> spf(static_cast<std::size_t>(limit) + 1ULL, 0);
    if (limit >= 1ULL) {
        spf[1] = 1;
    }

    for (u64 i = 2ULL; i <= limit; ++i) {
        if (spf[static_cast<std::size_t>(i)] == 0) {
            spf[static_cast<std::size_t>(i)] = static_cast<int>(i);
            if (i <= limit / i) {
                for (u64 j = i * i; j <= limit; j += i) {
                    if (spf[static_cast<std::size_t>(j)] == 0) {
                        spf[static_cast<std::size_t>(j)] = static_cast<int>(i);
                    }
                }
            }
        }
    }

    return spf;
}

void collect_divisors_recursive(const std::vector<std::pair<u64, int>>& factors,
                                const std::size_t index,
                                const u64 current,
                                std::vector<u64>& divisors) {
    if (index == factors.size()) {
        divisors.push_back(current);
        return;
    }

    const u64 prime = factors[index].first;
    const int exponent = factors[index].second;

    u64 value = current;
    for (int e = 0; e <= exponent; ++e) {
        collect_divisors_recursive(factors, index + 1ULL, value, divisors);
        if (e != exponent) {
            value *= prime;
        }
    }
}

void factor_2k_squared(const u64 k,
                       const std::vector<int>& spf,
                       std::vector<std::pair<u64, int>>& factors) {
    factors.clear();

    u64 n = k;
    int exponent_two_in_k = 0;
    while ((n & 1ULL) == 0ULL) {
        ++exponent_two_in_k;
        n >>= 1U;
    }

    factors.emplace_back(2ULL, 2 * exponent_two_in_k + 1);

    while (n > 1ULL) {
        const int p = spf[static_cast<std::size_t>(n)];
        int exp_in_k = 0;
        while (n % static_cast<u64>(p) == 0ULL) {
            n /= static_cast<u64>(p);
            ++exp_in_k;
        }
        factors.emplace_back(static_cast<u64>(p), 2 * exp_in_k);
    }
}

void build_pairs_for_k(const u64 k,
                       const i64 x_bound,
                       const std::vector<u64>& divisors,
                       std::vector<std::pair<i64, i64>>& plus_pairs,
                       std::vector<std::pair<i64, i64>>& minus_pairs,
                       std::vector<std::pair<i64, i64>>& plus_out,
                       std::vector<std::pair<i64, i64>>& plus_in,
                       std::vector<std::pair<i64, i64>>& minus_out,
                       std::vector<std::pair<i64, i64>>& minus_in) {
    plus_pairs.clear();
    minus_pairs.clear();
    plus_out.clear();
    plus_in.clear();
    minus_out.clear();
    minus_in.clear();

    const i64 kk = static_cast<i64>(k);
    const i64 n = static_cast<i64>(2ULL * k * k);
    const i64 sum_bound = 2LL * x_bound;

    plus_pairs.reserve(divisors.size() * 2ULL);
    minus_pairs.reserve(divisors.size() * 2ULL);

    for (const u64 d_u64 : divisors) {
        const i64 d = static_cast<i64>(d_u64);
        for (const int sign : {-1, 1}) {
            const i64 p = static_cast<i64>(sign) * d;
            const i64 q = -n / p;

            const i64 u_plus = p + kk;
            const i64 v_plus = q - kk;
            if (u_plus < v_plus &&
                -sum_bound <= u_plus && u_plus <= sum_bound &&
                -sum_bound <= v_plus && v_plus <= sum_bound) {
                plus_pairs.emplace_back(u_plus, v_plus);
            }

            const i64 u_minus = p - kk;
            const i64 v_minus = q + kk;
            if (u_minus < v_minus &&
                -sum_bound <= u_minus && u_minus <= sum_bound &&
                -sum_bound <= v_minus && v_minus <= sum_bound) {
                minus_pairs.emplace_back(u_minus, v_minus);
            }
        }
    }

    std::sort(plus_pairs.begin(), plus_pairs.end());
    plus_pairs.erase(std::unique(plus_pairs.begin(), plus_pairs.end()), plus_pairs.end());

    std::sort(minus_pairs.begin(), minus_pairs.end());
    minus_pairs.erase(std::unique(minus_pairs.begin(), minus_pairs.end()), minus_pairs.end());

    plus_out = plus_pairs;
    minus_out = minus_pairs;

    plus_in.reserve(plus_pairs.size());
    for (const auto& [u, v] : plus_pairs) {
        plus_in.emplace_back(v, u);
    }
    std::sort(plus_in.begin(), plus_in.end());

    minus_in.reserve(minus_pairs.size());
    for (const auto& [u, v] : minus_pairs) {
        minus_in.emplace_back(v, u);
    }
    std::sort(minus_in.begin(), minus_in.end());
}

u128 solve_range(const u64 thread_index,
                 const u64 thread_count,
                 const u64 k_limit,
                 const i64 x_bound,
                 const std::vector<int>& spf) {
    u128 subtotal = 0U;

    std::vector<std::pair<u64, int>> factors;
    std::vector<u64> divisors;

    std::vector<std::pair<i64, i64>> plus_pairs;
    std::vector<std::pair<i64, i64>> minus_pairs;
    std::vector<std::pair<i64, i64>> plus_out;
    std::vector<std::pair<i64, i64>> plus_in;
    std::vector<std::pair<i64, i64>> minus_out;
    std::vector<std::pair<i64, i64>> minus_in;

    factors.reserve(16ULL);
    divisors.reserve(2048ULL);

    for (u64 k = 1ULL + thread_index; k <= k_limit; k += thread_count) {
        factor_2k_squared(k, spf, factors);

        divisors.clear();
        collect_divisors_recursive(factors, 0ULL, 1ULL, divisors);

        build_pairs_for_k(k,
                          x_bound,
                          divisors,
                          plus_pairs,
                          minus_pairs,
                          plus_out,
                          plus_in,
                          minus_out,
                          minus_in);

        u128 count_a = 0U;
        u128 count_b = 0U;
        u128 count_c = 0U;

        for (const auto& [u, v] : plus_out) {
            count_a += static_cast<u128>(count_angle_a_for_pair(u, v, x_bound));
            count_c += static_cast<u128>(count_angle_c_for_pair(u, v, x_bound));
        }

        for (const auto& [u, v] : minus_out) {
            count_b += static_cast<u128>(count_angle_b_for_pair(u, v, x_bound));
        }

        u128 overlap_ab = 0U;
        {
            std::size_t i = 0ULL;
            std::size_t j = 0ULL;
            while (i < plus_out.size() && j < minus_out.size()) {
                const i64 key_plus = plus_out[i].first;
                const i64 key_minus = minus_out[j].first;

                if (key_plus < key_minus) {
                    ++i;
                } else if (key_minus < key_plus) {
                    ++j;
                } else {
                    const i64 x = key_plus;
                    const i64 y = plus_out[i].second;
                    const i64 z = minus_out[j].second;
                    if (valid_triangle_from_sums(x, y, z, x_bound)) {
                        ++overlap_ab;
                    }
                    ++i;
                    ++j;
                }
            }
        }

        u128 overlap_ac = 0U;
        {
            std::size_t i = 0ULL;
            std::size_t j = 0ULL;
            while (i < plus_in.size() && j < plus_out.size()) {
                const i64 key_in = plus_in[i].first;
                const i64 key_out = plus_out[j].first;

                if (key_in < key_out) {
                    ++i;
                } else if (key_out < key_in) {
                    ++j;
                } else {
                    const i64 y = key_in;
                    const i64 x = plus_in[i].second;
                    const i64 z = plus_out[j].second;
                    if (valid_triangle_from_sums(x, y, z, x_bound)) {
                        ++overlap_ac;
                    }
                    ++i;
                    ++j;
                }
            }
        }

        u128 overlap_bc = 0U;
        {
            std::size_t i = 0ULL;
            std::size_t j = 0ULL;
            while (i < plus_in.size() && j < minus_in.size()) {
                const i64 key_plus = plus_in[i].first;
                const i64 key_minus = minus_in[j].first;

                if (key_plus < key_minus) {
                    ++i;
                } else if (key_minus < key_plus) {
                    ++j;
                } else {
                    const i64 z = key_plus;
                    const i64 y = plus_in[i].second;
                    const i64 x = minus_in[j].second;
                    if (valid_triangle_from_sums(x, y, z, x_bound)) {
                        ++overlap_bc;
                    }
                    ++i;
                    ++j;
                }
            }
        }

        subtotal += count_a + count_b + count_c - overlap_ab - overlap_ac - overlap_bc;
    }

    return subtotal;
}

u128 solve_f(const u64 k_limit,
             const i64 x_bound,
             const bool allow_multithreading,
             const unsigned requested_threads,
             const std::vector<int>& spf) {
    const unsigned thread_count =
        choose_thread_count(allow_multithreading, requested_threads, static_cast<std::size_t>(k_limit));

    if (thread_count == 1U) {
        return solve_range(0ULL, 1ULL, k_limit, x_bound, spf);
    }

    std::vector<u128> partial(thread_count, 0U);
    std::vector<std::thread> workers;
    workers.reserve(thread_count);

    for (unsigned t = 0U; t < thread_count; ++t) {
        workers.emplace_back([&, t]() {
            partial[t] = solve_range(static_cast<u64>(t),
                                     static_cast<u64>(thread_count),
                                     k_limit,
                                     x_bound,
                                     spf);
        });
    }

    for (std::thread& worker : workers) {
        worker.join();
    }

    u128 total = 0U;
    for (const u128 value : partial) {
        total += value;
    }
    return total;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        std::cerr << "Usage: " << argv[0]
                  << " [--k=<value>] [--x=<value>] [--threads=<value>]"
                  << " [--single-thread] [--skip-checkpoints]\n";
        return 1;
    }

    u64 max_k_needed = options.k;
    if (options.run_checkpoints) {
        for (const Checkpoint& cp : kCheckpoints) {
            max_k_needed = std::max(max_k_needed, cp.k);
        }
    }

    const std::vector<int> spf = build_spf(max_k_needed);

    if (options.run_checkpoints) {
        for (const Checkpoint& cp : kCheckpoints) {
            const u128 actual =
                solve_f(cp.k, cp.x, options.allow_multithreading, options.requested_threads, spf);
            const u128 expected = static_cast<u128>(cp.expected);
            std::cout << "Checkpoint F(" << cp.k << ", " << cp.x << ") = "
                      << to_string_u128(actual) << " (expected " << cp.expected << ")\n";
            if (actual != expected) {
                std::cerr << "Checkpoint mismatch.\n";
                return 1;
            }
        }
    }

    const u128 answer =
        solve_f(options.k, options.x, options.allow_multithreading, options.requested_threads, spf);

    std::cout << "F(" << options.k << ", " << options.x << ") = " << to_string_u128(answer)
              << '\n';

    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;
import java.util.stream.LongStream;

public class Euler397 {
    static long floorDiv(long a, long b) {
        if (a >= 0)
            return a / b;
        return -(((-a) + b - 1) / b);
    }

    static long countInterval(long lo, long hi) {
        if (lo > hi)
            return 0;
        return hi - lo + 1;
    }

    static long countAngleA(long u, long v, long xBound) {
        long lo = Math.max(-xBound, v - xBound);
        long hi = Math.min(xBound, Math.min(u + xBound, floorDiv(u - 1, 2)));
        return countInterval(lo, hi);
    }

    static long countAngleB(long u, long v, long xBound) {
        long lo = Math.max(-xBound, Math.max(v - xBound, floorDiv(u, 2) + 1));
        long hi = Math.min(xBound, Math.min(u + xBound, floorDiv(v - 1, 2)));
        return countInterval(lo, hi);
    }

    static long countAngleC(long u, long v, long xBound) {
        long lo = Math.max(-xBound, Math.max(v - xBound, floorDiv(v, 2) + 1));
        long hi = Math.min(xBound, u + xBound);
        return countInterval(lo, hi);
    }

    static boolean validTriangle(long s1, long s2, long s3, long xBound) {
        if (!(s1 < s2 && s2 < s3))
            return false;
        long parity = s1 + s2 + s3;
        if ((parity & 1) != 0)
            return false;

        long a = (s1 + s2 - s3) / 2;
        long b = (s1 + s3 - s2) / 2;
        long c = (s2 + s3 - s1) / 2;

        if (a < -xBound || a > xBound)
            return false;
        if (b < -xBound || b > xBound)
            return false;
        if (c < -xBound || c > xBound)
            return false;

        return a < b && b < c;
    }

    static class Factor {
        long p;
        int e;

        Factor(long p, int e) {
            this.p = p;
            this.e = e;
        }
    }

    static class Pair implements Comparable<Pair> {
        long first, second;

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

        @Override
        public int compareTo(Pair o) {
            if (this.first != o.first)
                return Long.compare(this.first, o.first);
            return Long.compare(this.second, o.second);
        }

        @Override
        public boolean equals(Object o) {
            if (!(o instanceof Pair))
                return false;
            Pair p = (Pair) o;
            return this.first == p.first && this.second == p.second;
        }

        @Override
        public int hashCode() {
            return Long.hashCode(first) ^ Long.hashCode(second);
        }
    }

    static void getDivisors(List<Factor> factors, int idx, long current, List<Long> result) {
        if (idx == factors.size()) {
            result.add(current);
            return;
        }
        long p = factors.get(idx).p;
        int exp = factors.get(idx).e;
        long val = current;
        for (int i = 0; i <= exp; i++) {
            getDivisors(factors, idx + 1, val, result);
            if (i != exp)
                val *= p;
        }
    }

    static long solveForK(long k, long xBound, int[] spf) {
        long n_k = k;
        int exp2 = 0;
        while ((n_k & 1) == 0) {
            exp2++;
            n_k >>= 1;
        }

        List<Factor> factors = new ArrayList<>();
        factors.add(new Factor(2, 2 * exp2 + 1));

        while (n_k > 1) {
            int p = spf[(int) n_k];
            int expP = 0;
            while (n_k % p == 0) {
                n_k /= p;
                expP++;
            }
            factors.add(new Factor(p, 2 * expP));
        }

        List<Long> divisors = new ArrayList<>();
        getDivisors(factors, 0, 1L, divisors);

        long kk = k;
        long nVal = 2 * k * k;
        long sumBound = 2 * xBound;

        List<Pair> plusOut = new ArrayList<>();
        List<Pair> minusOut = new ArrayList<>();

        for (long d : divisors) {
            for (int sign : new int[] { -1, 1 }) {
                long p = sign * d;
                long q = -nVal / p;

                long uPlus = p + kk;
                long vPlus = q - kk;
                if (uPlus < vPlus && uPlus >= -sumBound && uPlus <= sumBound && vPlus >= -sumBound
                        && vPlus <= sumBound) {
                    plusOut.add(new Pair(uPlus, vPlus));
                }

                long uMinus = p - kk;
                long vMinus = q + kk;
                if (uMinus < vMinus && uMinus >= -sumBound && uMinus <= sumBound && vMinus >= -sumBound
                        && vMinus <= sumBound) {
                    minusOut.add(new Pair(uMinus, vMinus));
                }
            }
        }

        Collections.sort(plusOut);
        List<Pair> plusOutUnique = new ArrayList<>();
        if (!plusOut.isEmpty()) {
            plusOutUnique.add(plusOut.get(0));
            for (int i = 1; i < plusOut.size(); i++) {
                if (!plusOut.get(i).equals(plusOutUnique.get(plusOutUnique.size() - 1))) {
                    plusOutUnique.add(plusOut.get(i));
                }
            }
        }
        plusOut = plusOutUnique;

        Collections.sort(minusOut);
        List<Pair> minusOutUnique = new ArrayList<>();
        if (!minusOut.isEmpty()) {
            minusOutUnique.add(minusOut.get(0));
            for (int i = 1; i < minusOut.size(); i++) {
                if (!minusOut.get(i).equals(minusOutUnique.get(minusOutUnique.size() - 1))) {
                    minusOutUnique.add(minusOut.get(i));
                }
            }
        }
        minusOut = minusOutUnique;

        List<Pair> plusIn = new ArrayList<>(plusOut.size());
        for (Pair p : plusOut)
            plusIn.add(new Pair(p.second, p.first));
        Collections.sort(plusIn);

        List<Pair> minusIn = new ArrayList<>(minusOut.size());
        for (Pair p : minusOut)
            minusIn.add(new Pair(p.second, p.first));
        Collections.sort(minusIn);

        long countA = 0, countB = 0, countC = 0;
        for (Pair p : plusOut) {
            countA += countAngleA(p.first, p.second, xBound);
            countC += countAngleC(p.first, p.second, xBound);
        }
        for (Pair p : minusOut) {
            countB += countAngleB(p.first, p.second, xBound);
        }

        long overlapAB = 0;
        int i = 0, j = 0;
        while (i < plusOut.size() && j < minusOut.size()) {
            long keyP = plusOut.get(i).first;
            long keyM = minusOut.get(j).first;
            if (keyP < keyM)
                i++;
            else if (keyM < keyP)
                j++;
            else {
                long x = keyP;
                long y = plusOut.get(i).second;
                long z = minusOut.get(j).second;
                if (validTriangle(x, y, z, xBound))
                    overlapAB++;
                i++;
                j++;
            }
        }

        long overlapAC = 0;
        i = 0;
        j = 0;
        while (i < plusIn.size() && j < plusOut.size()) {
            long keyIn = plusIn.get(i).first;
            long keyOut = plusOut.get(j).first;
            if (keyIn < keyOut)
                i++;
            else if (keyOut < keyIn)
                j++;
            else {
                long y = keyIn;
                long x = plusIn.get(i).second;
                long z = plusOut.get(j).second;
                if (validTriangle(x, y, z, xBound))
                    overlapAC++;
                i++;
                j++;
            }
        }

        long overlapBC = 0;
        i = 0;
        j = 0;
        while (i < plusIn.size() && j < minusIn.size()) {
            long keyP = plusIn.get(i).first;
            long keyM = minusIn.get(j).first;
            if (keyP < keyM)
                i++;
            else if (keyM < keyP)
                j++;
            else {
                long z = keyP;
                long y = plusIn.get(i).second;
                long x = minusIn.get(j).second;
                if (validTriangle(x, y, z, xBound))
                    overlapBC++;
                i++;
                j++;
            }
        }

        return countA + countB + countC - overlapAB - overlapAC - overlapBC;
    }

    static String solve() {
        long kLimit = 1000000L;
        long xBound = 1000000000L;

        int[] spf = new int[(int) kLimit + 1];
        for (int i = 2; i <= kLimit; i++) {
            if (spf[i] == 0) {
                spf[i] = i;
                if ((long) i * i <= kLimit) {
                    for (int j = i * i; j <= kLimit; j += i) {
                        if (spf[j] == 0)
                            spf[j] = i;
                    }
                }
            }
        }

        long total = LongStream.rangeClosed(1, kLimit)
                .parallel()
                .map(k -> solveForK(k, xBound, spf))
                .sum();

        return Long.toString(total);
    }

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