Problem 224: Almost Right-angled Triangles II

Project Euler Problem 224 Solution

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

Problem Summary We want the number of integer-sided triangles $(a,b,c)$ with perimeter at most $75{,}000{,}000$ such that $$a^2+b^2=c^2-1.$$ Equivalently, $c^2=a^2+b^2+1$, so the longest side is only one unit of square-length beyond the right-triangle relation. By the law of cosines, the angle opposite $c$ satisfies $$\cos C=\frac{a^2+b^2-c^2}{2ab}=-\frac{1}{2ab},$$ which explains the title: the triangle is obtuse, but only barely. Mathematical Approach The fast solution does not scan all possible side lengths. Instead it walks through a ternary tree of integer triples that preserves the quadratic form $$Q(a,b,c)=c^2-a^2-b^2.$$ Basic arithmetic structure of the equation Once $a$ and $b$ are positive, the triangle inequality is automatic. Indeed, $$c^2=a^2+b^2+1 \lt a^2+b^2+2ab=(a+b)^2,$$ so $c \lt a+b$. The code still checks $a+b \gt c$, but mathematically every positive solution already gives a genuine triangle. The parity is also forced. Squares are $0$ or $1$ modulo $4$, so $$a^2+b^2+1\equiv c^2 \pmod 4$$ rules out the cases where one of $a,b$ is odd or both are odd. Therefore every solution has $$a\equiv b\equiv 0 \pmod 2,\qquad c\equiv 1 \pmod 2.$$ The smallest positive solution is then $(2,2,3)$. The three invariant-preserving transformations Write a triple as the column vector $t=(a,b,c)^T$....

Detailed mathematical approach

Problem Summary

A triangle with integer sides $a\le b\le c$ is called barely obtuse if $a^2+b^2=c^2-1$. We must count the barely obtuse triangles whose perimeter $a+b+c$ is at most $75\,000\,000$. ^[1] The smallest one is $(2,2,3)$, since $4+4=9-1$. By the law of cosines $c^2=a^2+b^2-2ab\cos\gamma$, the angle $\gamma$ opposite $c$ has $\cos\gamma=-1/(2ab)$, so it is obtuse, but only just: about $97.2^\circ$ for $(2,2,3)$ and $90.05^\circ$ for $(18,30,35)$. ^[2] The statement gives no checkpoint values, so Section 1 computes a few.

The whole solution in one sentence: starting from $(2,2,3)$, three fixed linear rules, the same ones that build the classical tree of Pythagorean triples, generate every triple $(a,b,c)$ of positive integers with $a^2+b^2=c^2-1$ exactly once, so we follow them until the perimeter exceeds the limit and count the triples with $a\le b$.

Mathematical Approach

After small cases and simplifications we show that three linear maps generate every solution from $(2,2,3)$ exactly once; then we count the triangles and bound the size of the search.

1. Small cases

Trying all pairs $a\le b$ by computer lists the barely obtuse triangles with perimeter at most $100$:

$(2,2,3)$ with perimeter $7$ and $(4,8,9)$ with perimeter $21$;
$(12,12,17)$ with perimeter $41$ and $(6,18,19)$ with perimeter $43$;
$(8,32,33)$ with perimeter $73$ and $(18,30,35)$ with perimeter $83$.

So there are $6$ of them. Three patterns stand out. All shorter sides are even and every $c$ is odd. The triangles $(2,2,3)$, $(4,8,9)$, $(6,18,19)$ and $(8,32,33)$ have the form $(2k,2k^2,2k^2+1)$, and indeed $4k^2+4k^4+1=(2k^2+1)^2$. Two of the six, $(2,2,3)$ and $(12,12,17)$, have two equal sides. The longest sides $3,9,17,19,33,35,51,\dots$ form the OEIS sequence of the numbers $n$ for which $n^2-1$ is a sum of two nonzero squares; its example $51^2-1=10^2+50^2$${}=22^2+46^2$${}=34^2+38^2$ shows three triangles with the same longest side. ^[3]

2. What the equation already guarantees

Let positive integers $a$, $b$, $c$ satisfy $a^2+b^2+1=c^2$.

The triangle conditions are automatic. $(a+b)^2-c^2=2ab-1\gt0$, so $c\lt a+b$, and $c^2-b^2=a^2+1\gt0$, so $c\gt b$; likewise $c\gt a$. Of the conditions $a\le b\le c$ and the triangle inequality, only $a\le b$ remains.
Parity. A square leaves the remainder $0$ or $1$ when divided by $4$, because $(2k)^2=4k^2$ and $(2k+1)^2=4k^2+4k+1$. If exactly one of $a,b$ were odd, $a^2+b^2+1$ would leave the remainder $2$; if both were odd, $3$. Neither is possible for $c^2$, so $a$ and $b$ are even, $c$ is odd, and $a,b\ge2$. ^[3]

From now on a solution is an ordered triple $(a,b,c)$ of positive integers with $a^2+b^2+1=c^2$, with no condition on the order of $a$ and $b$. A triangle with $a\lt b$ gives the two solutions $(a,b,c)$ and $(b,a,c)$; a triangle with $a=b$ gives one.

3. Three maps that preserve the equation

Put $Q(a,b,c)=c^2-a^2-b^2$. The solutions are the positive triples with $Q=1$; positive triples with $Q=0$ are the Pythagorean triples, called primitive when their entries have no common divisor above $1$. Let $M_2$ send $(a,b,c)$ to $(a',b',c')$ with

$$\begin{aligned}a'&=a+2b+2c,\\ b'&=2a+b+2c,\\ c'&=2a+2b+3c,\end{aligned}$$

and let $M_1(a,b,c)=M_2(a,-b,c)$ and $M_3(a,b,c)=M_2(-a,b,c)$. Written as $3\times3$ matrices, $M_1$, $M_2$ and $M_3$ are the matrices $A$, $B$ and $C$ of Berggren's tree of primitive Pythagorean triples. ^[4] Expanding the squares gives

$$\begin{aligned}&(2a+2b+3c)^2-(a+2b+2c)^2\\ &\quad-(2a+b+2c)^2=c^2-a^2-b^2,\end{aligned}$$

because the coefficients of $a^2$, $b^2$ and $c^2$ are $4-1-4=-1$, $4-4-1=-1$ and $9-4-4=1$, and those of $ab$, $ac$ and $bc$ are $8-4-4$, $12-4-8$ and $12-8-4$, all zero. So $Q(M_2t)=Q(t)$ for every triple $t$, and since changing the sign of $a$ or $b$ does not change $Q$, the maps $M_1$ and $M_3$ preserve $Q$ as well.

4. Children are larger

Call $M_1t$, $M_2t$ and $M_3t$ the children of a solution $t=(a,b,c)$. Subtracting $t$ from them and using $c-a\ge1$ and $c-b\ge1$ from Section 2:

$M_1$ adds $2(c-b)$ to $a$, and $2a+2(c-b)$ to $b$ and to $c$;
$M_2$ adds $2b+2c$, $2a+2c$ and $2a+2b+2c$ to $a$, $b$ and $c$;
$M_3$ adds $2b+2(c-a)$ to $a$ and to $c$, and $2(c-a)$ to $b$.

Every entry grows by at least $2$. Hence the children of a solution are again solutions, and along every branch of the tree all three sides and the perimeter strictly increase: once a triple exceeds the perimeter limit, so does everything below it, which is why the search may turn back there.

5. Every solution has exactly one parent

Now we go upwards. The inverse of $M_2$ sends $(a,b,c)$ to $(s_1,s_2,s_3)$ with

$$\begin{aligned}s_1&=a+2b-2c,\\ s_2&=2a+b-2c,\\ s_3&=3c-2a-2b,\end{aligned}$$

as substituting into the formulas of Section 3 confirms. Let $t=(a,b,c)$ be a solution and $s=M_2^{-1}t$. Then $Q(s)=1$, that is $s_3^2=s_1^2+s_2^2+1$, and:

$0\lt s_3\lt c$. The right inequality is $c\lt a+b$ from Section 2. The left one means $9c^2\gt4(a+b)^2$, that is $5a^2-8ab+5b^2+9\gt0$, which holds because $5a^2-8ab+5b^2$${}=4(a-b)^2+a^2+b^2$.
If $a\le b$, then $s_1\ge0$: indeed $(a+2b)^2-4c^2$${}=4a(b-a)+a^2-4\ge0$ because $a\ge2$, with equality only for $a=b=2$. In the same way $b\le a$ gives $s_2\ge0$. So $s_1$ and $s_2$ are never both negative.
If $s_1=0$, then $(s_3-s_2)(s_3+s_2)=1$ with $s_3\gt0$, so both factors equal $1$, $s_2=0$, $s_3=1$ and $t=M_2(0,0,1)=(2,2,3)$; likewise if $s_2=0$.

Consequently every solution $t\ne(2,2,3)$ falls into exactly one of three cases: if $s_1,s_2\gt0$ then $t=M_2s$; if $s_1\gt0\gt s_2$ then $t=M_1(s_1,-s_2,s_3)$; if $s_1\lt0\lt s_2$ then $t=M_3(-s_1,s_2,s_3)$. In each case $t$ is a child of a solution, its parent, whose third entry $s_3$ is smaller than $c$. No other solution $u$ has $t$ as a child: $t=M_iu$ gives $M_2^{-1}t=(u_1,-u_2,u_3)$, $u$ or $(-u_1,u_2,u_3)$ for $i=1,2,3$, so the signs of $s_1$ and $s_2$ determine both $i$ and $u$.

Going from parent to parent, the third entry decreases at every step, so after finitely many steps we arrive at $(2,2,3)$, the only solution without a parent: its image $M_2^{-1}(2,2,3)=(0,0,1)$ is not a solution, and all three maps send $(0,0,1)$ to $(2,2,3)$. With Section 4 and an induction on $c$ this proves the key claim: the tree with root $(2,2,3)$ and the three branches $M_1$, $M_2$, $M_3$ at every node contains every solution exactly once. An analogous descent shows that Berggren's tree with root $(3,4,5)$ contains every primitive Pythagorean triple exactly once, written with its odd leg first; ^[4] in our tree both legs are even and both orders occur. In 1963 Barning generated the solutions of $a^2+b^2=c^2-1$ with the same matrices from the root $(2,2,3)$. ^[5]

6. Worked example

Take $(18,30,35)$. Its image under $M_2^{-1}$ is $(8,-4,9)$, so $s_2\lt0$ and $(18,30,35)=M_1(8,4,9)$. Next, $M_2^{-1}(8,4,9)=(-2,2,3)$, so $(8,4,9)=M_3(2,2,3)$. The path from the root is therefore $M_3$ followed by $M_1$, that is $(18,30,35)=M_1M_3(2,2,3)$.

Going down instead, this is the whole tree up to perimeter $100$:

depth $0$: the root $(2,2,3)$, perimeter $7$;
depth $1$: $M_1$, $M_2$, $M_3$ give $(4,8,9)$, $(12,12,17)$, $(8,4,9)$, perimeters $21$, $41$, $21$;
depth $2$: $M_1$ and $M_3$ turn $(4,8,9)$ into $(6,18,19)$ and $(30,18,35)$, and $(8,4,9)$ into $(18,30,35)$ and $(18,6,19)$, perimeters $43$, $83$, $83$, $43$;
depth $3$: $M_1(6,18,19)=(8,32,33)$ and $M_3(18,6,19)=(32,8,33)$, perimeter $73$ each.

All other children are too large; for example $M_2(4,8,9)=(38,34,51)$ has perimeter $123$, and the children of $(12,12,17)$ have perimeters $119$, $239$ and $119$. These are $10$ solutions, and the $6$ with $a\le b$ are exactly the triangles of Section 1. The example also shows why the walk must not stop at a solution with $a\gt b$: $(8,4,9)$ is not counted, but its child $(18,30,35)$ is. In general a child under $M_1$ has $a'\lt b'$, because $b'-a'=a+b$, a child under $M_3$ has $a'\gt b'$, and $M_2$ reverses the order, because $a'-b'=b-a$.

7. Mirror images and equal sides

Let $S$ exchange the first two entries of a triple. Exchanging $a$ and $b$ in the input and in the output of $M_1$ gives $M_3$, while $M_2$ stays the same; that is, $SM_1S=M_3$ and $SM_2S=M_2$. Hence the mirror image of the solution at the end of a path is the solution at the end of the path with every $M_1$ replaced by $M_3$ and every $M_3$ by $M_1$: the mirror image of $(18,30,35)=M_1M_3(2,2,3)$ is $(30,18,35)=M_3M_1(2,2,3)$.

A solution with $a=b$ is its own mirror image, so by the uniqueness of paths its path contains neither $M_1$ nor $M_3$: the solutions with $a=b$ are exactly the triples $M_2^k(2,2,3)$. On a triple with $a=b$, $M_2$ sends $a$ to $3a+2c$ and $c$ to $4a+3c$, which multiplies $c+a\sqrt2$ by $3+2\sqrt2$. They are the solutions of Pell's equation $c^2-2a^2=1$, which all come from powers of the fundamental solution $3+2\sqrt2$: ^[6] $(2,2,3)$, $(12,12,17)$, $(70,70,99)$, $(408,408,577)$, and so on.

If $N(P)$ solutions have perimeter at most $P$ and $D(P)$ of them have $a=b$, the number of barely obtuse triangles is

$$T(P)=\frac{N(P)+D(P)}{2},$$

because every triangle with $a\lt b$ appears as two solutions. For $P=100$ this gives $(10+2)/2=6$. The programs count the solutions with $a\le b$ instead, which gives the same number.

8. How deep the tree goes

By Section 4, a step from a solution raises $\min(a,b)$ by at least $2$ and raises $c$ by at least $2\min(a,b)+2$. A solution at depth $d$, that is $d$ steps below the root, therefore has $\min(a,b)\ge2d+2$, and the step to depth $j$ adds at least $4j+2$ to $c$, so

$$c\ge3+\sum_{j=1}^{d}(4j+2)=2(d+1)^2+1.$$

Equality holds on the path that uses only $M_1$, which leaves $c-b=1$ unchanged: $M_1^d(2,2,3)=(2k,2k^2,2k^2+1)$ with $k=d+1$ is the family of Section 1. Since the perimeter exceeds $2c$, a solution at depth $d$ with perimeter at most $P$ has $4(d+1)^2+2\lt P$, so $d+1\lt\sqrt P/2$. For $P=75\,000\,000$ this gives $d\le4329$; the deepest solutions within the limit actually have depth $4328$ and lie on the paths that use only $M_1$ or only $M_3$.

9. A count without the tree

The equation can also be written as $a^2+1=c^2-b^2=(c-b)(c+b)$. For even $a$, every factorisation $a^2+1=uv$ with $u\lt v$ has odd factors and gives $b=(v-u)/2$ and $c=(v+u)/2$, with perimeter $a+v$, and every solution arises in this way from $u=c-b$ and $v=c+b$. For example $12^2+1=145=1\cdot145=5\cdot29$ gives $(12,72,73)$ and $(12,12,17)$, and $u=1$ always gives the family $(2k,2k^2,2k^2+1)$. Counting the factorisations with $b\ge a$ and $a+v\le P$ requires the factorisation of $a^2+1$ for every even $a\le P/3$. This is slower than the tree but independent of it, which makes it a good check.

How the Code Works

All three programs walk through the tree of Section 5 depth first, keeping the triples still to be visited on an explicit stack instead of using recursion. ^[7] In the C++ program ^[8] a Triple holds $a$, $b$, $c$ as 64-bit integers, the constant matrices kM1, kM2, kM3 are $M_1$, $M_2$, $M_3$, and transform multiplies a matrix by a triple. count_barely_obtuse pushes $(2,2,3)$ onto a std::vector used as the stack, with room for $2^{20}$ entries reserved in advance, and repeats until the stack is empty:

pop the top triple and skip it if its perimeter exceeds the limit;
add $1$ to the counter if $a\le b$ and $a+b\gt c$;
form the three children, drop any child with an entry $\le0$, and push each child whose perimeter is within the limit.

Three of these tests never reject anything, except the root when the limit is below $7$: a popped child was checked before it was pushed, $a+b\gt c$ holds by Section 2, and children are positive by Section 4. The test $a\le b$ picks one solution from each mirror pair (Section 7), while the walk continues below every node within the limit, as Section 6 requires.

Before the main count, run_checkpoints compares the tree count with brute_count for the perimeters $80$, $120$, $200$, $500$, $1000$ and $2000$. The brute force tries every $a\le P/3$ and $a\le b\le(P-a)/2$, takes $c$ as the integer square root of $a^2+b^2+1$ (a long double square root followed by integer corrections), and accepts the triple if $c^2=a^2+b^2+1$, $a+b\gt c$ and $a+b+c\le P$. Both methods give $5$, $8$, $12$, $29$, $55$ and $114$ triangles; a mismatch would end the program with exit status 2. The options --skip-checkpoints and --perimeter=N skip this step or change the limit.

The Python version does the same walk in solve, with the matrices as nested lists, a local transform and a list of tuples as the stack; it has no checkpoints and no options. The Java version keeps the stack in three long arrays stackA, stackB, stackC, one per coordinate, of length $10^6$ with an index top, and doubles them with Arrays.copyOf when they are full; it has no checkpoints either. At the given limit nothing can overflow: a child of a triple with perimeter at most $P$ has entries at most $2a+2b+3c\le3P$ and perimeter at most $5a+5b+7c\le7P=5.25\cdot10^8$, so even 32-bit integers would suffice; C++ and Java use 64-bit integers and Python unbounded ones. The C++ option parser does not check whether a --perimeter value fits in 64 bits.

Complexity and Verification

Every solution with perimeter at most $P$ is pushed and popped exactly once (Sections 4 and 5), and each one costs three products of a $3\times3$ matrix with a vector, so the running time is proportional to $N(P)$, the number of solutions itself. We do not prove a growth rate for $N(P)$, but the measured values $N(10^4)=1081$, $N(10^5)=11018$ and $N(10^6)=110239$ grow roughly in proportion to $P$. By contrast, the inner loop of brute_count tries about $(P-3a)/2$ values of $b$ for each $a\le P/3$, about $P^2/12\approx4.7\cdot10^{14}$ pairs in total for the given limit. On a current desktop computer the C++ and Java programs finish in a fraction of a second and the Python program in several seconds.

The stack holds at most two waiting siblings for each level of the current path plus the three newest children, so by Section 8 it never has more than $2d+3\lt\sqrt P+1$ entries, where $d$ is the largest depth; for the given limit the largest stack actually reached has $3535$ entries. The storage set aside is larger but fixed: about 25 MB for the $2^{20}$ triples reserved in C++ and 24 MB for the three Java arrays of $10^6$ long values, which never need to grow; the Python list grows with the stack.

The C++ checkpoints above are the program's own test. For this article we also checked that:

the tree and the factorisation method of Section 9 give exactly the same set of $1652$ triangles with perimeter at most $30\,000$, and hence agree for every smaller limit too;
every solution up to perimeter $30\,000$ has the parent, the path and the mirror image predicted by Sections 5 and 7, and every solution up to the full limit satisfies the bounds of Sections 4 and 8;
an independent C++ program that counts the factorisations of Section 9 with a sieve over $a^2+1$ finds $543$, $5512$ and $55123$ triangles for $P=10^4$, $10^5$ and $10^6$, as the tree does, and agrees with the tree at the full limit, where the three published programs also agree with one another.

Footnotes and References

The references below supply the background facts used above. The parent rule, the depth bound and the operation counts are derived explicitly in this article.

Project Euler 224 — Almost Right-angled Triangles II. The official statement calls an integer-sided triangle with sides a ≤ b ≤ c barely obtuse if a² + b² = c² − 1 and asks how many barely obtuse triangles have a perimeter of at most 75 000 000.
Wikipedia — Law of cosines. The law c² = a² + b² − 2ab cos γ, with γ the angle opposite c, generalises the Pythagorean theorem; an obtuse angle has a negative cosine.
OEIS A050795. Numbers n such that n² − 1 is a sum of two nonzero squares, that is, the longest sides of the barely obtuse triangles: 3, 9, 17, 19, 33, 35, 51, … A comment notes that a² + b² = c² − 1 modulo 4 forces a and b to be even and c to be odd, and the example gives the three representations of 51² − 1.
Wikipedia — Tree of primitive Pythagorean triples. The article gives Berggren's matrices A, B, C (here M₁, M₂, M₃) and proves with their inverses that every primitive triple, written with its odd leg first, is reached from (3, 4, 5) along exactly one path: exactly one inverse gives a triple with positive entries and a smaller hypotenuse.
F. J. M. Barning, Over Pythagorese en bijna-Pythagorese driehoeken en een generatieproces met behulp van unimodulaire matrices, Mathematisch Centrum Amsterdam, report ZW 1963-011 (in Dutch). Section 4 generates the triangles with x² + y² = z² ± 1 with the matrices of the Pythagorean tree, which leave x² + y² − z² invariant. For z² − 1 the tree starts at (2, 2, 3); completeness follows from the same descent as for Pythagorean triangles, mirror copies are avoided by not applying one matrix to triangles with equal legs, and these lie on the straight path of the matrix (x + 2y + 2z, 2x + y + 2z, 2x + 2y + 3z).
Wikipedia — Pell's equation. For x² − ny² = 1 all solutions follow from the fundamental solution by x_k + y_k√n = (x₁ + y₁√n)^k; for n = 2 the fundamental solution is (3, 2), and (17, 12) and (577, 408) are quoted as early examples.
Wikipedia — Depth-first search. The article describes the non-recursive form of depth-first search, which keeps the vertices still to be explored on an explicit stack with push and pop operations.
C++ — Euler224.cpp. The linked immutable C++ revision contains count_barely_obtuse, brute_count and run_checkpoints with the perimeters 80, 120, 200, 500, 1000 and 2000. It is the implementation record for this article.

Problem 224 source code

C++

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

namespace {

using i64 = long long;
using u64 = std::uint64_t;

struct Options {
    i64 perimeter_limit = 75'000'000;
    bool run_checkpoints = true;
};

struct Triple {
    i64 a;
    i64 b;
    i64 c;
};

bool parse_int_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_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--perimeter=", options.perimeter_limit)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.perimeter_limit >= 3;
}

Triple transform(const std::array<std::array<i64, 3>, 3>& m, const Triple& t) {
    const i64 x = m[0][0] * t.a + m[0][1] * t.b + m[0][2] * t.c;
    const i64 y = m[1][0] * t.a + m[1][1] * t.b + m[1][2] * t.c;
    const i64 z = m[2][0] * t.a + m[2][1] * t.b + m[2][2] * t.c;
    return {x, y, z};
}

u64 count_barely_obtuse(const i64 perimeter_limit) {
    static constexpr std::array<std::array<i64, 3>, 3> kM1{{
        {{1, -2, 2}},
        {{2, -1, 2}},
        {{2, -2, 3}},
    }};
    static constexpr std::array<std::array<i64, 3>, 3> kM2{{
        {{1, 2, 2}},
        {{2, 1, 2}},
        {{2, 2, 3}},
    }};
    static constexpr std::array<std::array<i64, 3>, 3> kM3{{
        {{-1, 2, 2}},
        {{-2, 1, 2}},
        {{-2, 2, 3}},
    }};

    std::vector<Triple> stack;
    stack.reserve(1 << 20);
    stack.push_back({2, 2, 3});

    u64 count = 0;

    while (!stack.empty()) {
        const Triple cur = stack.back();
        stack.pop_back();

        const i64 perimeter = cur.a + cur.b + cur.c;
        if (perimeter > perimeter_limit) {
            continue;
        }

        if (cur.a <= cur.b && cur.a + cur.b > cur.c) {
            ++count;
        }

        const Triple children[3] = {
            transform(kM1, cur),
            transform(kM2, cur),
            transform(kM3, cur),
        };

        for (const Triple& child : children) {
            if (child.a <= 0 || child.b <= 0 || child.c <= 0) {
                continue;
            }
            if (child.a + child.b + child.c <= perimeter_limit) {
                stack.push_back(child);
            }
        }
    }

    return count;
}

u64 brute_count(const int perimeter_limit) {
    u64 count = 0;
    for (int a = 1; a <= perimeter_limit / 3; ++a) {
        for (int b = a; b <= (perimeter_limit - a) / 2; ++b) {
            const i64 c2 = 1LL * a * a + 1LL * b * b + 1;
            i64 c = static_cast<i64>(std::sqrt(static_cast<long double>(c2)));
            while ((c + 1) * (c + 1) <= c2) {
                ++c;
            }
            while (c * c > c2) {
                --c;
            }
            if (c * c != c2) {
                continue;
            }
            if (a + b <= c) {
                continue;
            }
            if (a + b + c > perimeter_limit) {
                continue;
            }
            ++count;
        }
    }
    return count;
}

bool run_checkpoints() {
    for (int perimeter : {80, 120, 200, 500, 1000, 2000}) {
        const u64 slow = brute_count(perimeter);
        const u64 fast = count_barely_obtuse(perimeter);
        if (slow != fast) {
            std::cerr << "Checkpoint failed at perimeter=" << perimeter
                      << ": brute=" << slow << ", fast=" << fast << '\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 << count_barely_obtuse(options.perimeter_limit) << '\n';
    return 0;
}

Python

def solve():
    perimeter_limit = 75_000_000

    M1 = [[1, -2, 2], [2, -1, 2], [2, -2, 3]]
    M2 = [[1, 2, 2], [2, 1, 2], [2, 2, 3]]
    M3 = [[-1, 2, 2], [-2, 1, 2], [-2, 2, 3]]

    def transform(m, t):
        return (
            m[0][0]*t[0] + m[0][1]*t[1] + m[0][2]*t[2],
            m[1][0]*t[0] + m[1][1]*t[1] + m[1][2]*t[2],
            m[2][0]*t[0] + m[2][1]*t[1] + m[2][2]*t[2]
        )

    stack = [(2, 2, 3)]
    count = 0

    while stack:
        a, b, c = stack.pop()
        perimeter = a + b + c
        if perimeter > perimeter_limit:
            continue
        if a <= b and a + b > c:
            count += 1
        for m in (M1, M2, M3):
            child = transform(m, (a, b, c))
            if child[0] > 0 and child[1] > 0 and child[2] > 0:
                if child[0] + child[1] + child[2] <= perimeter_limit:
                    stack.append(child)

    return str(count)

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

Java

import java.util.*;

public class Euler224 {
    static class Triple {
        long a, b, c;

        Triple(long a, long b, long c) {
            this.a = a;
            this.b = b;
            this.c = c;
        }
    }

    public static String solve() {
        long perimeterLimit = 75000000;

        long[][] m1 = {
                { 1, -2, 2 },
                { 2, -1, 2 },
                { 2, -2, 3 }
        };
        long[][] m2 = {
                { 1, 2, 2 },
                { 2, 1, 2 },
                { 2, 2, 3 }
        };
        long[][] m3 = {
                { -1, 2, 2 },
                { -2, 1, 2 },
                { -2, 2, 3 }
        };

        long[][][] matrices = { m1, m2, m3 };

        long count = 0;

        // Custom stack for performance
        long[] stackA = new long[1000000];
        long[] stackB = new long[1000000];
        long[] stackC = new long[1000000];
        int top = 0;

        stackA[0] = 2;
        stackB[0] = 2;
        stackC[0] = 3;
        top++;

        while (top > 0) {
            top--;
            long a = stackA[top];
            long b = stackB[top];
            long c = stackC[top];

            if (a + b + c > perimeterLimit) {
                continue;
            }

            if (a <= b && a + b > c) {
                count++;
            }

            for (int i = 0; i < 3; i++) {
                long[][] m = matrices[i];
                long childA = m[0][0] * a + m[0][1] * b + m[0][2] * c;
                long childB = m[1][0] * a + m[1][1] * b + m[1][2] * c;
                long childC = m[2][0] * a + m[2][1] * b + m[2][2] * c;

                if (childA <= 0 || childB <= 0 || childC <= 0)
                    continue;
                if (childA + childB + childC <= perimeterLimit) {
                    if (top >= stackA.length) {
                        stackA = Arrays.copyOf(stackA, stackA.length * 2);
                        stackB = Arrays.copyOf(stackB, stackB.length * 2);
                        stackC = Arrays.copyOf(stackC, stackC.length * 2);
                    }
                    stackA[top] = childA;
                    stackB[top] = childB;
                    stackC[top] = childC;
                    top++;
                }
            }
        }

        return String.valueOf(count);
    }

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