Problem 390: Triangles with Non Rational Sides and Integral Area

Project Euler Problem 390 Solution

EulerSolve provides an optimized solution for Project Euler Problem 390, Triangles with Non Rational Sides and Integral Area, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The repository solutions work with the triangle family whose side lengths are $$\sqrt{b^2+1},\qquad \sqrt{c^2+1},\qquad \sqrt{b^2+c^2},$$ where $b$ and $c$ are positive integers and the brute-force validator counts only unordered pairs $b \le c$. If $A$ denotes the area, the goal is $$S(L)=\sum A,$$ summed over all triangles in this family whose area is an integer and satisfies $A \le L$. Mathematical Approach Step 1: Reduce the Geometry to an Integer Square Test Let $$x=\sqrt{b^2+1},\qquad y=\sqrt{c^2+1},\qquad z=\sqrt{b^2+c^2}.$$ Applying Heron's identity in the form $$16A^2=2x^2y^2+2y^2z^2+2z^2x^2-x^4-y^4-z^4$$ and substituting $x^2=b^2+1$, $y^2=c^2+1$, $z^2=b^2+c^2$ gives $$16A^2=4\left(b^2c^2+b^2+c^2\right).$$ Therefore $$A=\frac{1}{2}\sqrt{b^2c^2+b^2+c^2}.$$ This is exactly the relation used by the C++ brute-force checkpoint. Writing $$m^2=b^2c^2+b^2+c^2,$$ we have $A=m/2$. So the area is integral precisely when $m$ is an even integer square root of that expression. Step 2: Parity Forces Even Parameters Because $A$ is an integer, $m=2A$ is even, so $m^2 \equiv 0 \pmod 4$. But a square is only $0$ or $1$ modulo $4$. If either $b$ or $c$ were odd, then $$b^2c^2+b^2+c^2 \equiv 1 \text{ or } 3 \pmod 4,$$ which is impossible. Hence every valid solution has even $b$ and even $c$....

Detailed mathematical approach

Problem Summary

For positive integers $b$ and $c$ we are given the triangle with sides

$$\sqrt{1+b^2},\quad \sqrt{1+c^2},\quad \sqrt{b^2+c^2}.$$

The first two sides are always irrational, since $b^2\lt1+b^2\lt(b+1)^2$, yet the area can be a whole number: the sides $\sqrt5$, $\sqrt{65}$ and $\sqrt{68}$, which belong to $b=2$ and $c=8$, enclose the area $9$. Let $S(n)$ be the sum of the areas of all triangles of this family whose area is an integer not exceeding $n$. Exchanging $b$ and $c$ leaves the three sides unchanged, so every triangle contributes once. The statement gives $S(10^6)=18018206$ and asks for $S(10^{10})$. ^[1]

The whole solution in one sentence: the area equals $\tfrac12\sqrt{b^2c^2+b^2+c^2}$, so an integral area forces $b=2p$ and $c=2q$ with $A^2=4p^2q^2+p^2+q^2$, and three explicit transformations grow every solution of that equation out of a few trivial ones, which lets us visit the qualifying triangles one by one instead of searching for them.

Mathematical Approach

1. The triangle is the slanted face of a box corner

Put three points on the coordinate axes, at $(1,0,0)$, $(0,b,0)$ and $(0,0,c)$. The three distances between them are

$$\sqrt{1+b^2},\quad \sqrt{b^2+c^2},\quad \sqrt{1+c^2},$$

exactly the three given sides, so the triangle always exists and can be measured in space. Together with the origin the three points span a tetrahedron whose corner at the origin is right-angled, and the three faces meeting there have the areas $b/2$, $c/2$ and $bc/2$. De Gua's theorem, the three-dimensional Pythagoras, states that the squared area of the face opposite that corner is the sum of the squared areas of the other three, ^[2] so the area $A$ satisfies

$$A^2=\frac{b^2}{4}+\frac{c^2}{4}+\frac{b^2c^2}{4}.$$

Heron's formula confirms it. In its squared-sides form ^[3]

$$\begin{aligned}16A^2={}&2x^2y^2+2y^2z^2+2z^2x^2\\&-x^4-y^4-z^4,\end{aligned}$$

the substitution $x^2=1+b^2$, $y^2=1+c^2$, $z^2=b^2+c^2$ collapses the right-hand side to $4(b^2c^2+b^2+c^2)$. We keep the doubled area $m=2A$:

$$m^2=b^2c^2+b^2+c^2.$$

For $b=2$ and $c=8$ this gives $m^2=324$, $m=18$ and $A=9$, the value quoted in the statement.

2. An integral area forces both parameters to be even

The area is an integer exactly when $m^2$ is a perfect square and its root is even. Every square is $0$ or $1$ modulo $4$, and that alone fixes the parity of $b$ and $c$:

both odd: $m^2\equiv1+1+1\equiv3\pmod4$, so $m^2$ is not a square at all;
exactly one odd: the even parameter makes $b^2c^2$ and its own square vanish modulo $4$, leaving $m^2\equiv1\pmod4$, so a square root of $m^2$ is odd and the area a half-integer;
both even: $m^2\equiv0\pmod4$, the only surviving case.

The smallest witness of the middle case is $b=1$, $c=2$, with sides $\sqrt2$, $\sqrt5$, $\sqrt5$: here $m^2=9$ is a perfect square, but $m=3$ and the area is $3/2$. Such near misses are the rule: among the pairs $b\le c\le600$, $622$ make $m^2$ a square with an odd root and only $14$ give an integral area — in all $14$ both parameters are even, as the argument demands.

3. The equation in half parameters

Write $b=2p$ and $c=2q$ with integers $p,q\ge1$. Then $m^2=16p^2q^2+4p^2+4q^2$, and $m=2A$ turns this into the equation the programs actually work with:

$$A^2=4p^2q^2+p^2+q^2.$$

Three consequences follow.

$A\gt 2pq$, because $A^2$ exceeds $4p^2q^2$. Hence $bc=4pq\lt 2A\le 2n$ for every triangle that counts.
$p\ne q$. For $p=q$ the equation reads $A^2=2p^2(2p^2+1)$; the two factors are coprime, so both would have to be squares, and $2p^2$ never is. A direct test confirms this for all $p\le200000$.
Ordering the parameters bounds them, by $b\lt\sqrt{2n}$ and $c\lt2n/b$, but not usefully: a search over the pairs with $bc\le2n$ would still have to sift $238734365497$ of them for $n=10^{10}$, and no small cutoff on $c$ is admissible, because the triangle with $b=2$ and $c=4807526976$ has the integral area $5374978561$.

4. For a fixed parameter, a Pell equation

Fix $p$ and collect the terms containing $q$. With $D_p=4p^2+1$ the equation becomes

$$A^2-D_pq^2=p^2.$$

An equation $u^2-Dv^2=N$ with a fixed non-square $D\gt0$ is called a generalised Pell equation, and $D_p$ is never a square because $(2p)^2\lt D_p\lt(2p+1)^2$; for $p=1$ ours reads $A^2-5q^2=1$, solved by the $(q,A)=(4,9)$ of our first triangle. For every $p$ it has the trivial solution $q=0$, $A=p$, which is no triangle — it would need $c=0$ — but is the seed of everything that follows. Classically all solutions arise from finitely many representatives by multiplication with the units of $\mathbb Z[\sqrt{D_p}]$; ^[4] we shall not need that classification, since the descent of Section 9 proves directly that the forest below misses nothing.

5. The unit that drives the recurrence

The steps come from the units of $\mathbb Z[\sqrt{D_p}]$, the numbers $x+y\sqrt{D_p}$ with $x^2-D_py^2=\pm1$: a unit of norm $+1$ preserves the right-hand side of our equation, one of norm $-1$ flips its sign and has to be squared first. A unit with $y=0$ is $\pm1$; for $y=1$ we need $x^2=D_p\pm1$, and of those two values $4p^2+2$ is no square while $4p^2$ is, which gives

$$(2p)^2-D_p\cdot1^2=-1,$$

and a larger $y$ gives a larger unit, so $2p+\sqrt{D_p}$ is the smallest unit above $1$: the negative Pell equation is solvable for every $D_p=4p^2+1$. Squaring a solution of the negative equation solves the positive one, ^[4] and

$$(2p+\sqrt{D_p})^2=(8p^2+1)+4p\sqrt{D_p},$$

where indeed $(8p^2+1)^2-D_p(4p)^2=1$, checked for all $p\le50$. The continued fraction agrees: rationalising gives $\sqrt{D_p}-2p=1/(\sqrt{D_p}+2p)$, so the next complete quotient $\sqrt{D_p}+2p$ has integer part $4p$ and the same fractional part again, whence $\sqrt{D_p}=[2p;\overline{4p}]$ — verified for $p=1,2,3,5,7,12$. ^[5]

Abbreviate $\alpha=8p^2+1$, $\beta=4p$ and $\gamma=4pD_p$. Multiplying $A+q\sqrt{D_p}$ by the unit $\alpha+\beta\sqrt{D_p}$ yields the step

$$q'=\alpha q+\beta A,\qquad A'=\gamma q+\alpha A.$$

Three identities are immediate: $\alpha^2-D_p\beta^2=1$, $\gamma=\beta D_p$ and $\gamma^2-D_p\alpha^2=-D_p$. With them

$$A'^2-D_pq'^2=A^2-D_pq^2,$$

because the mixed term vanishes with $\gamma=\beta D_p$. The step therefore carries a solution into another solution with the same $p$: a triangle into a triangle.

6. The first triangle of a seed, and how many seeds there are

Apply the step to the trivial solution $q=0$, $A=p$:

$$\begin{aligned}q_1&=4p^2,\\A_1&=(8p^2+1)p=8p^3+p.\end{aligned}$$

Every $p$ thus produces the triangle $b=2p$, $c=8p^2$ with the integral area $8p^3+p$. For $p=1,2,3,4$ these are $(b,c)=(2,8)$, $(4,32)$, $(6,72)$, $(8,128)$ with the areas $9$, $66$, $219$ and $516$; the first is the triangle of the statement.

Section 8 shows that every further step increases the area, so the tree grown from seed $p$ has no area below $8p^3+p$ and a seed with $8p^3+p\gt n$ can be skipped. For $n=10^{10}$ the condition $8p^3+p\le n$ holds exactly for $p\le1077$, because $8\cdot1077^3+1077=9993949341$ while the next seed already exceeds $10^{10}$. That number is at the same time the largest integral area not exceeding $10^{10}$ anywhere in the family.

7. Two roles for each coordinate, two directions for each role

The equation is symmetric in $p$ and $q$, so one solution satisfies two Pell equations at once:

$$A^2-(4p^2+1)q^2=p^2,$$

$$A^2-(4q^2+1)p^2=q^2.$$

Either coordinate may therefore be the fixed parameter of a step. Only the squares occur, so the sign of the partner is free as well, and flipping it reverses the direction of travel: stepping $(p,-q,A)$ gives the area $\alpha A-\gamma q$ and the partner $-(\alpha q-\beta A)$, which is the inverse step applied to $(p,q,A)$ up to the sign of the partner, and that sign leaves the triangle unchanged. Hence the three successors the programs push for a newly found solution $(p,q',A')$:

$(p,q',A')$ — keep the old parameter $p$ and step on;
$(q',p,A')$ — make the new value $q'$ the parameter and step forward;
$(q',-p,A')$ — make $q'$ the parameter and step the other way.

The fourth conceivable move, stepping backwards on the $p$-branch, only leads back to where the solution came from.

8. Every step increases the area

Call $(p,q,A)$ a state when $p\ge1$, $A\gt0$ and $A^2=4p^2q^2+p^2+q^2$, and call the state tidy when in addition $q\ge0$ or $|q|\lt p$.

Lemma. The step carries a tidy state into a solution with $q'\gt p$ and $A'\gt A$.

If $q\ge0$, then $A\ge p$ because $A^2=D_pq^2+p^2\ge p^2$, hence $q'=\alpha q+\beta A\ge4p^2\ge4p\gt p$ and $A'=\gamma q+\alpha A\ge\alpha A\ge9A$. Now let $q=-s$ with $1\le s\lt p$. Since $\alpha-1=8p^2$ and $\gamma=4p(4p^2+1)$, the claim $A'\gt A$ reads $\alpha A-\gamma s\gt A$ and is equivalent to $2pA\gt(4p^2+1)s$; both sides are positive, so we may square and use $A^2=D_ps^2+p^2$, which turns it into

$$4p^4\gt s^2(4p^2+1).$$

This holds because $s\le p-1$ gives $s^2(4p^2+1)\le4p^4-8p^3+5p^2-2p+1$. In the same way $q'\gt p$ reads $\beta A-\alpha s\gt p$ and becomes, after squaring,

$$s^2+2p(8p^2+1)s+p^2\lt16p^4,$$

whose left-hand side is at most $16p^4-16p^3+4p^2-4p+1$ for $s\le p-1$. Both polynomial inequalities are strict for every $p\ge1$, which proves the lemma.

The seeds $(p,0,p)$ are tidy, and all three successors of a tidy state are tidy again: the first two have a non-negative partner, and the third has the partner $-p$ with $p\lt q'$ by the lemma. By induction every state the traversal creates is tidy, and three things follow: the parameter of a state is always the smaller of the two coordinates of the solution it produces; the area grows along every edge, by a factor of at least $\alpha\ge9$ for the first two kinds of successor; and the traversal terminates, because the area passes any bound. At $n=10^{10}$ the smallest ratio between a child's area and its parent's is exactly $9$, at the seeds, and the next smallest is $121/9\approx13.44$.

9. Every triangle exactly once

Let $(p,q,A)$ be any solution with $1\le p\lt q$. The inverse step on the parameter $p$,

$$q^-=\alpha q-\beta A,\qquad A^-=\alpha A-\gamma q,$$

is a solution again, since the identities of Section 5 work both ways, and it is smaller:

$A^-\gt0$: from $\alpha^2-16p^2D_p=1$ we get $\alpha^2A^2-\gamma^2q^2=D_pq^2+\alpha^2p^2$, which is positive;
$A^-\lt A$: as in Section 8 this is equivalent to $4p^4\lt q^2(4p^2+1)$, true because $q\gt p$;
$|q^-|\lt q$, because $(q^-)^2=((A^-)^2-p^2)/D_p$ while $q^2=(A^2-p^2)/D_p$.

If $q^-=0$ we have reached the seed $(p,0,p)$, because then $A^-=p$. Otherwise replace $q^-$ by $|q^-|$ and put the smaller coordinate first; the result is a solution with two positive and, by Section 3, unequal coordinates, so the same rule applies to it. Call it the parent. For $(4,15,121)$, for instance, the inverse step gives $(4,-1,9)$, normalised to $(1,4,9)$. Every parent has a strictly smaller positive integral area, so a chain of parents cannot go on for ever and must end in a seed. ^[6]

Conversely the three successors of Section 7 are exactly the solutions whose parent is the given one. A solution $Z$ with parent $Y$ is the forward step of a state whose parameter is the smaller coordinate of $Z$ and whose solution is $Y$, so that parameter is one of $Y$'s two coordinates and the partner is the other one up to sign: four candidates. With $Y$'s smaller coordinate as the parameter, the two signs give the two neighbours of $Y$ on that branch, one with a larger and one with a smaller area, and only the larger one can be $Z$ — the first successor. With $Y$'s larger coordinate both signs are admissible and give the second and the third. The three are distinct: the first keeps $Y$'s smaller coordinate smallest, the others promote $Y$'s larger coordinate, and their two areas differ by $2\gamma p\gt0$.

The forest rooted at the seeds with $8p^3+p\le n$ therefore holds every triangle with integral area at most $n$ exactly once, so adding the area of every solution met is legitimate; the argument has the shape that generates all Markov triples from $(1,1,1)$. ^[7]

10. A worked example: $S(1000)=1092$

Here $8p^3+p\le1000$ holds for $p\le4$, so there are four seeds, with first areas $9$, $66$, $219$ and $516$. For $p=1$ we have $\alpha=9$, $\beta=4$, $\gamma=20$, and the walk starts at the state $(1,0,1)$:

$(1,0,1)$ gives $q'=4$, $A'=9$: the triangle $(b,c)=(2,8)$ of the statement;
$(1,4,9)$ gives $q'=9\cdot4+4\cdot9=72$ and $A'=20\cdot4+9\cdot9=161$: the triangle $(2,144)$;
$(4,1,9)$, now with $\alpha=129$, $\beta=16$, $\gamma=1040$, gives $q'=129+144=273$ and $A'=1040+1161=2201$, an area beyond $1000$, so this branch is cut;
$(4,-1,9)$ gives $q'=-129+144=15$ and $A'=-1040+1161=121$: the triangle $(8,30)$, confirmed by $8^2\cdot30^2+8^2+30^2=58564=242^2$ and $242/2=121$.

The seeds $p=2,3,4$ contribute only their first triangles $(4,32)$, $(6,72)$ and $(8,128)$ below $1000$. The six areas $9$, $66$, $121$, $161$, $219$, $516$ add up to $1092$, and a direct search over all pairs $b\le c$ with $bc\le2000$ finds exactly these six triangles. The negative partner is indispensable: the third smallest triangle of the family, $(8,30)$ with area $121$, is reached only through it.

How the Code Works

The C++ program ^[8] keeps the limit $n$, by default $10^{10}$, and the self-check flag in Options; parse_arguments reads --limit= and --skip-checkpoints and refuses a limit beyond the signed 64-bit range, which bounds the stored fields of a state but not the trial arithmetic: the widest intermediate needs $99$ bits at the published limit, while a limit near $2^{63}$ can push a coefficient out of the 128-bit range. seed_first_area(p) computes the $8p^3+p$ of Section 6 in 128 bits, and max_seed doubles a bound past the limit, then binary-searches the largest admissible $p$: $1077$ here.

A State holds the parameter, the partner and the area, exactly the data of Section 8. next_state computes the three coefficients of Section 5 — its local names a, b, c stand for $\alpha$, $\beta$, $\gamma$ and have nothing to do with the statement's side parameters — forms $q'$ and $A'$ as 128-bit products, rejects the state when $A'\le0$ or $A'\gt n$, and otherwise fills in the three successors of Section 7. sum_for_seed walks depth-first over a reserved std::vector<State> stack, adding $A'$ for every accepted state and pushing its three children, so each triangle is counted once by Section 9. solve determines $p_{\max}$, collects the contribution of every seed and adds the partial sums; to_string_u128 renders the total in decimal, digits collected in reverse and then reversed, for main to print.

The same file carries the validator of Sections 1 and 2: brute_sum_unordered runs over the pairs $b\le c\le\lfloor2n/b\rfloor+1$, slightly more than $bc\le2n$, forms $b^2c^2+b^2+c^2$, takes its integer square root with the binary search isqrt_u128, and adds the area only when the root is exact and even and the area itself is at most $n$. run_checkpoints verifies the statement's triangle ($m=18$, area $9$), compares the validator with the fast walk at $n=1000$, and reproduces $S(10^6)=18018206$.

Two guards in next_state never fire: the lemma of Section 8 gives $q'\gt0$, and $q'^2=(A'^2-p^2)/D_p$ keeps $q'$ below $A'/2$, inside the signed 64-bit range whenever $A'\le n$. The Python file is smaller: max_seed, sum_for_seed and a solve with the limit fixed at $10^{10}$, the recurrence written out inside the walk, unbounded integers instead of 128-bit types ^[9] and no parser, validator, checkpoints or guard on $q'$, which by the lemma changes nothing. Java also fixes the limit and packages no checks; it keeps State objects in an ArrayList, computes each trial step in BigInteger, narrows with longValue() after the limit test, and accumulates in long, ample for $1218$ areas below $10^{10}$. Its maxSeed doubles while hi stays at most $2000000$, a cap that would come too late in general, since $8p^3+p$ overflows a long already at $p=1048576$; at the published limit the doubling stops at $p=2048$ and nothing overflows.

Complexity and Verification

Let $T(n)$ be the number of triangles. The walk evaluates one state per seed and three per accepted solution, $p_{\max}+3T(n)$ in all, each with a bounded number of integer operations, so with $p_{\max}\sim(n/8)^{1/3}$ the operation count is $\Theta(n^{1/3}+T(n))$. At $n=10^{10}$ that is $1077+3\cdot1218=4731$ evaluations, an accepted depth of at most $8$ and at most $17$ pending states on the stack — the bound $2d+1$ for depth $d$, since every step leaves two deferred siblings. The operands are not bounded by $n^2$: a pending state only satisfies $p\lt n/2$, so its $\gamma=4p(4p^2+1)$ can reach the order of $n^3$, and at $n=1000$ the pending state $(72,1,161)$ already has $\gamma=5972256$. That is why the C++ computes in 128 bits and why one operation costs $O(\log n)$ bit work. The measured $T$ for $n=10,10^2,\dots,10^{10}$ is $1,2,6,17,35,74,145,293,589,1218$, close to a doubling per decade; that is an observation, not an asymptotic law, and the only bound proved here is $T(n)\ge p_{\max}$, since every admissible seed contributes its first triangle. On one desktop machine each program needs a few hundredths of a second, the C++ including its checkpoints, whose $\Theta(n\log n)$ validator is why they run at $n=1000$ and not higher.

Beyond the packaged checkpoints, the following was verified for this article.

A direct search over the pairs $b\le c$ with $bc\le2n$ reproduces $S$ at $n=10,10^2,\dots,10^6$: $9$, $75$, $1092$, $39696$, $794446$ and $18018206$, with $1,2,6,17,35,74$ triangles. The published Python agrees, and the C++ adds $S(10^7)=334160265$ and $S(10^8)=6808392896$.
An independent enumeration that shares no code with the forest: for every $p$ with $2p^2\lt n$ it tests $0\le s\lt p$ for an integral $A$ with $A^2-D_ps^2=p^2$, then multiplies both representatives $(s,A)$ and $(-s,A)$ by the unit while the area stays at most $n$, keeping the results with $q\gt p$. Both signs matter, just as in Section 7. It yields the same triangles as the walk at $n=10^6,10^7,10^8,10^9$ and, at $n=10^{10}$, the same $1218$ triples $(p,q,A)$, compared through a hash of the sorted list.
The structural claims at $n=10^{10}$: the invariant $A^2=4p^2q^2+p^2+q^2$ holds at every state, every emitted pair is distinct, the parameter is always the smaller coordinate, every edge raises the area, and the parent chain of each of the $1218$ solutions descends, in at most $8$ links, to one of the $1077$ seeds.
The area formula of Section 1 against the cross product of the edge vectors for random $b$ and $c$, and the parity count of Section 2.
The three published programs print the same value at the full limit, compared through their digests so that the value is not reproduced here.

Footnotes and References

The references supply the background facts used above. The lemma of Section 8, the parent map of Section 9 and the proof that the forest meets every triangle exactly once are derived in this article.

Project Euler 390 — Triangles with Non Rational Sides and Integral Area. The official statement gives the triangle with sides √5, √65 and √68, of area 9 and with b = 2, c = 8, defines S(n) as the sum of the areas of all triangles with sides √(1+b²), √(1+c²) and √(b²+c²) whose area is an integer not exceeding n, states S(10⁶) = 18018206 and asks for S(10¹⁰).
Wikipedia — De Gua's theorem. For a tetrahedron with a right-angled corner, the squared area of the opposite face equals the sum of the squared areas of the three faces meeting at that corner; the article presents it as the three-dimensional Pythagorean theorem. Section 1 applies it at the origin.
Wikipedia — Heron's formula. The article gives the area as a quarter of the root of 2(a²b² + a²c² + b²c²) − (a⁴ + b⁴ + c⁴), the form used in Section 1 to confirm the area independently of De Gua's theorem.
Wikipedia — Pell's equation. All solutions of x² − ny² = 1 are powers of the fundamental solution; the generalised equation x² − ny² = N has a finite set of solutions of which every solution is a Pell multiple; and squaring a solution of x² − ny² = −1 solves the positive equation. Sections 4 and 5 use all three.
Wikipedia — Periodic continued fraction. The continued fraction of the root of a non-square integer is eventually periodic of the shape [a₀; a₁, …, a₁, 2a₀] with a palindromic block, so a period of length one forces the single partial denominator 2a₀. For D = 4p² + 1 the expansion [2p; 4p, 4p, …] is derived in Section 5.
Wikipedia — Proof by infinite descent. The method rests on the impossibility of an infinite strictly decreasing sequence of positive integers: "there cannot be an infinity of ever-smaller natural numbers". Section 9 descends on the area.
Wikipedia — Markov number. A new Markov triple is obtained from (x, y, z) by replacing z with 3xy − z, and the graph joining each normalised triple to its neighbours is connected and starts at (1, 1, 1). Section 9 builds the analogue for our equation.
C++ — Euler390.cpp. The linked immutable revision contains max_seed, next_state, sum_for_seed, the validator brute_sum_unordered and the checkpoints. It is the implementation record for this article.
Python — numeric types documents that integers have unlimited precision, so the Python version needs no wide-integer type. Java — BigInteger documents the arbitrary-precision integers of the Java steps and warns that longValue() keeps only the low-order 64 bits of a value too large for a long; the bound on q′ above makes that narrowing safe.

Problem 390 source code

C++

#include <algorithm>
#include <atomic>
#include <cstdint>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>
#include <thread>
#include <vector>

namespace {

using i64 = std::int64_t;
using u64 = std::uint64_t;
using i128 = __int128_t;
using u128 = __uint128_t;

struct Options {
    u64 limit = 10000000000ULL;
    unsigned int threads = std::max(1U, std::thread::hardware_concurrency());
    bool run_checkpoints = true;
};

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

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

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

    value = parsed;
    return true;
}

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

    value = static_cast<unsigned int>(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 (parse_u64_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }
        if (parse_uint_after_prefix(arg, "--threads=", options.threads)) {
            continue;
        }

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

    if (options.threads == 0U) {
        options.threads = 1U;
    }

    if (options.limit > static_cast<u64>(std::numeric_limits<i64>::max())) {
        std::cerr << "--limit must be <= " << std::numeric_limits<i64>::max()
                  << " to keep intermediate signed states safe.\n";
        return false;
    }

    return true;
}

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

    std::string digits;
    while (value > 0U) {
        const int digit = static_cast<int>(value % 10U);
        digits.push_back(static_cast<char>('0' + digit));
        value /= 10U;
    }

    std::reverse(digits.begin(), digits.end());
    return digits;
}

u64 isqrt_u128(const u128 n) {
    u128 lo = 0U;
    u128 hi = static_cast<u128>(std::numeric_limits<u64>::max());

    while (lo < hi) {
        const u128 mid = (lo + hi + 1U) >> 1;
        if (mid <= n / mid) {
            lo = mid;
        } else {
            hi = mid - 1U;
        }
    }

    return static_cast<u64>(lo);
}

u128 seed_first_area(const i64 p) {
    const u128 pp = static_cast<u128>(p);
    return 8U * pp * pp * pp + pp;
}

i64 max_seed(const u64 limit) {
    i64 lo = 0;
    i64 hi = 1;

    while (seed_first_area(hi) <= static_cast<u128>(limit)) {
        if (hi > std::numeric_limits<i64>::max() / 2) {
            break;
        }
        hi *= 2;
    }

    while (lo < hi) {
        const i64 mid = lo + (hi - lo + 1) / 2;
        if (seed_first_area(mid) <= static_cast<u128>(limit)) {
            lo = mid;
        } else {
            hi = mid - 1;
        }
    }

    return lo;
}

struct State {
    i64 p;
    i64 q;
    u64 area;
};

bool next_state(const State& current,
                const u64 limit,
                State& same_p,
                State& swap_branch,
                State& sign_branch) {
    const i128 p = static_cast<i128>(current.p);
    const i128 p2 = p * p;

    const i128 a = 8 * p2 + 1;
    const i128 b = 4 * p;
    const i128 c = 4 * p * (4 * p2 + 1);

    const i128 q1 = a * static_cast<i128>(current.q) + b * static_cast<i128>(current.area);
    const i128 s1 = c * static_cast<i128>(current.q) + a * static_cast<i128>(current.area);

    if (s1 <= 0 || static_cast<u128>(s1) > static_cast<u128>(limit)) {
        return false;
    }
    if (q1 <= 0 || q1 > static_cast<i128>(std::numeric_limits<i64>::max())) {
        return false;
    }

    const i64 nq = static_cast<i64>(q1);
    const u64 ns = static_cast<u64>(s1);

    same_p = State{current.p, nq, ns};
    swap_branch = State{nq, current.p, ns};
    sign_branch = State{nq, -current.p, ns};
    return true;
}

u128 sum_for_seed(const i64 seed, const u64 limit) {
    u128 sum = 0U;
    std::vector<State> stack;
    stack.reserve(256);
    stack.push_back(State{seed, 0, static_cast<u64>(seed)});

    while (!stack.empty()) {
        const State current = stack.back();
        stack.pop_back();

        State same_p{};
        State swap_branch{};
        State sign_branch{};
        if (!next_state(current, limit, same_p, swap_branch, sign_branch)) {
            continue;
        }

        sum += static_cast<u128>(same_p.area);

        stack.push_back(same_p);
        stack.push_back(swap_branch);
        stack.push_back(sign_branch);
    }

    return sum;
}

u128 solve(const u64 limit, const unsigned int requested_threads) {
    const i64 pmax = max_seed(limit);
    if (pmax <= 0) {
        return 0U;
    }

    const unsigned int workers = std::max(
        1U,
        std::min<unsigned int>(requested_threads, static_cast<unsigned int>(pmax)));

    std::atomic<i64> next_seed(1);
    std::vector<u128> partials(workers, 0U);
    std::vector<std::thread> threads;
    threads.reserve(workers);

    for (unsigned int t = 0; t < workers; ++t) {
        threads.emplace_back([&, t]() {
            u128 local = 0U;
            while (true) {
                const i64 seed = next_seed.fetch_add(1, std::memory_order_relaxed);
                if (seed > pmax) {
                    break;
                }
                local += sum_for_seed(seed, limit);
            }
            partials[t] = local;
        });
    }

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

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

u128 brute_sum_unordered(const u64 limit) {
    if (limit == 0U) {
        return 0U;
    }

    u128 total = 0U;
    const u64 max_b = 2U * limit;

    for (u64 b = 1U; b <= max_b; ++b) {
        const u64 max_c = max_b / b + 1U;
        if (max_c < b) {
            continue;
        }

        for (u64 c = b; c <= max_c; ++c) {
            const u128 m2 =
                static_cast<u128>(b) * b * c * c + static_cast<u128>(b) * b + static_cast<u128>(c) * c;
            const u64 m = isqrt_u128(m2);
            if (static_cast<u128>(m) * m != m2) {
                continue;
            }
            if ((m & 1U) != 0U) {
                continue;
            }

            const u64 area = m / 2U;
            if (area <= limit) {
                total += area;
            }
        }
    }

    return total;
}

void run_checkpoints(const Options& options) {
    {
        constexpr u64 b = 2U;
        constexpr u64 c = 8U;
        const u128 m2 = static_cast<u128>(b) * b * c * c +
                        static_cast<u128>(b) * b +
                        static_cast<u128>(c) * c;
        const u64 m = isqrt_u128(m2);
        if (static_cast<u128>(m) * m != m2 || m != 18U) {
            throw std::runtime_error("Checkpoint failed: sample triangle should have m=18.");
        }
        if (m / 2U != 9U) {
            throw std::runtime_error("Checkpoint failed: sample triangle area should be 9.");
        }
    }

    {
        constexpr u64 small_limit = 1000U;
        const u128 brute = brute_sum_unordered(small_limit);
        const u128 fast = solve(small_limit, 1U);
        if (brute != fast) {
            throw std::runtime_error(
                "Checkpoint failed: recurrence and brute-force mismatch at N=1000.");
        }
    }

    {
        constexpr u64 statement_limit = 1000000U;
        constexpr u64 statement_value = 18018206U;
        const u128 sample = solve(statement_limit, options.threads);
        if (sample != static_cast<u128>(statement_value)) {
            throw std::runtime_error("Checkpoint failed: S(10^6) mismatch.");
        }
    }
}

}  // namespace

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

    try {
        if (options.run_checkpoints) {
            run_checkpoints(options);
        }

        const u128 answer = solve(options.limit, options.threads);
        std::cout << "S(" << options.limit << ") = " << to_string_u128(answer)
                  << '\n';
    } catch (const std::exception& ex) {
        std::cerr << ex.what() << '\n';
        return 1;
    }

    return 0;
}

Python

def max_seed(limit):
    lo = 0
    hi = 1
    def seed_first_area(p):
        return 8 * p * p * p + p
        
    while seed_first_area(hi) <= limit:
        hi *= 2
        
    while lo < hi:
        mid = (lo + hi + 1) // 2
        if seed_first_area(mid) <= limit:
            lo = mid
        else:
            hi = mid - 1
    return lo

def sum_for_seed(seed, limit):
    total = 0
    stack = [(seed, 0, seed)]
    
    while stack:
        p, q, area = stack.pop()
        
        p2 = p * p
        a = 8 * p2 + 1
        b = 4 * p
        c = 4 * p * (4 * p2 + 1)
        
        q1 = a * q + b * area
        s1 = c * q + a * area
        
        if s1 > 0 and s1 <= limit:
            total += s1
            stack.append((p, q1, s1))
            stack.append((q1, p, s1))
            stack.append((q1, -p, s1))
            
    return total

def solve():
    limit = 10000000000
    pmax = max_seed(limit)
    if pmax <= 0: return "0"
    
    ans = 0
    for seed in range(1, pmax + 1):
        ans += sum_for_seed(seed, limit)
        
    return str(ans)

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

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;

public class Euler390 {
    static long seedFirstArea(long p) {
        return 8 * p * p * p + p;
    }

    static long maxSeed(long limit) {
        long lo = 0;
        long hi = 1;
        while (true) {
            if (hi > 2000000)
                break;
            if (seedFirstArea(hi) <= limit)
                hi *= 2;
            else
                break;
        }

        while (lo < hi) {
            long mid = lo + (hi - lo + 1) / 2;
            if (seedFirstArea(mid) <= limit) {
                lo = mid;
            } else {
                hi = mid - 1;
            }
        }
        return lo;
    }

    static class State {
        long p;
        long q;
        long area;

        State(long p, long q, long area) {
            this.p = p;
            this.q = q;
            this.area = area;
        }
    }

    static long sumForSeed(long seed, long limit) {
        long total = 0;
        List<State> stack = new ArrayList<>();
        stack.add(new State(seed, 0, seed));

        BigInteger limitBi = BigInteger.valueOf(limit);

        while (!stack.isEmpty()) {
            State current = stack.remove(stack.size() - 1);
            BigInteger p = BigInteger.valueOf(current.p);
            BigInteger q = BigInteger.valueOf(current.q);
            BigInteger area = BigInteger.valueOf(current.area);

            BigInteger p2 = p.multiply(p);
            BigInteger a = BigInteger.valueOf(8).multiply(p2).add(BigInteger.ONE);
            BigInteger b = BigInteger.valueOf(4).multiply(p);
            BigInteger c = BigInteger.valueOf(4).multiply(p).multiply(
                    BigInteger.valueOf(4).multiply(p2).add(BigInteger.ONE));

            BigInteger q1 = a.multiply(q).add(b.multiply(area));
            BigInteger s1 = c.multiply(q).add(a.multiply(area));

            if (s1.compareTo(BigInteger.ZERO) > 0 && s1.compareTo(limitBi) <= 0) {
                long s1Long = s1.longValue();
                total += s1Long;

                long q1Long = q1.longValue();
                long pLong = current.p;

                stack.add(new State(pLong, q1Long, s1Long));
                stack.add(new State(q1Long, pLong, s1Long));
                stack.add(new State(q1Long, -pLong, s1Long));
            }
        }
        return total;
    }

    static String solve() {
        long limit = 10000000000L;
        long pmax = maxSeed(limit);
        if (pmax <= 0)
            return "0";

        long ans = 0;
        for (long seed = 1; seed <= pmax; seed++) {
            ans += sumForSeed(seed, limit);
        }
        return Long.toString(ans);
    }

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