Problem 251: Cardano Triplets

Project Euler Problem 251 Solution

EulerSolve provides an optimized solution for Project Euler Problem 251, Cardano Triplets, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A Cardano triplet is a triple of positive integers $(a,b,c)$ satisfying $$\sqrt[3]{a+b\sqrt{c}}+\sqrt[3]{a-b\sqrt{c}}=1.$$ The task is to count all such triplets with $a+b+c\le L$. The implementation does not search directly in $(a,b,c)$-space; it first converts the radical identity into a Diophantine equation and then enumerates the resulting factorization patterns. Mathematical Approach Let $$x=\sqrt[3]{a+b\sqrt{c}},\qquad y=\sqrt[3]{a-b\sqrt{c}}.$$ Then $x+y=1$. The whole solution comes from extracting arithmetic consequences of this identity. From the Cube-Root Identity to an Integer Equation Using $$x^3+y^3=(x+y)^3-3xy(x+y),$$ we get $$2a=x^3+y^3=1-3xy,$$ so $$xy=\frac{1-2a}{3}.$$ Multiplying the two cube-root arguments also gives $$x^3y^3=(a+b\sqrt{c})(a-b\sqrt{c})=a^2-b^2c,$$ hence $$a^2-b^2c=\left(\frac{1-2a}{3}\right)^3.$$ After rearranging, $$27b^2c=27a^2+(2a-1)^3=(a+1)^2(8a-1).$$ This factorization is the key algebraic simplification behind the code. Why $a=3u-1$ The left-hand side is divisible by $27$, so $(a+1)^2(8a-1)$ must also be divisible by $27$....

Detailed mathematical approach

Problem Summary

A triple of positive integers $(a,b,c)$ is called a Cardano triplet if

$$\sqrt[3]{a+b\sqrt c}+\sqrt[3]{a-b\sqrt c}=1.$$

The statement gives the example $(2,1,5)$ and the checkpoint that exactly $149$ Cardano triplets satisfy $a+b+c\le1000$. We must count the Cardano triplets with $a+b+c\le110\,000\,000$. ^[1] In the example the second number under a cube root, $2-\sqrt5$, is negative, so the cube roots are the real ones, as in $\sqrt[3]{-8}=-2$. Throughout, $L$ denotes the bound on $a+b+c$.

The whole solution in one sentence: the condition holds exactly when $a=3u-1$ and $b^2c=u^2(8u-3)$ for a positive integer $u$, so we run through every $u$ that the budget allows, factor $u$ and $8u-3$, and count the ways to write $u^2(8u-3)$ as a square $b^2$ times a number $c$ with $b+c\le L-a$.

Mathematical Approach

We turn the cube roots into an equation in integers and show that nothing is lost on the way, then count the solutions for one value of $a$ at a time, and finally find how far $a$ can go.

1. A first example: the golden ratio

Let $\varphi=(1+\sqrt5)/2$ and $\psi=(1-\sqrt5)/2$. Squaring gives $\varphi^2=\varphi+1$, hence $\varphi^3=\varphi^2+\varphi$${}=2\varphi+1=2+\sqrt5$, and in the same way $\psi^3=2\psi+1=2-\sqrt5$. Therefore

$$\sqrt[3]{2+\sqrt5}+\sqrt[3]{2-\sqrt5}=\varphi+\psi=1,$$

which is the example $(2,1,5)$. Here $2-\sqrt5\approx-0.236$, and its cube root $\psi\approx-0.618$ is the real cube root: every real number has exactly one real cube root, and the cube root of a negative number is negative. ^[2] A numerical test must therefore use a real cube root such as C's cbrt, because pow(x, 1.0/3) cannot raise a negative base to a fractional power and returns NaN. ^[3]

2. From cube roots to an equation in integers

Let $x$ and $y$ be the two real cube roots in the condition, so that $x^3+y^3=2a$. The product $xy$ is a real number whose cube is $(a+b\sqrt c)(a-b\sqrt c)=a^2-b^2c$, so $xy=\sqrt[3]{a^2-b^2c}$ by the uniqueness of real cube roots. Now suppose that $x+y=1$. The identity $(x+y)^3=x^3+y^3+3xy(x+y)$ turns into $1=2a+3xy$, so $xy=(1-2a)/3$. Cubing this and multiplying by $27$ gives $27a^2-27b^2c=(1-2a)^3$, that is,

$$27b^2c=27a^2+(2a-1)^3=(a+1)^2(8a-1),$$

where both right-hand sides expand to $8a^3+15a^2+6a-1$. For $(2,1,5)$ both sides equal $27\cdot5=9\cdot15=135$.

3. The shape of $a$

The left-hand side is a multiple of $27$, so $3$ divides $(a+1)^2(8a-1)$. The two factors are tied together modulo $3$: $(8a-1)-2(a+1)=3(2a-1)$, so $3$ divides $8a-1$ exactly when it divides $a+1$. Hence $3\mid a+1$, and we write $a=3u-1$ with a positive integer $u$. Then $a+1=3u$ and $8a-1=3(8u-3)$, the right-hand side is $27u^2(8u-3)$, and the equation becomes

$$b^2c=u^2m,\qquad m=8u-3.$$

So every Cardano triplet has $a\in\{2,5,8,11,\dots\}$. The example has $u=1$, $m=5$ and $b^2c=5=u^2m$.

4. Every solution of the equation is a triplet

The converse is true as well, and the cube roots can be written down explicitly. With $m=8u-3$, expanding the cube gives

$$\begin{aligned}\Bigl(\frac{1\pm\sqrt m}{2}\Bigr)^3&=\frac{1+3m}{8}\pm\frac{3+m}{8}\sqrt m\\&=(3u-1)\pm u\sqrt m.\end{aligned}$$

Now let $b$ and $c$ be positive integers with $b^2c=u^2m$. Then $b\sqrt c=u\sqrt m$, because both are positive and have the same square, so $a\pm b\sqrt c$ are exactly the two cubes above. Their real cube roots are $(1\pm\sqrt m)/2$, and these add up to $1$. Together with Sections 2 and 3 this proves the key fact.

Claim. $(a,b,c)$ is a Cardano triplet if and only if $a=3u-1$ and $b^2c=u^2(8u-3)$ for some positive integer $u$.

Since $m\ge5$, the second root $(1-\sqrt m)/2$ is always negative; equivalently, $a^2-b^2c=(xy)^3=(1-2u)^3\lt0$. The name comes from Cardano's formula, first published in his Ars Magna of 1545: if $t^3+pt+q=0$ with real $p$ and $q$ and $D=q^2/4+p^3/27\gt0$, the real root is ^[4]

$$t=\sqrt[3]{-\tfrac q2+\sqrt D}+\sqrt[3]{-\tfrac q2-\sqrt D}.$$

For $p=2a-1$ and $q=-2a$, Section 2 shows that $D=(a+1)^2(8a-1)/27$, which equals $b^2c$ for a triplet. So the left side of the condition is Cardano's formula for the real root of $t^3+(2a-1)t-2a$${}=(t-1)(t^2+t+2a)$, and since the quadratic factor has the negative discriminant $1-8a$, that root is $1$.

5. Small cases

For each $u$ the triplets come from the numbers $b$ with $b^2\mid u^2m$, and then $c=u^2m/b^2$. The first eight values of $u$:

$u=1$: $u^2m=5$, so $b=1$ and $(2,1,5)$ has sum $8$.
$u=2$: $52=2^2\cdot13$, $b\in\{1,2\}$: $(5,2,13)$ and $(5,1,52)$, with sums $20$ and $58$.
$u=3$: $189=3^3\cdot7$, $b\in\{1,3\}$: $(8,3,21)$ and $(8,1,189)$, with sums $32$ and $198$.
$u=4$: $464=2^4\cdot29$, $b\in\{1,2,4\}$: $(11,4,29)$, $(11,2,116)$ and $(11,1,464)$, with sums $44$, $129$ and $476$.
$u=5$: $925=5^2\cdot37$: $(14,5,37)$ and $(14,1,925)$, with sums $56$ and $940$.
$u=6$: $1620=2^2\cdot3^4\cdot5$, $b\in\{1,2,3,6,9,18\}$: $(17,18,5)$, $(17,9,20)$, $(17,6,45)$, $(17,3,180)$, $(17,2,405)$ and $(17,1,1620)$, with sums $40$, $46$, $68$, $200$, $424$ and $1638$.
$u=7$: $2597=7^2\cdot53$: sums $80$ and $2618$. $u=8$: $3904=2^6\cdot61$, $b\in\{1,2,4,8\}$: sums $92$, $271$, $1001$ and $3928$ for $b=8,4,2,1$.

The six triplets for $u=6$ are the six ways of writing $6\sqrt{45}$ as $b\sqrt c$: $6\sqrt{45}=18\sqrt5=9\sqrt{20}$${}=3\sqrt{180}=2\sqrt{405}=\sqrt{1620}$. The choice $b=u$, $c=m$ always works, so $(3u-1,u,8u-3)$ is a triplet with sum $12u-4$ for every $u$, and there are infinitely many Cardano triplets. For $u=9,\dots,14$ it is the cheapest triplet, with sums from $104$ to $164$, and Section 8 shows that no larger $u$ can reach a sum of at most $100$. Counting the sums up to $100$ in the list gives $1+2+1+1+1+3+1+1=11$ triplets for $L=100$.

6. Splitting $b$ into a part of $u$ and a part of $m$

The program never factors the large number $u^2m$. It factors $u$ and $m$ separately and assembles $b$ from two pieces. Let $g=\gcd(b,u)$ and write $u=gh$ and $b=gb_1$, so that $\gcd(h,b_1)=1$. Dividing $b^2c=u^2m$ by $g^2$ leaves $b_1^2c=h^2m$. Since $b_1^2$ divides $h^2m$ and is coprime to $h^2$, it divides $m$ by the generalised Euclid's lemma; ^[5] we write $m=b_1^2r$, and then $c=h^2r$. Conversely, for every divisor $h$ of $u$ and every $b_1$ with $b_1^2\mid m$ and $\gcd(h,b_1)=1$, put

$$b=\frac uh\,b_1,\qquad c=h^2\,\frac{m}{b_1^2}.$$

Then $b^2c=u^2m$, and $\gcd(b,u)=\frac uh\gcd(b_1,h)=\frac uh$, so the first construction returns exactly this $h$ and $b_1$. The triplets for one $u$ therefore correspond one to one to the pairs $(h,b_1)$ with $h\mid u$, $b_1^2\mid m$ and $\gcd(h,b_1)=1$.

For $u=6$ and $m=45=3^2\cdot5$ the candidates are $h\in\{1,2,3,6\}$ and $b_1\in\{1,3\}$. With $b_1=1$ the four values of $h$ give $b=6,3,2,1$ and $c=45,180,405,1620$. With $b_1=3$, so $r=5$, the values $h=1$ and $h=2$ give $(b,c)=(18,5)$ and $(9,20)$, while $h=3$ and $h=6$ would give $(6,45)$ and $(3,180)$ a second time. The condition $\gcd(h,b_1)=1$ removes exactly these two duplicates and leaves the six triplets of Section 5. It can only fail through the prime $3$: $\gcd(u,8u-3)=\gcd(u,3)$, so a prime that divides both $h\mid u$ and $b_1\mid m$ divides $3$. When $3\nmid u$, every pair is admissible.

7. The budget for one value of $u$

Put $R=L-a$, the part of the budget left for $b+c$. For a fixed $b_1$, with $r=m/b_1^2$, every pair has $b=\frac uh b_1\ge b_1$ and $c=h^2r\ge r$. This gives three cheap tests before the exact check $b+c\le R$:

if $b_1+r\gt R$, no $h$ can work;
$h^2r=c\le R-b\le R-b_1$, so $h\le\lfloor\sqrt{(R-b_1)/r}\rfloor$;
$ub_1/h=b\le R-c\le R-r$, so $h\ge\lceil ub_1/(R-r)\rceil$.

These bounds only discard pairs that cannot fit; the exact test decides the rest. For $L=100$ and $u=6$ we have $R=100-17=83$. For $b_1=1$ and $r=45$ the first test passes, since $46\le83$, but $h\le\lfloor\sqrt{82/45}\rfloor=1$, so only $(17,6,45)$ remains. For $b_1=3$ and $r=5$ the window is $1\le h\le\lfloor\sqrt{80/5}\rfloor=4$; it contains the divisors $1$, $2$ and $3$ of $6$, the divisor $3$ fails the gcd test, and $h=1,2$ give $(17,18,5)$ and $(17,9,20)$. These are the three triplets with $u=6$ and sum at most $100$.

8. How large $u$ can be

Write $b+c=\frac b2+\frac b2+c$. The inequality of arithmetic and geometric means for these three positive numbers gives ^[6]

$$b+c\ge3\sqrt[3]{\frac{b^2c}{4}}=3\sqrt[3]{\frac{u^2(8u-3)}{4}},$$

with equality exactly when $b/2=c$. Every triplet with parameter $u$ therefore has $a+b+c\ge B(u)$, where

$$B(u)=3u-1+3\sqrt[3]{u^2(8u-3)/4}$$

increases with $u$, because $3u-1$, $u^2$ and $8u-3$ all increase and the cube root preserves order. Only the values of $u$ up to the largest $u_{\max}$ with $B(u)\le L$ can contribute. In integers, $B(u)\le L$ means $27u^2(8u-3)\le4(L-3u+1)^3$ with $L-3u+1\ge0$. Since $B(u)/u$ tends to $3+3\sqrt[3]2\approx6.78$, $u_{\max}$ is close to $L/6.78$: it is $14$ for $L=100$, $147$ for $L=1000$ and $16\,224\,755$ for $L=110\,000\,000$. The bound is attained: for $u=16$ we have $u^2m=32000=4\cdot20^3$, and $(47,40,20)$ has $b=2c$ and sum exactly $B(16)=47+3\cdot20=107$. For small $L$ the bound is generous, though: the largest $u$ that really occurs is $8$ for $L=100$ and $141$ for $L=1000$, with the triplet $(422,423,125)$.

9. Why every triplet is counted exactly once

Sections 2 to 4 match the Cardano triplets one to one with the pairs $(u,b)$ such that $b^2\mid u^2(8u-3)$, Section 6 matches these with the triples $(u,h,b_1)$, and Sections 7 and 8 discard only candidates whose sum exceeds $L$. So counting the triples $(u,h,b_1)$ with $u\le u_{\max}$ that pass the gcd test and the exact test $b+c\le R$ counts every Cardano triplet with $a+b+c\le L$ exactly once.

How the Code Works

The C++ program ^[7] starts in main, where parse_arguments reads the optional arguments --limit= (default $110\,000\,000$) and --skip-checkpoints. Unless the checks are skipped, run_checkpoints calls solve for the limits $8$, $100$ and $1000$ and requires $1$, $11$ and $149$ triplets; on a mismatch it writes a message to standard error and the program exits with status $2$. Then main prints the count for the limit.

solve returns $0$ for a limit below $8$, the sum of the cheapest triplet $(2,1,5)$. max_u_bound finds $u_{\max}$ of Section 8 by binary search over $0\le u\le\lfloor(L+1)/3\rfloor$, which is valid because $B$ increases; can_have_solution evaluates $B(u)$ in long double with std::cbrt and accepts $u$ when $B(u)\le L+10^{-9}$. The small slack matters because the bound can be attained exactly, as for $(47,40,20)$. Two tables follow: sieve_primes lists the primes up to $\lfloor\sqrt{8u_{\max}-3}\rfloor+1$ with the sieve of Eratosthenes, and build_spf stores the smallest prime factor of every number up to $u_{\max}$ with the linear sieve. ^[8]

The main loop handles every $u$ from $1$ to $u_{\max}$ and adds up the counts. For each $u$ it sets $a=3u-1$ and remaining $=R=L-a$. factorize_with_spf splits $u$ into prime powers by repeated division by the stored smallest prime factor, and generate_divisors turns them into the list of all divisors $h$: for each prime $p$ with exponent $e$, the divisors found so far are multiplied by $p,p^2,\dots,p^e$ and appended. For $m=8u-3$, factorize_square_part uses trial division ^[9] by the listed primes while $p^2$ does not exceed the unfactored part $x$, and records $(p,\lfloor e/2\rfloor)$ for every prime with exponent $e\ge2$. The loop ends either when $p^2$ exceeds $x$ or when the list is used up, and since the list reaches $\sqrt m$, in both cases every prime up to $\sqrt x$ has been tried. So the rest $x$ is $1$ or a prime that occurs only once in $m$, and no square factor is missed. A second call of generate_divisors turns the recorded pairs into all $b_1$ with $b_1^2\mid m$.

For each $b_1$ the loop computes m_over $=r=m/b_1^2$, skips $b_1$ when $b_1+r\gt R$, and computes the window of Section 7: h_max with floor_sqrt_i64, a floating-point square root corrected by exact integer comparisons, and h_min as $\lceil ub_1/(R-r)\rceil$ by integer division. Every divisor $h$ of $u$ inside the window with std::gcd(h, b1) == 1 gives $b=(u/h)\,b_1$ and $c=h^2r$, and the triplet is counted when $b+c\le R$. The divisor list is not sorted, so every divisor is tested against the window. All quantities fit easily in 64-bit integers, since $ub_1$ stays below $2\cdot10^{11}$ and $c\le R$; the 128-bit types in floor_sqrt_i64 and in the final comparison are a safety margin.

The Java program repeats these steps with the same structure (sievePrimes, buildSPF, factorizeWithSPF, factorizeSquarePart, generateDivisors, maxUBound) and its own Euclidean gcd. It differs in three details: the limit is the constant LIMIT and no checkpoints are run; $B(u)$ is evaluated in double; and floorSqrt and the final test $b+c\le R$ compare BigInteger values, which is safe but not needed at these sizes. The Python file does not implement the method. Its solve takes the problem number from its own file name, compiles ../solutionsCpp/Euler251.cpp with clang++ or g++ (-std=c++17 -O2) into a hidden executable next to the source whenever that executable is missing or older than the source, runs it without arguments, and parse_output reads the result: it would prefer a value after answer: or after an equals sign, but the C++ program prints only the count, so the last non-empty line is returned. The Python result is therefore the C++ result, checkpoints included; a failed checkpoint makes the C++ program exit with status $2$, and the Python call then raises an error.

Complexity and Verification

For $L=110\,000\,000$ we have $U=u_{\max}=16\,224\,755$. The linear sieve takes $O(U)$ time, and its table of $U+1$ four-byte integers, about 65 MB, is the largest data structure, so the memory is $O(L)$; the trial-division list holds the $1375$ primes up to $11\,393$. For each $u$, factoring $u$ costs $O(\log u)$, and trial division of $m$ makes at most $1375$ tests, so in the worst case trial division costs $O(U\sqrt L/\log L)=O(L^{3/2}/\log L)$. In practice it stops much earlier, because the loop ends once $p^2$ exceeds the unfactored part: over all $u$ it makes about $4.2\cdot10^9$ tests, about $260$ per $u$. With $\tau(u)$ the number of divisors of $u$ and $s(m)$ the number of $d$ with $d^2\mid m$, the divisor loops make at most $\sum_{u\le U}\tau(u)\,s(8u-3)$ iterations, each with at most one gcd of $O(\log L)$ steps. At the target the divisor lists of all $u$ contain $271\,869\,755$ entries, in line with Dirichlet's estimate $U\ln U+(2\gamma-1)U$${}\approx2.72\cdot10^8$, where $\gamma$ is the Euler–Mascheroni constant, ^[10] and the inner loop tests about $2.4\cdot10^8$ pairs $(b_1,h)$, an order of magnitude less work than the trial division. Compiled with -O2 on a current desktop computer, the C++ program finishes within a few seconds.

Before the main count the C++ program checks the limits $8$, $100$ and $1000$. For this article we also verified:

A brute force over $a$ and $b$, taking $c$ from the equation of Section 2, gives $1$, $11$ and $149$ triplets for these limits, $807$ for $L=5000$ and $3312$ for $L=20\,000$. For all $4\,455\,100$ triples with $a+b+c\le300$, we evaluated both real cube roots in double precision and accepted a sum that differs from $1$ by less than $10^{-9}$. Exactly the $40$ triplets of the equation were accepted: their errors are below $10^{-15}$, while every other triple misses $1$ by more than $2\cdot10^{-6}$. This is a numerical cross-check of the Claim; Sections 2 to 4 give the proof.
The compiled C++ program, run with --limit=, agrees with the brute force for every $L\le1200$ and for $510$ further limits up to $20\,000$.
A Python transcription of the enumeration over $(u,h,b_1)$ produces exactly the brute-force list for $L=5000$, with no triplet counted twice.
An independent program agrees with the C++ program at $L=10^6$, $10^7$, $3\cdot10^7$ and $110\,000\,000$. It enumerates all $b$ with $b^2\mid u^2(8u-3)$ directly, without the gcd split or the window for $h$, uses the exact integer cutoff of Section 8, and factors $8u-3$ with a smallest-prime-factor table instead of trial division.
The Java program and the Python bridge print the same value as the C++ program, and the floating-point cutoff of max_u_bound equals the exact integer cutoff at the target and for every $L\le200\,000$.

These additional checks are not part of the published solutions.

Footnotes and References

The references below supply the background facts used above; the reduction to $b^2c=u^2(8u-3)$, the splitting of $b$, the bounds and all counts are derived or computed in this article.

Project Euler 251 — Cardano Triplets. The official statement calls a triple of positive integers (a, b, c) a Cardano triplet if the cube roots of a + b√c and a − b√c add up to 1, gives the example (2, 1, 5), states that 149 Cardano triplets have a + b + c ≤ 1000 and asks how many have a + b + c ≤ 110 000 000.
Wikipedia — Cube root. Every real number has exactly one real cube root, and the cube root of a negative number is negative: the real cube roots of 8 and −8 are 2 and −2. For negative numbers the real cube root is not the principal cube root, which is complex.
cppreference — std::cbrt. The function returns the cube root of its argument, for example −7 for −343. The page notes that std::cbrt(num) is not equivalent to std::pow(num, 1.0 / 3), because std::pow cannot raise a negative base to a fractional exponent; its example prints −nan for pow(−343, 1.0 / 3).
Wikipedia — Cubic equation. If t³ + pt + q = 0 has real coefficients with q²/4 + p³/27 > 0, Cardano's formula gives its real root as the sum of the cube roots of −q/2 + √(q²/4 + p³/27) and −q/2 − √(q²/4 + p³/27), and the other two roots are non-real. The article attributes the method to Scipione del Ferro and Tartaglia and its first publication to Cardano's Ars Magna in 1545.
Wikipedia — Euclid's lemma. The generalisation states that if an integer n divides a product ab and is coprime with a, then n divides b. Section 6 applies it to b₁², which divides h²m and is coprime with h².
Wikipedia — Inequality of arithmetic and geometric means. The arithmetic mean of a list of non-negative real numbers is at least their geometric mean, with equality if and only if all the numbers are equal. Section 8 applies it to b/2, b/2 and c.
C++ — Euler251.cpp. The linked immutable C++ revision contains sieve_primes, build_spf, factorize_square_part, max_u_bound, solve and the checkpoints for the limits 8, 100 and 1000. It is the implementation record for this article.
cp-algorithms — Linear sieve. The linear sieve finds the smallest prime factor of every number up to n, together with the list of primes, in O(n) time, and the stored factors give the factorisation of any number in that range in time proportional to the length of the factorisation. The page compares it with the classic sieve of Eratosthenes, which takes O(n log log n) time.
Wikipedia — Trial division. Trial division tests candidate factors in increasing order, the effort is reduced by using only primes as candidates, and the candidates need only go up to √n, because a composite n has a factor no larger than √n.
Wikipedia — Divisor summatory function. The sum of the numbers of divisors of all n ≤ x is x log x + x(2γ − 1) + O(√x), a result first established by Dirichlet in 1849; the complexity section uses it for the total length of the divisor lists.

Problem 251 source code

C++

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

namespace {

using i64 = std::int64_t;
using i128 = __int128_t;

struct Options {
    i64 limit = 110000000;
    int threads = 0;
    bool run_checkpoints = true;
};

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

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

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

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

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

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

std::vector<int> sieve_primes(const int limit) {
    std::vector<bool> is_prime(static_cast<std::size_t>(limit + 1), true);
    is_prime[0] = false;
    if (limit >= 1) {
        is_prime[1] = false;
    }

    for (int i = 2; static_cast<i64>(i) * i <= limit; ++i) {
        if (!is_prime[static_cast<std::size_t>(i)]) {
            continue;
        }
        for (int j = i * i; j <= limit; j += i) {
            is_prime[static_cast<std::size_t>(j)] = false;
        }
    }

    std::vector<int> primes;
    for (int i = 2; i <= limit; ++i) {
        if (is_prime[static_cast<std::size_t>(i)]) {
            primes.push_back(i);
        }
    }
    return primes;
}

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

    for (int i = 2; i <= n; ++i) {
        if (spf[static_cast<std::size_t>(i)] == 0) {
            spf[static_cast<std::size_t>(i)] = i;
            primes.push_back(i);
        }
        for (int p : primes) {
            const i64 v = static_cast<i64>(p) * static_cast<i64>(i);
            if (v > n || p > spf[static_cast<std::size_t>(i)]) {
                break;
            }
            spf[static_cast<std::size_t>(v)] = p;
        }
    }

    return spf;
}

void factorize_with_spf(int x, const std::vector<int>& spf, std::vector<std::pair<int, int>>& factors) {
    factors.clear();
    while (x > 1) {
        const int p = spf[static_cast<std::size_t>(x)];
        int e = 0;
        do {
            x /= p;
            ++e;
        } while (x > 1 && spf[static_cast<std::size_t>(x)] == p);
        factors.push_back({p, e});
    }
}

void factorize_square_part(i64 m, const std::vector<int>& primes, std::vector<std::pair<int, int>>& square_factors) {
    square_factors.clear();
    i64 x = m;

    for (int p : primes) {
        if (static_cast<i64>(p) * static_cast<i64>(p) > x) {
            break;
        }
        if ((x % p) != 0) {
            continue;
        }

        int e = 0;
        do {
            x /= p;
            ++e;
        } while ((x % p) == 0);

        if (e >= 2) {
            square_factors.push_back({p, e / 2});
        }
    }
}

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

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

i64 floor_sqrt_i64(const i64 x) {
    if (x <= 0) {
        return 0;
    }
    i64 r = static_cast<i64>(std::sqrt(static_cast<long double>(x)));
    while (static_cast<i128>(r + 1) * static_cast<i128>(r + 1) <= x) {
        ++r;
    }
    while (static_cast<i128>(r) * static_cast<i128>(r) > x) {
        --r;
    }
    return r;
}

bool can_have_solution(const i64 u, const i64 limit) {
    const long double uu = static_cast<long double>(u);
    const long double mm = 8.0L * uu - 3.0L;
    const long double lower_bound = 3.0L * uu - 1.0L + 3.0L * std::cbrt((uu * uu * mm) / 4.0L);
    return lower_bound <= static_cast<long double>(limit) + 1e-9L;
}

i64 max_u_bound(const i64 limit) {
    i64 lo = 0;
    i64 hi = (limit + 1) / 3;
    while (lo < hi) {
        const i64 mid = lo + (hi - lo + 1) / 2;
        if (can_have_solution(mid, limit)) {
            lo = mid;
        } else {
            hi = mid - 1;
        }
    }
    return lo;
}

i64 solve(const i64 limit, int threads) {
    if (limit < 8) {
        return 0;
    }

    const i64 max_u = max_u_bound(limit);
    if (max_u <= 0) {
        return 0;
    }

    const int sqrt_max_m = static_cast<int>(std::sqrt(static_cast<long double>(8 * max_u - 3))) + 1;
    const std::vector<int> primes = sieve_primes(sqrt_max_m);
    const std::vector<int> spf = build_spf(static_cast<int>(max_u));

    if (threads <= 0) {
        threads = static_cast<int>(std::thread::hardware_concurrency());
        if (threads <= 0) {
            threads = 1;
        }
    }
    threads = std::max(1, std::min(threads, static_cast<int>(max_u)));

    std::vector<i64> local_counts(static_cast<std::size_t>(threads), 0);
    std::vector<std::thread> workers;
    workers.reserve(static_cast<std::size_t>(threads));

    for (int t = 0; t < threads; ++t) {
        workers.emplace_back([&, t]() {
            std::vector<std::pair<int, int>> factors_u;
            std::vector<std::pair<int, int>> square_factors_m;
            std::vector<i64> divisors_h;
            std::vector<i64> divisors_b1;
            factors_u.reserve(16);
            square_factors_m.reserve(8);
            divisors_h.reserve(256);
            divisors_b1.reserve(64);

            i64 subtotal = 0;
            for (i64 u = static_cast<i64>(t) + 1; u <= max_u; u += threads) {
                const i64 a = 3 * u - 1;
                const i64 remaining = limit - a;
                if (remaining <= 1) {
                    continue;
                }

                factorize_with_spf(static_cast<int>(u), spf, factors_u);
                generate_divisors(factors_u, divisors_h);

                const i64 m = 8 * u - 3;
                factorize_square_part(m, primes, square_factors_m);
                generate_divisors(square_factors_m, divisors_b1);

                for (const i64 b1 : divisors_b1) {
                    const i64 b1_sq = b1 * b1;
                    const i64 m_over = m / b1_sq;

                    if (b1 >= remaining || b1 + m_over > remaining) {
                        continue;
                    }

                    const i64 h_max = floor_sqrt_i64((remaining - b1) / m_over);
                    if (h_max <= 0) {
                        continue;
                    }

                    i64 h_min = 1;
                    const i64 denom = remaining - m_over;
                    if (denom > 0) {
                        h_min = (u * b1 + denom - 1) / denom;
                    }

                    for (const i64 h : divisors_h) {
                        if (h < h_min || h > h_max) {
                            continue;
                        }
                        if (std::gcd(h, b1) != 1) {
                            continue;
                        }

                        const i64 b = (u / h) * b1;
                        const i64 c = h * h * m_over;
                        if (static_cast<i128>(b) + static_cast<i128>(c) <= remaining) {
                            ++subtotal;
                        }
                    }
                }
            }
            local_counts[static_cast<std::size_t>(t)] = subtotal;
        });
    }

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

    i64 total = 0;
    for (const i64 v : local_counts) {
        total += v;
    }
    return total;
}

bool run_checkpoints() {
    if (solve(8, 1) != 1) {
        std::cerr << "Checkpoint failed for limit=8" << '\n';
        return false;
    }
    if (solve(100, 1) != 11) {
        std::cerr << "Checkpoint failed for limit=100" << '\n';
        return false;
    }
    if (solve(1000, 1) != 149) {
        std::cerr << "Checkpoint failed for limit=1000" << '\n';
        return false;
    }
    return true;
}

}  // namespace

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

    std::cout << solve(options.limit, options.threads) << '\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 ""

    answer_candidates = []
    equal_candidates = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answer_candidates.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equal_candidates.append(m2.group(1).strip())

    if answer_candidates:
        return answer_candidates[-1]
    if equal_candidates:
        return equal_candidates[-1]
    return lines[-1]


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 = subprocess.check_output([str(binary)], text=True)
    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.List;
import java.util.concurrent.*;

public class Euler251 {
    static final long LIMIT = 110000000L;

    static List<Integer> sievePrimes(int limit) {
        boolean[] isPrime = new boolean[limit + 1];
        for (int i = 2; i <= limit; i++)
            isPrime[i] = true;

        for (int i = 2; (long) i * i <= limit; ++i) {
            if (!isPrime[i])
                continue;
            for (int j = i * i; j <= limit; j += i) {
                isPrime[j] = false;
            }
        }

        List<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= limit; ++i) {
            if (isPrime[i])
                primes.add(i);
        }
        return primes;
    }

    static int[] buildSPF(int n) {
        int[] spf = new int[n + 1];
        if (n >= 1)
            spf[1] = 1;
        List<Integer> primes = new ArrayList<>();

        for (int i = 2; i <= n; ++i) {
            if (spf[i] == 0) {
                spf[i] = i;
                primes.add(i);
            }
            for (int p : primes) {
                long v = (long) p * i;
                if (v > n || p > spf[i])
                    break;
                spf[(int) v] = p;
            }
        }
        return spf;
    }

    static class Factor {
        int p, e;

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

    static void factorizeWithSPF(int x, int[] spf, List<Factor> factors) {
        factors.clear();
        while (x > 1) {
            int p = spf[x];
            int e = 0;
            do {
                x /= p;
                e++;
            } while (x > 1 && spf[x] == p);
            factors.add(new Factor(p, e));
        }
    }

    static void factorizeSquarePart(long m, List<Integer> primes, List<Factor> squareFactors) {
        squareFactors.clear();
        long x = m;
        for (int p : primes) {
            if ((long) p * p > x)
                break;
            if (x % p != 0)
                continue;
            int e = 0;
            do {
                x /= p;
                e++;
            } while (x % p == 0);
            if (e >= 2) {
                squareFactors.add(new Factor(p, e / 2));
            }
        }
    }

    static void generateDivisors(List<Factor> factors, List<Long> divisors) {
        divisors.clear();
        divisors.add(1L);
        for (Factor f : factors) {
            int baseSize = divisors.size();
            long pe = 1;
            for (int k = 1; k <= f.e; ++k) {
                pe *= f.p;
                for (int i = 0; i < baseSize; ++i) {
                    divisors.add(divisors.get(i) * pe);
                }
            }
        }
    }

    static long floorSqrt(long x) {
        if (x <= 0)
            return 0;
        long r = (long) Math.sqrt(x);
        while (BigIntegerUtil.multiply(r + 1, r + 1).compareTo(BigIntegerUtil.valueOf(x)) <= 0)
            r++;
        while (BigIntegerUtil.multiply(r, r).compareTo(BigIntegerUtil.valueOf(x)) > 0)
            r--;
        return r;
    }

    static boolean canHaveSolution(long u, long limit) {
        double uu = (double) u;
        double mm = 8.0 * uu - 3.0;
        double lowerBound = 3.0 * uu - 1.0 + 3.0 * Math.cbrt((uu * uu * mm) / 4.0);
        return lowerBound <= (double) limit + 1e-9;
    }

    static long maxUBound(long limit) {
        long lo = 0;
        long hi = (limit + 1) / 3;
        while (lo < hi) {
            long mid = lo + (hi - lo + 1) / 2;
            if (canHaveSolution(mid, limit)) {
                lo = mid;
            } else {
                hi = mid - 1;
            }
        }
        return lo;
    }

    static long gcd(long a, long b) {
        while (b != 0) {
            long t = b;
            b = a % b;
            a = t;
        }
        return a;
    }

    public static String solve() {
        if (LIMIT < 8)
            return "0";

        long maxU = maxUBound(LIMIT);
        if (maxU <= 0)
            return "0";

        int sqrtMaxM = (int) Math.sqrt(8.0 * maxU - 3.0) + 1;
        List<Integer> primes = sievePrimes(sqrtMaxM);
        int[] spf = buildSPF((int) maxU);

        int threads = Math.max(1, Runtime.getRuntime().availableProcessors());
        threads = (int) Math.min(threads, Math.max(1, maxU));

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

        for (int t = 0; t < threads; ++t) {
            final int tId = t;
            final int tCount = threads;
            futures.add(executor.submit(() -> {
                List<Factor> factorsU = new ArrayList<>(16);
                List<Factor> squareFactorsM = new ArrayList<>(8);
                List<Long> divisorsH = new ArrayList<>(256);
                List<Long> divisorsB1 = new ArrayList<>(64);

                long subtotal = 0;
                for (long u = tId + 1; u <= maxU; u += tCount) {
                    long a = 3 * u - 1;
                    long remaining = LIMIT - a;
                    if (remaining <= 1)
                        continue;

                    factorizeWithSPF((int) u, spf, factorsU);
                    generateDivisors(factorsU, divisorsH);

                    long m = 8 * u - 3;
                    factorizeSquarePart(m, primes, squareFactorsM);
                    generateDivisors(squareFactorsM, divisorsB1);

                    for (long b1 : divisorsB1) {
                        long b1Sq = b1 * b1;
                        long mOver = m / b1Sq;

                        if (b1 >= remaining || b1 + mOver > remaining)
                            continue;

                        long hMax = floorSqrt((remaining - b1) / mOver);
                        if (hMax <= 0)
                            continue;

                        long hMin = 1;
                        long denom = remaining - mOver;
                        if (denom > 0) {
                            hMin = (u * b1 + denom - 1) / denom;
                        }

                        for (long h : divisorsH) {
                            if (h < hMin || h > hMax)
                                continue;
                            if (gcd(h, b1) != 1)
                                continue;

                            long b = (u / h) * b1;
                            long c = h * h * mOver;
                            if (BigIntegerUtil.valueOf(b).add(BigIntegerUtil.valueOf(c))
                                    .compareTo(BigIntegerUtil.valueOf(remaining)) <= 0) {
                                subtotal++;
                            }
                        }
                    }
                }
                return subtotal;
            }));
        }

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

        return String.valueOf(total);
    }

    static class BigIntegerUtil {
        static java.math.BigInteger valueOf(long val) {
            return java.math.BigInteger.valueOf(val);
        }

        static java.math.BigInteger multiply(long a, long b) {
            return valueOf(a).multiply(valueOf(b));
        }
    }

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