Problem 223: Almost Right-angled Triangles I

Project Euler Problem 223 Solution

EulerSolve provides an optimized solution for Project Euler Problem 223, Almost Right-angled Triangles I, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We want to count integer-sided triangles $(a,b,c)$ with $a \le b \le c$, perimeter $a+b+c \le P$, and $$a^2+b^2=c^2+1.$$ For the Project Euler instance, $P=25{,}000{,}000$. The equation says that the triangle is exactly one unit away from being right-angled, so a brute-force search over all triples is completely impractical. The successful idea is to fix the smallest side $a$, rewrite the equation as a factorization of $a^2-1$, and count only the divisor pairs that can still satisfy the perimeter and ordering constraints. Mathematical Approach The central observation is that once $a$ is fixed, the other two sides can be recovered from a factor pair of $a^2-1$. That turns the problem from a quadratic search in $(b,c)$ into a bounded divisor-counting problem. Fix the smallest side first Because $a \le b \le c$, every valid triangle satisfies $$3a \le a+b+c \le P.$$ Therefore $$1 \le a \le \left\lfloor \frac{P}{3}\right\rfloor.$$ The implementations loop over this range of possible smallest sides and solve the rest of the triangle from there....

Detailed mathematical approach

Problem Summary

A triangle with integer sides $a\le b\le c$ is called barely acute if $a^2+b^2=c^2+1$. We must count the barely acute triangles whose perimeter $a+b+c$ is at most $25\,000\,000$, counting each triangle once, with its sides sorted. ^[1] The statement gives no checkpoint values. We write $N(P)$ for the number of barely acute triangles with perimeter at most $P$; the task asks for $N(25\,000\,000)$, and the smaller values of $N(P)$ quoted below were computed for this article.

The whole solution in one sentence: apart from the triangles $(1,b,b)$, the equation says $(c-b)(c+b)=a^2-1$, so the triangles with smallest side $a$ correspond to the divisors of $a^2-1=(a-1)(a+1)$ that lie in a short interval and have the right parity, and a table of smallest prime factors lets us list these divisors quickly for every $a$ up to $P/3$.

Mathematical Approach

We look at small triangles first and then turn the triangles with a fixed smallest side into divisors of a single number.

1. Small cases

A short search lists all barely acute triangles with perimeter at most $50$. There are $32$ of them: the $24$ triangles $(1,1,1)$, $(1,2,2)$, …, $(1,24,24)$, and the eight triangles $(5,5,7)$, $(4,7,8)$, $(8,9,12)$, $(7,11,13)$, $(6,17,18)$, $(11,13,17)$, $(10,15,18)$ and $(9,19,21)$, with the perimeters $17$, $19$, $29$, $31$, $41$, $41$, $43$ and $49$. For example, $8^2+9^2=145=12^2+1$. Three things stand out: the smallest side $1$ gives a whole family, no triangle has the smallest side $2$ or $3$, and for $a\ge2$ the longest side is always longer than the middle one. The smallest side $8$ also has the triangle $(8,31,32)$, but its perimeter $71$ is too large here.

2. Every solution is an acute triangle

Three positive lengths form a triangle exactly when the longest one is shorter than the sum of the other two. Here this holds automatically:

$$(a+b)^2=a^2+b^2+2ab=c^2+1+2ab\gt c^2,$$

so $a+b\gt c$, and every solution of the equation with $a\le b\le c$ is a genuine triangle. The law of cosines, $c^2=a^2+b^2-2ab\cos\gamma$ with $\gamma$ the angle opposite the longest side, ^[2] gives $\cos\gamma=1/(2ab)$. This is positive, so the largest angle $\gamma$, and with it every angle, is acute. As $ab$ grows, $\cos\gamma$ tends to $0$ and $\gamma$ approaches $90^\circ$ from below: for $(8,9,12)$ we get $\cos\gamma=1/144$ and $\gamma\approx89.60^\circ$, and for $(17,71,73)$ already $\gamma\approx89.98^\circ$, while the equilateral triangle $(1,1,1)$ has $\cos\gamma=1/2$. The fast method therefore never has to test the triangle inequality.

3. The family with smallest side 1

For $a=1$ the equation becomes $1+b^2=c^2+1$, so $b=c$, and every $b\ge1$ gives the barely acute triangle $(1,b,b)$. Its perimeter $1+2b$ is at most $P$ exactly when $b\le(P-1)/2$, so this family contributes

$$\left\lfloor\frac{P-1}{2}\right\rfloor$$

triangles, $24$ for $P=50$. It contains every triangle with $b=c$, because for $a\ge2$ the equation gives $c^2-b^2=a^2-1\gt0$, so $c\gt b$.

4. A difference of two squares

Fix $a\ge2$ and write the equation as $c^2-b^2=a^2-1$. A difference of two squares factors as $c^2-b^2=(c-b)(c+b)$, ^[3] so with

$$v=c-b,\qquad u=c+b$$

we get two positive integers with $uv=a^2-1$. Their sum $u+v=2c$ is even, so $u$ and $v$ have the same parity. Conversely, if $a^2-1=uv$ with positive integers $v\lt u$ of the same parity, then $b=(u-v)/2$ and $c=(u+v)/2$ are positive integers with $c^2-b^2=uv=a^2-1$. The factor $v$ alone determines the triangle, because $u=(a^2-1)/v$. For $a=8$ the number $63$ has the factor pairs $1\cdot63$, $3\cdot21$ and $7\cdot9$, all odd, and they give $(b,c)=(31,32)$, $(9,12)$ and $(1,8)$. The last pair is the triangle $(1,8,8)$ of Section 3 with its sides in the wrong order, because $b=1\lt a$: the ordering $a\le b$ is a real condition.

5. The two conditions become an interval

Since $v\ge1$, the inequality $c=b+v\gt b$ holds automatically. The ordering $a\le b$ means $u-v\ge2a$. Multiplying by $v\gt0$ and using $uv=a^2-1$ turns it into $a^2-1\ge v^2+2av$, that is $(v+a)^2\le2a^2-1$. As $v$ is an integer, this is

$$v\le v_{\max}=\left\lfloor\sqrt{2a^2-1}\right\rfloor-a.$$

The perimeter is $a+b+c=a+u$, so the bound $a+u\le P$ means $(a^2-1)/v\le P-a$. Only $a\le P/3$ can occur, because $3a\le a+b+c$, so $P-a$ is positive, and the bound becomes

$$v\ge v_{\min}=\left\lceil\frac{a^2-1}{P-a}\right\rceil.$$

Every step above is an equivalence, and together with Section 4 this proves the key claim. Claim. For $a\ge2$, the barely acute triangles with smallest side $a$ and perimeter at most $P$ correspond one to one to the divisors $v$ of $a^2-1$ with $v_{\min}\le v\le v_{\max}$ for which $v$ and $(a^2-1)/v$ have the same parity. By Section 2 each of them is a genuine triangle. For $a=8$ we get $v_{\max}=\lfloor\sqrt{127}\rfloor-8=3$, which removes the pair $7\cdot9$. For $a=2$ we get $v_{\max}=\lfloor\sqrt7\rfloor-2=0$, and for $a=3$ the only candidate $v=1$ has the cofactor $8$ of the other parity; this is why Section 1 found no triangle with smallest side $2$ or $3$. Since $\sqrt{2a^2-1}\lt\sqrt2\,a$, we also have $v_{\max}\lt(\sqrt2-1)\,a\lt0.42\,a$: only small divisors matter.

6. Parity

If $a$ is even, $a^2-1$ is odd, so all its divisors are odd and the parity condition holds by itself. If $a$ is odd, $a^2-1=a^2-1^2$ is a difference of two odd squares and therefore a multiple of $8$; ^[3] indeed $a-1$ and $a+1$ are consecutive even numbers, and one of them is a multiple of $4$. Two odd numbers cannot have an even product, so for odd $a$ both $v$ and $(a^2-1)/v$ must be even.

7. Worked example: a = 17 and P = 200

Here $a^2-1=288=16\cdot18=2^5\cdot3^2$, which has $18$ divisors. The interval of Section 5 is given by

$$\begin{aligned}v_{\min}&=\left\lceil\frac{288}{183}\right\rceil=2,\\v_{\max}&=\lfloor\sqrt{577}\rfloor-17=7.\end{aligned}$$

The divisors up to $7$ are $1$, $2$, $3$, $4$ and $6$. The divisor $1$ lies below the interval: its cofactor $288$ would give the perimeter $17+288=305$. The divisor $3$ is odd while its cofactor $96$ is even. The divisors $2$, $4$ and $6$ have the even cofactors $144$, $72$ and $48$, and they give

$$(17,71,73),\qquad(17,34,38),\qquad(17,21,27),$$

with the perimeters $161$, $89$ and $65$. For example, $17^2+21^2=730=27^2+1$. For $P=100$ the interval starts at $\lceil288/83\rceil=4$, and only the last two triangles remain. In total $N(100)=72$ and $N(200)=167$, and $49$ and $99$ of these triangles belong to the family $(1,b,b)$.

8. How large the smallest side can be

The bound $a\le P/3$ of Section 5 is far from tight. Since $v\le\sqrt{2a^2-1}-a$ and

$$\bigl(\sqrt{2a^2-1}-a\bigr)\bigl(\sqrt{2a^2-1}+a\bigr)=a^2-1,$$

we get

$$u=\frac{a^2-1}{v}\ge\sqrt{2a^2-1}+a,$$

so every triangle with smallest side $a$ has a perimeter of at least $2a+\sqrt{2a^2-1}$. As $\sqrt{2a^2-1}\gt\sqrt2\,a-1$ for $a\ge1$, only smallest sides with $(2+\sqrt2)\,a-1\lt P$ can occur, that is $a\lt(P+1)/(2+\sqrt2)\approx0.293\,(P+1)$. For $P=25\,000\,000$ the interval of Section 5 is nonempty exactly for $3\le a\le7\,322\,328$, while $\lfloor P/3\rfloor=8\,333\,333$.

Equality in this perimeter bound holds exactly when $v=\sqrt{2a^2-1}-a$, that is, when $u-v=2a$ and so $b=a$. For a given $a$, such a triangle exists exactly when $2a^2-1$ is a perfect square $s^2$, and it is then $(a,a,s)$; the other triangles with the same smallest side, such as $(29,37,47)$ next to $(29,29,41)$, have larger perimeters. The equation $s^2-2a^2=-1$ is the negative Pell equation for $2$, and the odd powers of its fundamental solution $1+\sqrt2$ give an infinite tower of solutions, such as $(1+\sqrt2)^3=7+5\sqrt2$ and $(1+\sqrt2)^5=41+29\sqrt2$. ^[4] The numbers $a$ for which $2a^2-1$ is a square are $1,5,29,169,985,5741,\dots$, ^[5] so the barely acute triangles with $a=b$ are $(1,1,1)$, $(5,5,7)$, $(29,29,41)$, $(169,169,239)$, $(985,985,1393)$ and so on. For these $a$ the bound $v_{\max}$ is attained exactly, so a square root that comes out slightly too small would lose a triangle; the code avoids this with exact integer corrections.

9. Factoring a² − 1 with a sieve

To list divisors we need the prime factorization of $a^2-1=(a-1)(a+1)$. Both factors are at most $\lfloor P/3\rfloor+1$. A common divisor of them divides their difference $2$, so they share at most the prime $2$, and they share it exactly when $a$ is odd. Their factorizations can therefore be merged: a prime of only one factor is copied, and for the prime $2$ of an odd $a$ the two exponents are added. For $a=17$, $16=2^4$ and $18=2\cdot3^2$ give $2^5\cdot3^2$.

The factorizations come from a table of smallest prime factors, where $\operatorname{spf}(m)$ is the smallest prime that divides $m$. Dividing $m$ by $\operatorname{spf}(m)$ again and again lists the prime factors of $m$ in increasing order, with at most $\log_2m$ divisions; for $18$ we get $18\to9\to3\to1$ and the primes $2$, $3$, $3$. The table is filled by the linear sieve. A loop over $i$ records every $i$ without an entry as a new prime, and then sets $\operatorname{spf}(p\,i)=p$ for every prime $p\le\operatorname{spf}(i)$ with $p\,i$ inside the table. Every composite $m$ is written exactly once, as $p\cdot(m/p)$ with $p=\operatorname{spf}(m)$, so the sieve takes time proportional to the size of the table. ^[6]

10. Listing only the small divisors

If $a^2-1=p_1^{e_1}\cdots p_r^{e_r}$, it has $\tau(a^2-1)=(e_1+1)\cdots(e_r+1)$ divisors, ^[7] and they can be listed prime by prime: start with the list $1$, and for each prime $p_k$ multiply every divisor listed so far by $p_k,p_k^2,\dots,p_k^{e_k}$ and append the products. By unique factorization every divisor appears exactly once. Since only divisors up to $v_{\max}$ are needed, a product larger than $v_{\max}$ is not appended. Nothing is lost this way: if $d\le v_{\max}$ is a divisor whose largest prime is $p_k$, with exponent $j$, then $d/p_k^{\,j}\le d\le v_{\max}$ was kept earlier, and $d$ is formed from it.

For $a=17$ and $v_{\max}=7$, the prime $2$ turns the list $1$ into $1,2,4$, since $8$, $16$ and $32$ are too large, and the prime $3$ adds $3$ and $6$ but not $12$, $9$, $18$ or $36$. The final list $1,2,4,3,6$ holds exactly the divisors of Section 7. Every product formed is a divisor, so at most $\tau(a^2-1)-1$ products are computed, and fewer in practice, because a divisor that was dropped is never multiplied again.

11. The complete count

Adding the family of Section 3 to the divisor counts of Section 5 gives

$$N(P)=\left\lfloor\frac{P-1}{2}\right\rfloor+\sum_{a=2}^{\lfloor P/3\rfloor}D_a(P),$$

where $D_a(P)$ is the number of divisors $v$ described in the claim of Section 5. For $P=50$ the eight values of $a$ from $4$ to $11$ contribute one triangle each, and $24+8=32$. For $P=10^2,10^3,\dots,10^6$ the formula gives $N(P)=72$, $1046$, $13\,656$, $168\,921$ and $2\,012\,468$. The family $(1,b,b)$ supplies $49$, $499$, $4999$, $49\,999$ and $499\,999$ of them, so the other triangles, $23$ of $72$ at first, make up about three quarters at $P=10^6$.

12. A different route: a tree of triangles

The published programs do not use it, but a second generator gives an independent check. The second of Berggren's three matrices for the tree of primitive Pythagorean triples ^[8] sends $(a,b,c)$ to

$$(a+2b+2c,\;2a+b+2c,\;2a+2b+3c),$$

and the first and the third do the same after changing the sign of $b$ and of $a$, respectively. Expanding the squares shows that all three preserve $a^2+b^2-c^2$, and each of them maps a triangle with $a\le b\le c$ to one with positive sides and a larger perimeter; for example, they send $(1,1,1)$ to $(1,3,3)$, $(5,5,7)$ and $(3,1,3)$. A separate program written for this article starts from $(1,1,1)$ and $(1,2,2)$, applies the three matrices, puts the two smaller sides in order and, at a triangle with $a=b$, skips the third matrix, whose result is then the mirror image of the first. For $P\le10^6$ it produces every barely acute triangle exactly once. This is a finite check only: preserving the equation does not prove that the generator reaches every triangle for every $P$.

How the Code Works

The C++ program ^[9] starts in main with validate, described in the next section, and then counts for $P=25\,000\,000$, or for a bound given as the first command-line argument. count_barely_acute sets $a_{\max}=\lfloor P/3\rfloor$, takes the family of Section 3 as a1_count, $\lfloor(P-1)/2\rfloor$, and builds the table of smallest prime factors up to $a_{\max}+1$ with build_spf, the linear sieve of Section 9, which also collects the primes in a list and sets spf[1] to $1$. The values $a=2,\dots,a_{\max}$ are then handed out in consecutive chunks of $2048$, every value is counted exactly once, and the chunk counts are added to a1_count in unsigned 64-bit integers.

count_for_a handles one value of $a$. It forms $n=a^2-1$ in unsigned 64-bit arithmetic, computes $v_{\min}$ by the ceiling division $(n+P-a-1)/(P-a)$, and estimates $v_{\max}$ as $\lfloor\sqrt{2a^2-1}-a\rfloor$ with a long double square root. An estimate of $0$ or less ends the call at once; for $a\le8\,333\,333$ this happens only for $a=2$, where $v_{\max}=0$ is exact. Otherwise two loops correct the estimate by exact comparisons of $v^2+2av$ with $n$, the inequality of Section 5, in the 128-bit integer type __int128, a compiler extension: ^[10] they raise it while $v_{\max}+1$ still satisfies the inequality and lower it while $v_{\max}$ does not. If then $v_{\max}\lt v_{\min}$, the call returns before any factoring (Section 8). Otherwise factorize_int factors $a-1$ and $a+1$ with the table, merge_factors merges the two sorted lists of prime–exponent pairs and adds the exponents of a common prime (Section 9), and the divisor loop of Section 10 fills the vector divisors with the divisors up to $v_{\max}$. Every product formed there divides $n$, so it stays below $7\cdot10^{13}$ and fits easily into 64 bits.

The last loop tests each listed divisor $v$: it skips $v\lt v_{\min}$, computes $u=n/v$, skips the divisors with an odd $u-v$ (Section 6), and finally tests $u\ge v+2a$ and $a+u\le P$. By Section 5 these two tests are equivalent to $v\le v_{\max}$ and $v\ge v_{\min}$, so they never reject a divisor; they only guard against a mistake in the bounds. Every divisor that passes is one triangle.

The Python program has the same structure. get_spf, factorize_int and merge_factors match their C++ counterparts, get_divisors is the pruned divisor loop, and worker carries out the steps of count_for_a for every $a$ of one chunk, with math.sqrt for the estimate and Python's exact integers for the corrections. solve cuts $2\le a\le\lfloor P/3\rfloor$ into the same chunks of $2048$ values, adds the results of worker to $\lfloor(P-1)/2\rfloor$ and returns the count as a string. The Java program inlines count_for_a into solve, which counts every chunk in a task of its own. It keeps the divisors in a long array of $256$ entries that doubles when it is full, and it does the corrections in long arithmetic, which is safe because $v^2+2av$ stays far below $2^{63}$ here. Two differences matter only when something goes wrong: neither Python nor Java runs the self-checks of the C++ program, and Java's solve silently ignores an exception thrown by a task, whose chunk would then be missing from the total.

Complexity and Verification

The sieve fills a table of $\lfloor P/3\rfloor+2=8\,333\,335$ entries in time proportional to its size, and while it runs it also keeps the list of the $560\,690$ primes up to $8\,333\,334$. The table needs about $33$ MB as 32-bit integers in C++ and Java, and more than twice as much in CPython, whose list holds 8-byte references to integer objects. Then every $a\le P/3$ costs a square root and a few comparisons. A value of $a$ with a nonempty interval also needs the two factorizations, whose $e_1+\dots+e_r$ divisions in all are fewer than $\tau(a^2-1)=(e_1+1)\cdots(e_r+1)$, and the divisor loop, which forms at most $\tau(a^2-1)-1$ products. So the running time is

$$O\Bigl(P+\sum_{a}\tau\bigl(a^2-1\bigr)\Bigr),$$

where $a$ runs over the values with a nonempty interval, at the target the $7\,322\,326$ values $3\le a\le7\,322\,328$. Apart from the table, only small buffers and, in Python and Java, one entry for each of the $4070$ chunks are needed. Over these values of $a$ the sum of $\tau(a^2-1)$ is about $1.24\cdot10^9$, roughly $170$ divisors per value. The pruning of Section 10 limits the divisor loop to about $1.08\cdot10^9$ products, and about $5.6\cdot10^8$ divisors reach the final tests. The largest divisor count among these values is $\tau(a^2-1)=9216$, for $a=5\,395\,391$. On the arm64 computer used for this article, with the C++ program compiled by Clang 22 with -O2, Java 17 and Python 3.11, the C++ and Java programs finish in under a second and the Python program in about $11$ seconds.

Before its main computation, the C++ program compares its count with brute_count for $P=80$, $120$, $200$, $300$ and $500$. brute_count tries every pair $a\le b$ with $a\le P/3$ and $b\le(P-a)/2$, takes the exact integer square root of $a^2+b^2-1$ and checks the remaining conditions, including the triangle inequality. It then counts $P=50\,000$ with two different settings of the second argument of count_barely_acute, which only change how the chunks are scheduled, and requires the same result; a mismatch stops the program with an error message. All three programs print the same value. Further checks written for this article, which are not part of the published solutions, confirm the method:

a brute force and the divisor method produce the same set of $4932$ triangles for $P=4000$ and the same count for every $P\le1000$, and a brute force in C gives $N(10^5)=168\,921$;
a separate Python implementation of the divisor method and the tree of Section 12 reproduce the values of Section 11 up to $P=10^6$, and the tree agrees with the C++ program for $P=10^7$ and for $P=25\,000\,000$;
on this computer, where long double has the same 53-bit precision as double, the floating-point estimate of $v_{\max}$ was exact for every $a\le8\,333\,333$, so the correction loops never changed it; they guarantee the right value without relying on this.

Footnotes and References

The references below supply the background facts used above. The correspondence between triangles and divisors, the interval for $v$ and the pruned divisor listing are derived explicitly in this article; the operation totals were measured for it.

Project Euler 223 — Almost Right-angled Triangles I. The official statement calls an integer-sided triangle with sides a ≤ b ≤ c barely acute if a² + b² = c² + 1, and asks how many barely acute triangles have a perimeter of at most 25 000 000.
Wikipedia — Law of cosines. For a triangle with sides a, b, c and the angle γ opposite c, c² = a² + b² − 2ab cos γ; for a right angle, cos γ = 0 and the law reduces to Pythagoras' theorem. Section 2 uses it to obtain cos γ = 1/(2ab).
Wikipedia — Difference of two squares. Every difference of two squares factors as a² − b² = (a + b)(a − b), and the difference of two odd squares is a multiple of 8.
Wikipedia — Pell's equation. The negative Pell equation x² − ny² = −1 is solvable for n = 2, with the fundamental solution (1, 1), and the odd powers of the fundamental solution give an infinite tower of further solutions.
OEIS A001653 — Numbers k such that 2k² − 1 is a square. The sequence begins 1, 5, 29, 169, 985, 5741, 33461, …; by Section 8 these are exactly the values of a for which some barely acute triangle has b = a.
cp-algorithms — Linear Sieve. The algorithm computes the smallest prime factor of every number up to n together with the list of primes. It goes through every composite number exactly once, so it runs in O(n) time, and it needs an array of n entries plus the list of primes.
Wikipedia — Divisor function. If n is the product of the prime powers pᵢ^aᵢ of distinct primes pᵢ, the number of divisors of n is σ₀(n) = ∏ (aᵢ + 1); this article writes τ(n) for σ₀(n).
Wikipedia — Tree of primitive Pythagorean triples. The article gives Berggren's three matrices A, B and C, which turn a Pythagorean triple into a larger one and generate every primitive triple exactly once from (3, 4, 5). Section 12 applies them to barely acute triangles.
C++ — Euler223.cpp. The linked immutable C++ revision contains build_spf, factorize_int, merge_factors, count_for_a, count_barely_acute, brute_count and validate. It is the implementation record for this article.
GCC manual — 128-bit Integers. The types __int128 and unsigned __int128 are an extension, supported for targets that have an integer mode wide enough to hold 128 bits.

Problem 223 source code

C++

#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <thread>
#include <vector>

namespace {

constexpr int kDefaultLimit = 25'000'000;

struct Factor {
    int p;
    int e;
};

std::vector<int> build_spf(int n) {
    std::vector<int> spf(n + 1, 0);
    std::vector<int> primes;
    primes.reserve(n / 10);

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

void factorize_int(int x, const std::vector<int>& spf, std::vector<Factor>& out) {
    out.clear();
    while (x > 1) {
        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);
        out.push_back({p, e});
    }
}

void merge_factors(const std::vector<Factor>& a,
                   const std::vector<Factor>& b,
                   std::vector<Factor>& out) {
    out.clear();
    out.reserve(a.size() + b.size());

    std::size_t i = 0;
    std::size_t j = 0;
    while (i < a.size() || j < b.size()) {
        if (j == b.size() || (i < a.size() && a[i].p < b[j].p)) {
            out.push_back(a[i++]);
        } else if (i == a.size() || b[j].p < a[i].p) {
            out.push_back(b[j++]);
        } else {
            out.push_back({a[i].p, a[i].e + b[j].e});
            ++i;
            ++j;
        }
    }
}

std::uint64_t count_for_a(int a,
                          int perimeter_limit,
                          const std::vector<int>& spf,
                          std::vector<Factor>& left_factors,
                          std::vector<Factor>& right_factors,
                          std::vector<Factor>& merged_factors,
                          std::vector<std::uint64_t>& divisors) {
    const std::uint64_t n = static_cast<std::uint64_t>(a) * static_cast<std::uint64_t>(a) - 1ULL;
    if (n == 0) return 0;

    const std::uint64_t denom = static_cast<std::uint64_t>(perimeter_limit - a);
    if (denom == 0) return 0;

    const std::uint64_t v_min = (n + denom - 1) / denom;

    const long double root = std::sqrt(static_cast<long double>(2) * a * a - 1.0L);
    std::int64_t v_max = static_cast<std::int64_t>(std::floor(root - a));
    if (v_max <= 0) return 0;

    while ((static_cast<__int128>(v_max + 1) * (v_max + 1) + static_cast<__int128>(2) * a * (v_max + 1)) <=
           static_cast<__int128>(n)) {
        ++v_max;
    }
    while ((static_cast<__int128>(v_max) * v_max + static_cast<__int128>(2) * a * v_max) >
           static_cast<__int128>(n)) {
        --v_max;
    }

    if (v_max < static_cast<std::int64_t>(v_min)) return 0;

    factorize_int(a - 1, spf, left_factors);
    factorize_int(a + 1, spf, right_factors);
    merge_factors(left_factors, right_factors, merged_factors);

    divisors.clear();
    divisors.push_back(1);
    for (const Factor& f : merged_factors) {
        const std::size_t base_size = divisors.size();
        std::uint64_t mul = 1;
        for (int e = 1; e <= f.e; ++e) {
            mul *= static_cast<std::uint64_t>(f.p);
            for (std::size_t i = 0; i < base_size; ++i) {
                const std::uint64_t candidate = divisors[i] * mul;
                if (candidate <= static_cast<std::uint64_t>(v_max)) {
                    divisors.push_back(candidate);
                }
            }
        }
    }

    std::uint64_t count = 0;
    for (std::uint64_t v : divisors) {
        if (v < v_min) continue;

        const std::uint64_t u = n / v;
        if (((u - v) & 1ULL) != 0ULL) continue;
        if (u < v + static_cast<std::uint64_t>(2 * a)) continue;
        if (static_cast<std::uint64_t>(a) + u > static_cast<std::uint64_t>(perimeter_limit)) continue;

        ++count;
    }

    return count;
}

std::uint64_t count_barely_acute(int perimeter_limit, int threads) {
    if (perimeter_limit < 3) return 0;

    const int a_max = perimeter_limit / 3;
    const std::uint64_t a1_count = static_cast<std::uint64_t>((perimeter_limit - 1) / 2);
    const std::vector<int> spf = build_spf(a_max + 1);

    if (threads < 1) threads = 1;
    if (threads > a_max) threads = a_max;

    std::atomic<int> next_a{2};
    constexpr int kChunk = 2048;

    std::vector<std::uint64_t> partial(static_cast<std::size_t>(threads), 0);
    std::vector<std::thread> pool;
    pool.reserve(static_cast<std::size_t>(threads));

    for (int t = 0; t < threads; ++t) {
        pool.emplace_back([&, t]() {
            std::vector<Factor> left_factors;
            std::vector<Factor> right_factors;
            std::vector<Factor> merged_factors;
            std::vector<std::uint64_t> divisors;

            left_factors.reserve(8);
            right_factors.reserve(8);
            merged_factors.reserve(16);
            divisors.reserve(256);

            std::uint64_t local = 0;
            while (true) {
                const int start = next_a.fetch_add(kChunk, std::memory_order_relaxed);
                if (start > a_max) break;
                const int end = std::min(a_max, start + kChunk - 1);

                for (int a = start; a <= end; ++a) {
                    local += count_for_a(a,
                                         perimeter_limit,
                                         spf,
                                         left_factors,
                                         right_factors,
                                         merged_factors,
                                         divisors);
                }
            }

            partial[static_cast<std::size_t>(t)] = local;
        });
    }

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

    std::uint64_t total = a1_count;
    for (std::uint64_t v : partial) total += v;
    return total;
}

std::uint64_t brute_count(int perimeter_limit) {
    std::uint64_t count = 0;
    for (int a = 1; a <= perimeter_limit / 3; ++a) {
        for (int b = a; b <= (perimeter_limit - a) / 2; ++b) {
            const long long c2 = 1LL * a * a + 1LL * b * b - 1LL;
            if (c2 <= 0) continue;

            long long c = static_cast<long long>(std::sqrt(static_cast<long double>(c2)));
            while ((c + 1) * (c + 1) <= c2) ++c;
            while (c * c > c2) --c;

            if (c < b) continue;
            if (c * c != c2) continue;
            if (a + b + c > perimeter_limit) continue;
            if (a + b <= c) continue;

            ++count;
        }
    }
    return count;
}

bool validate() {
    const std::vector<int> brute_limits = {80, 120, 200, 300, 500};
    for (int limit : brute_limits) {
        const std::uint64_t slow = brute_count(limit);
        const std::uint64_t fast = count_barely_acute(limit, 1);
        if (slow != fast) {
            std::cerr << "Validation failed at perimeter=" << limit
                      << ": brute=" << slow << ", fast=" << fast << "\n";
            return false;
        }
    }

    const int thread_check_limit = 50'000;
    const std::uint64_t single = count_barely_acute(thread_check_limit, 1);

    unsigned hw = std::thread::hardware_concurrency();
    if (hw == 0) hw = 2;
    const int multi_threads = static_cast<int>(std::min<unsigned>(hw, 8));
    const std::uint64_t multi = count_barely_acute(thread_check_limit, multi_threads);
    if (single != multi) {
        std::cerr << "Thread consistency failed at perimeter=" << thread_check_limit
                  << ": single=" << single << ", multi=" << multi << "\n";
        return false;
    }

    return true;
}

}  // namespace

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

    int perimeter_limit = kDefaultLimit;
    int threads = static_cast<int>(std::thread::hardware_concurrency());
    if (threads <= 0) threads = 1;

    if (argc > 1) perimeter_limit = std::max(3, std::atoi(argv[1]));
    if (argc > 2) threads = std::max(1, std::atoi(argv[2]));

    const std::uint64_t answer = count_barely_acute(perimeter_limit, threads);
    std::cout << answer << '\n';
    return 0;
}

Python

import math
import multiprocessing

def get_spf(limit):
    spf = [0] * (limit + 1)
    primes = []
    for i in range(2, limit + 1):
        if spf[i] == 0:
            spf[i] = i
            primes.append(i)
        for p in primes:
            if p * i > limit or p > spf[i]:
                break
            spf[p * i] = p
    spf[1] = 1
    return spf

def factorize_int(x, spf):
    factors = []
    while x > 1:
        p = spf[x]
        e = 0
        while x > 1 and spf[x] == p:
            x //= p
            e += 1
        factors.append((p, e))
    return factors

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

def get_divisors(factors, v_max):
    divisors = [1]
    for p, e in factors:
        base_size = len(divisors)
        mul = 1
        for _ in range(1, e + 1):
            mul *= p
            for i in range(base_size):
                candidate = divisors[i] * mul
                if candidate <= v_max:
                    divisors.append(candidate)
    return divisors

def worker(args):
    start_a, end_a, perimeter_limit, spf = args
    local_count = 0
    
    for a in range(start_a, end_a + 1):
        n = a * a - 1
        if n == 0: continue
        
        denom = perimeter_limit - a
        if denom <= 0: continue
        
        v_min = (n + denom - 1) // denom
        
        root = math.sqrt(2.0 * a * a - 1.0)
        v_max = int(math.floor(root - a))
        if v_max <= 0: continue
        
        while (v_max + 1) * (v_max + 1) + 2 * a * (v_max + 1) <= n:
            v_max += 1
        while v_max * v_max + 2 * a * v_max > n:
            v_max -= 1
            
        if v_max < v_min: continue
        
        left_factors = factorize_int(a - 1, spf)
        right_factors = factorize_int(a + 1, spf)
        merged_factors = merge_factors(left_factors, right_factors)
        
        divisors = get_divisors(merged_factors, v_max)
        
        for v in divisors:
            if v < v_min: continue
            u = n // v
            if (u - v) % 2 != 0: continue
            if u < v + 2 * a: continue
            if a + u > perimeter_limit: continue
            local_count += 1
            
    return local_count

def solve(perimeter_limit=25000000):
    if perimeter_limit < 3: return "0"
    
    a_max = perimeter_limit // 3
    a1_count = (perimeter_limit - 1) // 2
    
    threads = multiprocessing.cpu_count() or 1
    threads = min(threads, a_max)
    
    spf = get_spf(a_max + 1)
    
    chunk = 2048
    tasks = []
    start = 2
    while start <= a_max:
        end = min(a_max, start + chunk - 1)
        tasks.append((start, end, perimeter_limit, spf))
        start += chunk
        
    total = a1_count
    
    if threads > 1 and len(tasks) > 1:
        with multiprocessing.Pool(threads) as pool:
            results = pool.map(worker, tasks)
            total += sum(results)
    else:
        for t in tasks:
            total += worker(t)
            
    return str(total)

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

Java

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

public class Euler223 {
    static class Factor {
        int p, e;

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

    static int[] buildSpf(int n) {
        int[] spf = new int[n + 1];
        List<Integer> primes = new ArrayList<>(n / 10);
        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;
            }
        }
        if (n >= 1)
            spf[1] = 1;
        return spf;
    }

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

    static void mergeFactors(List<Factor> a, List<Factor> b, List<Factor> out) {
        out.clear();
        int i = 0, j = 0;
        while (i < a.size() || j < b.size()) {
            if (j == b.size() || (i < a.size() && a.get(i).p < b.get(j).p)) {
                out.add(a.get(i++));
            } else if (i == a.size() || b.get(j).p < a.get(i).p) {
                out.add(b.get(j++));
            } else {
                out.add(new Factor(a.get(i).p, a.get(i).e + b.get(j).e));
                i++;
                j++;
            }
        }
    }

    public static String solve() {
        int perimeterLimit = 25000000;
        if (perimeterLimit < 3)
            return "0";

        int aMax = perimeterLimit / 3;
        long a1Count = (perimeterLimit - 1) / 2;
        int[] spf = buildSpf(aMax + 1);

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

        int chunk = 2048;
        List<int[]> tasks = new ArrayList<>();
        for (int start = 2; start <= aMax; start += chunk) {
            int end = Math.min(aMax, start + chunk - 1);
            tasks.add(new int[] { start, end });
        }

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

        for (int[] task : tasks) {
            final int startA = task[0];
            final int endA = task[1];

            futures.add(executor.submit(() -> {
                long local = 0;
                List<Factor> leftFactors = new ArrayList<>(8);
                List<Factor> rightFactors = new ArrayList<>(8);
                List<Factor> mergedFactors = new ArrayList<>(16);
                long[] divisors = new long[256];

                for (int a = startA; a <= endA; ++a) {
                    long n = (long) a * a - 1;
                    if (n == 0)
                        continue;

                    long denom = perimeterLimit - a;
                    if (denom <= 0)
                        continue;

                    long vMin = (n + denom - 1) / denom;
                    double root = Math.sqrt(2.0 * a * a - 1.0);
                    long vMax = (long) Math.floor(root - a);
                    if (vMax <= 0)
                        continue;

                    while (true) {
                        long nextV = vMax + 1;
                        // Avoid overflow using BigInteger or checking sizes if necessary, but max a ~
                        // 8.3M, a^2 ~ 6e13. n is 6e13.
                        // (vMax+1)^2 + 2*a*(vMax+1)
                        // vMax won't exceed sqrt(n) ~ 8.3M
                        // (8.3M)^2 ~ 6.9e13, fits in signed long (9e18).
                        long test = nextV * nextV + 2L * a * nextV;
                        if (test <= n)
                            vMax++;
                        else
                            break;
                    }
                    while (true) {
                        long test = vMax * vMax + 2L * a * vMax;
                        if (test > n)
                            vMax--;
                        else
                            break;
                    }

                    if (vMax < vMin)
                        continue;

                    factorizeInt(a - 1, spf, leftFactors);
                    factorizeInt(a + 1, spf, rightFactors);
                    mergeFactors(leftFactors, rightFactors, mergedFactors);

                    divisors[0] = 1;
                    int divSize = 1;

                    for (Factor f : mergedFactors) {
                        int baseSize = divSize;
                        long mul = 1;
                        for (int e = 1; e <= f.e; ++e) {
                            mul *= f.p;
                            for (int i = 0; i < baseSize; ++i) {
                                long cand = divisors[i] * mul;
                                if (cand <= vMax) {
                                    if (divSize == divisors.length) {
                                        divisors = Arrays.copyOf(divisors, divisors.length * 2);
                                    }
                                    divisors[divSize++] = cand;
                                }
                            }
                        }
                    }

                    for (int i = 0; i < divSize; ++i) {
                        long v = divisors[i];
                        if (v < vMin)
                            continue;

                        long u = n / v;
                        if (((u - v) & 1) != 0)
                            continue;
                        if (u < v + 2L * a)
                            continue;
                        if (a + u > perimeterLimit)
                            continue;
                        local++;
                    }
                }
                return local;
            }));
        }

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

        return String.valueOf(total);
    }

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