Problem 210: Obtuse Angled Triangles

Project Euler Problem 210 Solution

EulerSolve provides an optimized solution for Project Euler Problem 210, Obtuse Angled Triangles, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $r$ be divisible by 4 and set $a=r/4$. The fixed vertices are $O=(0,0)$ and $C=(a,a)$. The third vertex $B=(x,y)$ ranges over all lattice points in the taxicab diamond $$|x|+|y|\le r.$$ The task is to count how many non-degenerate triangles $OBC$ are obtuse. The geometric difficulty is that the allowed points form an $L_1$-ball rather than a Euclidean circle, while the obtuse-angle test itself is Euclidean. The implementations resolve this by splitting the count according to which vertex is obtuse and then converting the only nonlinear case into a lattice-point count inside a circle with a parity restriction. Mathematical Approach Separating the three obtuse-angle cases A triangle is obtuse exactly when one of its angles has negative dot product. With $B=(x,y)$ and $C=(a,a)$, the three angle tests are $$\overrightarrow{OB}\cdot\overrightarrow{OC}=a(x+y),$$ $$\overrightarrow{CO}\cdot\overrightarrow{CB}=2a^2-a(x+y),$$ $$\overrightarrow{BO}\cdot\overrightarrow{BC}=x^2+y^2-a(x+y).$$ Therefore the triangle is obtuse at $$O \iff x+y \lt 0,$$ $$C \iff x+y \gt 2a=\frac r2,$$ $$B \iff x^2+y^2-a(x+y) \lt 0.$$ These three regions can be counted separately, because a non-degenerate triangle cannot have more than one obtuse angle. The only overlap that matters is the collinear line $y=x$, which will be removed at the end....

Detailed mathematical approach

Problem Summary

For a positive integer $r$, the set $S(r)$ consists of the points $(x,y)$ with integer coordinates and $|x|+|y|\le r$: the lattice points of a square standing on one corner, the closed disk of radius $r$ in taxicab geometry. ^[2] Let $O=(0,0)$ and $C=(r/4,r/4)$. Then $N(r)$ is the number of points $B$ of $S(r)$ for which the triangle $OBC$ has an obtuse angle, that is, its largest angle $\alpha$ satisfies $90^\circ\lt\alpha\lt180^\circ$. The upper bound excludes the points $B$ of the line through $O$ and $C$: they give flat triangles with a $180^\circ$ angle, or no triangle at all for $B=O$ and $B=C$. The statement gives $N(4)=24$ and $N(8)=100$ and asks for $N(1\,000\,000\,000)$. ^[1] All radii in this article are multiples of $4$, as $10^9$ is, so $C=(a,a)$ with the integer $a=r/4$ is a lattice point.

The whole solution in one sentence: apart from the points of the line $OC$, which are subtracted at the end, the angle at $O$ or at $C$ is obtuse exactly when $B$ lies beyond one of the lines $x+y=0$ and $x+y=r/2$, which gives two short formulas, and the angle at $B$ is obtuse exactly when $B$ lies strictly inside the circle with diameter $OC$, so $N(r)$ is a polynomial in $r$ plus the number of lattice points in one disk, which we count column by column in a number of steps proportional to $r$.

Mathematical Approach

We classify a tiny diamond by hand, then turn every angle condition into an inequality for the coordinates of $B$.

1. A small case by hand

For $r=4$ we have $a=1$ and $C=(1,1)$. The diamond $S(4)$ has $41$ points. Five of them, $(-2,-2)$, $(-1,-1)$, $O$, $C$ and $(2,2)$, lie on the line $y=x$ through $O$ and $C$ and give no proper triangle. The other $36$ points give $24$ obtuse, $10$ right-angled and $2$ acute triangles, which confirms $N(4)=24$. For instance, $B=(-1,0)$ lies behind $O$ as seen from $C$, and the angle at $O$ is obtuse; $B=(1,0)$ gives a right angle at $B$; and $B=(-1,2)$ gives an acute triangle. For a billion we need an exact test.

2. The dot-product test

For two vectors, $\mathbf p\cdot\mathbf q=p_1q_1+p_2q_2$, and also $\mathbf p\cdot\mathbf q=|\mathbf p|\,|\mathbf q|\cos\theta$, where $\theta$ is the angle between them. ^[4] The cosine is negative exactly for $90^\circ\lt\theta\le180^\circ$. So the angle of a triangle at a vertex is obtuse exactly when the dot product of the two vectors from that vertex to the other two vertices is negative; it is right when the product is $0$ and acute when the product is positive. With $B=(x,y)$ and $C=(a,a)$ the three products are

$$\begin{aligned}\overrightarrow{OB}\cdot\overrightarrow{OC}&=a(x+y),\\\overrightarrow{CO}\cdot\overrightarrow{CB}&=a(2a-x-y),\\\overrightarrow{BO}\cdot\overrightarrow{BC}&=x^2+y^2-a(x+y).\end{aligned}$$

So for $B$ off the line $y=x$, where $OBC$ is a proper triangle, the angle at $O$ is obtuse exactly when $x+y\lt0$, because $a\gt0$; the angle at $C$ exactly when $x+y\gt2a=r/2$; and the angle at $B$ exactly when $x^2+y^2\lt a(x+y)$. For the three points of Section 1, where $a=1$, the products are $-1$, $3$, $2$ for $B=(-1,0)$; $1$, $1$, $0$ for $B=(1,0)$; and $1$, $1$, $4$ for $B=(-1,2)$.

3. Flat triangles and three separate counts

The angles of a proper triangle add up to $180^\circ$, so at most one of them is obtuse. Hence for every $B$ off the line $y=x$ at most one of the three products is negative, and $B$ is counted exactly when one of them is. On the line, a point $B=(t,t)$ with $t\ne0$ and $t\ne a$ is collinear with $O$ and $C$, and the angles are $0^\circ$, $0^\circ$ and $180^\circ$: the product at the middle point is negative, because $\cos180^\circ=-1$, and the other two are positive. So one test fires there although these points must not be counted. For $B=O$ or $B=C$ one of the vectors is zero, and no product is negative.

Let $N_O$, $N_C$ and $N_B$ be the numbers of points of $S(r)$ with a negative product at $O$, at $C$ and at $B$. By the argument above, no point belongs to two of these sets. The line $y=x$ meets $S(r)$ in the $r+1$ points with $|t|\le r/2$, and apart from $O$ and $C$ each of them lies in exactly one set: the $r/2$ points with $t\lt0$ in the first, the $a$ points with $a\lt t\le2a$ in the second and the $a-1$ points with $0\lt t\lt a$ in the third. These are $r-1$ points together, and therefore

$$N(r)=N_O+N_C+N_B-(r-1).$$

For $r=4$ the three subtracted points are $(-1,-1)$ and $(-2,-2)$, counted in $N_O$, and $(2,2)$, counted in $N_C$.

4. The size of the diamond

The row $y=j$ of $S(r)$, for $|j|\le r$, contains the $2(r-|j|)+1$ points with $|x|\le r-|j|$. Adding over all rows gives

$$|S(r)|=(2r+1)+2\sum_{j=-r}^{r}(r-|j|)=2r^2+2r+1,$$

because the sum equals $r+2\,(0+1+\dots+(r-1))=r^2$. These are the centred square numbers $r^2+(r+1)^2$, which the OEIS lists with exactly this lattice-point meaning. ^[3] For $r=4$ we get the $41$ points of Section 1, and for $r=8$ we get $145$.

5. The angle at O: half of the diamond

The point reflection $(x,y)\mapsto(-x,-y)$ maps $S(r)$ onto itself and changes the sign of $x+y$, so it pairs the points with $x+y\lt0$ with those with $x+y\gt0$. The remaining points lie on the line $x+y=0$; they are $(t,-t)$ with $|t|\le r/2$, so there are $r+1$ of them. Hence

$$N_O=\frac{(2r^2+2r+1)-(r+1)}{2}=r^2+\frac r2.$$

For $r=4$ this is $(41-5)/2=18$, and for $r=8$ it is $68$.

6. The angle at C: a band of diagonals

For the condition $x+y\gt r/2$ we turn the diamond by $45^\circ$, using $p=x+y$ and $q=x-y$. If $x$ and $y$ have the same sign, then $|x+y|=|x|+|y|\ge|x-y|$; otherwise the roles of $x+y$ and $x-y$ are exchanged. So $|x|+|y|=\max(|p|,|q|)$, and $S(r)$ becomes the square $|p|\le r$, $|q|\le r$. Since $x=(p+q)/2$ and $y=(p-q)/2$, the integer points correspond exactly to the integer pairs with $p\equiv q\pmod2$. On the diagonal $p=s$ the admissible $q$ are the integers of $[-r,r]$ with the parity of $s$: $r+1$ of them for even $s$ and $r$ for odd $s$. There are $r/2$ diagonals $s=r/2+1,\dots,r$, and since $r/2$ is even they start with an odd $s$ and alternate, so $r/4$ of them are odd and $r/4$ are even:

$$N_C=\frac r4\,(r+1)+\frac r4\,r=\frac{r(2r+1)}4.$$

For $r=8$ the diagonals $s=5,6,7,8$ hold $8+9+8+9=34$ points.

7. The angle at B: the circle of Thales

Let $M=(a/2,a/2)$ be the midpoint of $OC$, and let $\mathbf h=M-O=C-M$. Then $O-B=(M-B)-\mathbf h$ and $C-B=(M-B)+\mathbf h$, and multiplying out gives

$$\overrightarrow{BO}\cdot\overrightarrow{BC}=|BM|^2-|\mathbf h|^2,\qquad|\mathbf h|^2=\frac{a^2}2.$$

So, for $B$ off the line $y=x$, the angle at $B$ is obtuse when $B$ lies strictly inside the circle with centre $M$ and diameter $OC$, right when $B$ lies on it and acute when $B$ lies outside, by Thales's theorem and its companion statements. ^[5] In coordinates, $x^2+y^2-a(x+y)$${}=(x-a/2)^2+(y-a/2)^2-a^2/2$. For $r=4$ the circle passes through $O$, $(1,0)$, $C$ and $(0,1)$, which explains the right angle at $B=(1,0)$ in Section 1 and the one at $(0,1)$, and no lattice point lies strictly inside it, so $N_B=0$.

The disk needs no help from the diamond. For every point inside it, $|x|+|y|\le|x-a/2|+|y-a/2|+a$${}\le\sqrt2\,|BM|+a\lt2a=r/2$, where the middle step uses $(|s|+|t|)^2\le2(s^2+t^2)$, which is $(|s|-|t|)^2\ge0$ rearranged. So $N_B$ is simply the number of lattice points strictly inside the circle, including the $a-1$ points of the line $y=x$ that Section 3 subtracts.

8. Whole numbers for the circle

Multiplying the condition $(x-a/2)^2+(y-a/2)^2\lt a^2/2$ by $4$ and writing $u=2x-a$ and $v=2y-a$ gives $u^2+v^2\lt2a^2$. The new coordinates are integers with the same parity as $a$, and conversely every pair of integers with $u\equiv v\equiv a\pmod2$ comes from the lattice point $x=(u+a)/2$, $y=(v+a)/2$. Because $u^2+v^2$ is an integer, the strict inequality is the same as $u^2+v^2\le L$ with $L=2a^2-1$. So $N_B$ is the number of integer pairs $(u,v)$ with $u\equiv v\equiv a\pmod2$ and $u^2+v^2\le L$. The disk is now centred at the origin, so it is symmetric under $u\mapsto-u$ and $v\mapsto-v$.

9. How many points the disk holds

Back in the coordinates $(x,y)$, the disk has centre $M$, radius $\rho=a/\sqrt2$ and area $\pi a^2/2$, and $N_B$ is close to this area, as in the Gauss circle problem: on average every unit square contains one lattice point. ^[7] Precisely, the unit squares centred at the lattice points inside the disk have disjoint interiors and lie within distance $\rho+\sqrt2/2$ of $M$. They also cover every point at distance less than $\rho-\sqrt2/2$ from $M$, because such a point lies in the unit square of a lattice point at distance at most $\sqrt2/2$, which is then inside the disk. Comparing areas gives

$$\frac{\pi(a-1)^2}2\le N_B\le\frac{\pi(a+1)^2}2.$$

For $r=4000$, that is $a=1000$, the count is $N_B=1570777$, against the area $1570796.3$. For $r=10^9$ the disk holds roughly $10^{17}$ lattice points, far too many to visit one by one; the next three sections count them in about $3.5\cdot10^8$ steps.

10. Counting the disk column by column

Fix a column $u\equiv a\pmod2$ with $u^2\le L$. If some integer $v$ with the parity of $a$ satisfies $v^2\le L-u^2$, let $V(u)$ be the largest one, which is nonnegative, and the admissible values of $v$ in this column are $-V(u),-V(u)+2,\dots,V(u)$, which are $V(u)+1$ numbers. Otherwise the column is empty, and we set $V(u)=-1$, so that $V(u)+1=0$ still counts its points. The columns $u$ and $-u$ hold the same number of points, so

$$N_B=\sum_{u}m(u)\,\bigl(V(u)+1\bigr),$$

where $u$ runs over the nonnegative columns, and $m(0)=1$ while $m(u)=2$ for $u\gt0$. For $r=8$ we have $a=2$ and $L=7$. The column $u=0$ allows $v^2\le7$, so $V(0)=2$ and $v$ takes the values $-2$, $0$, $2$; the column $u=2$ allows $v^2\le3$, so $V(2)=0$. Thus $N_B=3+2\cdot1=5$. Going back with $x=(u+2)/2$ and $y=(v+2)/2$, these are the points $(1,0)$, $(0,1)$, $(1,1)$, $(2,1)$ and $(1,2)$; the middle one lies on the line $y=x$ and is one of the $r-1=7$ subtracted points.

11. Empty columns and Pell's equation

A column can be empty. For even $a$ this never happens, because $v=0$ is admissible whenever $u^2\le L$. For odd $a$ every admissible $v$ is odd, so the column $u$ contains $v=\pm1$ unless $u^2=L$, that is, unless

$$u^2-2a^2=-1.$$

This is the negative Pell equation for $2$. Its positive solutions $(u,a)$ are $(1,1)$, $(7,5)$, $(41,29)$, $(239,169)$ and so on, and they come from the odd powers of $1+\sqrt2$; for instance $(1+\sqrt2)^3=7+5\sqrt2$. ^[6] For $r=4$ the only column $u=1$ is empty, which gives $N_B=0$ once more. For $r=20$ we have $a=5$ and $L=49$: the columns $u=1,3,5$ have $V(u)=5,5,3$, and the column $u=7$ is empty, $V(7)=-1$, so $N_B=2\,(6+6+4+0)=32$. The target value $a=r/4$ is even, but a program that also checks small radii must treat these columns correctly.

12. A pointer that only moves down

Each $V(u)$ could be found with its own integer square root, but a simpler sweep suffices. As $u$ grows, $L-u^2$ shrinks, so $V(u)$ never increases. We therefore keep one value $v$, start with $V$ of the first column, and in each further column lower $v$ in steps of $2$ until $u^2+v^2\le L$ holds again. If even the smallest candidate fails, $v=0$ for even $a$ or $v=1$ for odd $a$, the column is empty, and so is every later column, whose $u$ is larger; the sweep stops there. There are about $\sqrt L/2$ columns, and $v$ falls from about $\sqrt L$ towards $0$ in steps of $2$, so it is lowered at most about $\sqrt L/2$ times. The whole sweep takes about $\sqrt L\approx\sqrt2\,a$ steps, each a few multiplications and comparisons.

13. The complete formula

Collecting Sections 3, 5 and 6 gives

$$\begin{aligned}N(r)&=r^2+\frac r2+\frac{r(2r+1)}4+N_B-(r-1)\\&=\frac{3r^2}2-\frac r4+1+N_B.\end{aligned}$$

For $r=8$, the pieces $N_O=68$, $N_C=34$ and $N_B=5$ give $N(8)=68+34+5-7=100$, as in the statement. For $r=4$ the second line gives $24-1+1+0=24$. With $N_B=12$, $21$ and $32$ for $r=12$, $16$ and $20$, it gives $N(12)=226$, $N(16)=402$ and $N(20)=628$. All terms are integers because $4$ divides $r$.

How the Code Works

The C++ program ^[8] reads the radius from the option --r=, with $10^9$ as default, and parse_arguments rejects $r=0$ and every $r$ that is not a multiple of $4$, since Sections 3, 5 and 6 need $r/4$ to be an integer. isqrt_u64 returns $\lfloor\sqrt n\rfloor$ for $n\ge1$. It takes the floating-point square root as a first guess $g$ and corrects it with integer comparisons, raising $g$ while $g+1\le\lfloor n/(g+1)\rfloor$ and lowering it while $g\gt\lfloor n/g\rfloor$; the divisions avoid the overflow that squaring a 64-bit guess could cause. For $n=0$ the guess is $0$ and the second test divides by zero, which C++ leaves undefined. ^[9] This happens only for $r=4$, whose only column $u=1$ leaves $L-1=0$, so every default run meets it in its first checkpoint, before the main computation. On the arm64 computer used for this article the division gave $0$ and the checks passed, but the undefined-behaviour sanitizer flags it, and other machines need not get past this checkpoint.

count_circle_points_with_parity(a) prepares Section 8: limit is $L=2a^2-1$, parity is $a\bmod2$, the last column max_u is $\lfloor\sqrt L\rfloor$, lowered by one if its parity is wrong, and the first column first_u is $0$ or $1$. count_circle_u_range(limit, parity, begin_u, end_u) is the sweep of Section 12 over the columns from begin_u to end_u. It starts with $v=\lfloor\sqrt{L-u^2}\rfloor$ for the first column, lowered by one if its parity is wrong; when that is impossible because the root is $0$ and $a$ is odd, it sets $v$ to the sentinel kNoValue, the largest 64-bit number, which plays the role of $V=-1$. In every column it lowers $v$ by $2$ while $u^2+v^2\gt L$, switches to the sentinel when $v\le1$ still fails, stops at the sentinel, and otherwise adds $v+1$, doubled unless $u=0$. When there are at least $2\cdot10^6$ columns, count_circle_points_with_parity by default cuts them into consecutive blocks, each block runs its own sweep from its own integer square root, and the block counts are added; every column lies in exactly one block, so the total is the same.

solve_obtuse_count(r) assembles Section 13 with $a=r/4$: count_o is $r^2+r/2$, count_c is $r(2r+1)/4$, count_b comes from the sweep and degenerate is $r-1$, all held as signed 128-bit integers of type __int128_t, a compiler extension. ^[9] brute_force_count(r) is independent of Sections 4 to 12: it visits every point of the square $|x|,|y|\le r$, skips the points outside the diamond and those with $x=y$, and counts $B$ when one of the three products of Section 2 is negative. Unless --skip-checkpoints is given, run_checkpoints runs before the main computation. It checks $N(4)=24$ and $N(8)=100$ and compares the fast count with the brute force for $r=12,16,\dots,32$; the case $r=20$ contains the empty column of Section 11. Its last check computes $N(4\cdot10^6)$ once with the blocks forbidden and once with the default settings. With $707\,107$ columns, fewer than $2\cdot10^6$, the second call does not cut the range either, so this check repeats one computation instead of testing the blocks. Finally the program prints $N(r)$ and the elapsed time.

The Python and Java programs compute only $r=10^9$, without checkpoints, in one sweep. Python's solve contains the same count_circle_u_range, with math.isqrt as exact integer square root and $-1$ as the marker for "no value", and its integers never overflow. Java's main inlines the sweep in long arithmetic; its isqrt corrects the Math.sqrt guess by squaring, which is safe because the roots stay below $4\cdot10^8$. ^[10] All sizes fit into 64 bits: $L\lt1.25\cdot10^{17}$, $u^2+v^2\le2L$, $r(2r+1)\approx2\cdot10^{18}$, and $N(r)$ is at most $|S(r)|=2r^2+2r+1$, about $2\cdot10^{18}$, all below $2^{63}-1\approx9.22\cdot10^{18}$. This bound matters for Java, whose integer operations do not report overflow.

Complexity and Verification

Counting each operation on 64-bit numbers as one step, the closed formulas and the correction cost $O(1)$. The sweep visits about $\sqrt L/2\approx a/\sqrt2$ columns and lowers $v$ at most about as often, so it takes $O(a)=O(r)$ time and $O(1)$ memory; each block of the C++ program adds one integer square root and one stored partial count. For $r=10^9$ there are $176\,776\,696$ columns, and a single sweep lowers $v$ $176\,766\,455$ times. On the arm64 computer used for this article, the C++ program, compiled with -O2, and the Java program each finished the main computation in less than $0.2$ seconds, and the interpreted Python loop took about $30$ seconds.

Besides the checks built into the C++ program, we verified the following with programs that are not part of the published solutions:

For every $r=4,8,\dots,200$, the formula of Section 13 agrees with two brute forces: the dot-product test of Section 2, and a test with a cross product for collinearity and squared side lengths, since expanding $|\mathbf p-\mathbf q|^2$ shows that the angle between $\mathbf p$ and $\mathbf q$ is obtuse exactly when $|\mathbf p-\mathbf q|^2\gt|\mathbf p|^2+|\mathbf q|^2$. No point has two negative products.
The sweep agrees with a count that takes one integer square root per column for every $a\lt300$, and the bounds of Section 9 hold for every $a\le3000$.
For $r=4\cdot10^6$ the formula gives $N_B=1\,570\,796\,323\,109$ and $N(r)=25\,570\,795\,323\,110$, the value printed by the last checkpoint.
The three programs print the same value of $N(10^9)$; the C++ program, run with --skip-checkpoints so that the division of the checkpoint $r=4$ does not occur, gives it both with its default blocks and with the blocks forbidden, and the sanitizer reports nothing in that run.

Footnotes and References

The references below supply the background facts used above. The split into two half-planes and a disk, the correction for flat triangles, the column sweep and the bounds for the number of lattice points are derived explicitly in this article.

Project Euler 210 — Obtuse Angled Triangles. The official statement defines S(r) as the set of integer points (x, y) with |x| + |y| ≤ r, O = (0, 0), C = (r/4, r/4), and N(r) as the number of points B of S(r) for which the triangle OBC has an obtuse angle, its largest angle α satisfying 90° < α < 180°. It gives N(4) = 24 and N(8) = 100 and asks for N(1 000 000 000).
Wikipedia — Taxicab geometry. The taxicab distance between (x₁, y₁) and (x₂, y₂) is |x₁ − x₂| + |y₁ − y₂|, and in two dimensions a taxicab circle is a square oriented diagonally to the coordinate axes.
OEIS A001844 — Centered square numbers. The entry gives a(n) = 2n² + 2n + 1 = n² + (n + 1)², with the terms 1, 5, 13, 25, 41, …, and the comment that a(n) is the number of integer solutions (x, y) of |x| + |y| ≤ n.
Wikipedia — Dot product. The dot product is the sum of the products of corresponding coordinates, and it equals ‖a‖ ‖b‖ cos θ, where θ is the angle between a and b; two nonzero vectors are orthogonal exactly when their dot product is 0.
Wikipedia — Thales's theorem. If A, B and C are distinct points on a circle and AC is a diameter, then the angle ABC is a right angle. A related result states that the angle ABC is greater than 90° if B lies inside the circle and less than 90° if B lies outside.
Wikipedia — Pell's equation. The negative Pell equation x² − ny² = −1 is solvable for n = 2, with the fundamental solution (1, 1), and odd powers of the fundamental solution give an infinite tower of further solutions. OEIS A001653, the numbers k such that 2k² − 1 is a square, lists 1, 5, 29, 169, 985, …, the values of a with an empty column in Section 11.
Wikipedia — Gauss circle problem. The number of lattice points in a circle of radius r centred at the origin is πr² plus an error term that is small in comparison, because on average each unit square contains one lattice point; Gauss proved a bound on the error that grows linearly in r.
C++ — Euler210.cpp. The linked immutable C++ revision contains isqrt_u64, count_circle_u_range, count_circle_points_with_parity, solve_obtuse_count, brute_force_count and run_checkpoints. It is the implementation record for this article.
cppreference — Arithmetic operators. Built-in integer division by zero has undefined behaviour. The GCC manual — 128-bit integers describes the extension types __int128 and unsigned __int128 for targets with a 128-bit integer mode.
Python — math.isqrt returns the floor of the exact square root of a nonnegative integer. The Java Language Specification, §4.2.1 and §4.2.2: the values of long range from −9223372036854775808 to 9223372036854775807, and the integer operators do not indicate overflow or underflow in any way.

Problem 210 source code

C++

#include <algorithm>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <vector>

namespace {

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

constexpr u64 kDefaultR = 1'000'000'000ULL;
constexpr u64 kCheckpointR1 = 4ULL;
constexpr u64 kCheckpointExpected1 = 24ULL;
constexpr u64 kCheckpointR2 = 8ULL;
constexpr u64 kCheckpointExpected2 = 100ULL;
constexpr u64 kThreadConsistencyR = 4'000'000ULL;

struct Options {
    u64 r = kDefaultR;
    bool allow_multithreading = true;
    bool run_checkpoints = true;
    unsigned requested_threads = 0;
};

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

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

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

    value = parsed;
    return true;
}

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

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

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

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }

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

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

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

    if (options.r == 0ULL) {
        std::cerr << "--r must be >= 1.\n";
        return false;
    }
    if ((options.r % 4ULL) != 0ULL) {
        std::cerr << "This solver requires r to be divisible by 4.\n";
        return false;
    }

    return true;
}

u64 isqrt_u64(u64 n) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(n)));
    while ((r + 1ULL) <= n / (r + 1ULL)) {
        ++r;
    }
    while (r > n / r) {
        --r;
    }
    return r;
}

unsigned choose_thread_count(bool allow_multithreading,
                             unsigned requested_threads,
                             std::size_t workload) {
    if (!allow_multithreading || workload < 2'000'000ULL) {
        return 1;
    }

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

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

i128 count_circle_u_range(u64 limit, int parity, u64 begin_u, u64 end_u) {
    if (begin_u > end_u) {
        return 0;
    }

    constexpr u64 kNoValue = std::numeric_limits<u64>::max();

    u64 v = isqrt_u64(limit - begin_u * begin_u);
    if (static_cast<int>(v & 1ULL) != parity) {
        if (v == 0ULL) {
            v = kNoValue;
        } else {
            --v;
        }
    }

    i128 total = 0;
    for (u64 u = begin_u; u <= end_u; u += 2ULL) {
        if (v == kNoValue) {
            break;
        }

        const u64 u2 = u * u;
        while (u2 + v * v > limit) {
            if (v <= 1ULL) {
                v = kNoValue;
                break;
            }
            v -= 2ULL;
        }
        if (v == kNoValue) {
            break;
        }

        const i128 count_v = static_cast<i128>(v + 1ULL);
        const i128 multiplicity = (u == 0ULL ? static_cast<i128>(1) : static_cast<i128>(2));
        total += multiplicity * count_v;
    }

    return total;
}

u64 count_circle_points_with_parity(u64 a,
                                    bool allow_multithreading,
                                    unsigned requested_threads) {
    const u64 limit = 2ULL * a * a - 1ULL;
    const int parity = static_cast<int>(a & 1ULL);

    u64 max_u = isqrt_u64(limit);
    if (static_cast<int>(max_u & 1ULL) != parity) {
        --max_u;
    }

    const u64 first_u = static_cast<u64>(parity);
    if (first_u > max_u) {
        return 0ULL;
    }

    const u64 step_count = (max_u - first_u) / 2ULL + 1ULL;
    const unsigned threads =
        choose_thread_count(allow_multithreading, requested_threads, static_cast<std::size_t>(step_count));

    i128 total = 0;
    if (threads == 1U) {
        total = count_circle_u_range(limit, parity, first_u, max_u);
    } else {
        std::vector<std::thread> workers;
        std::vector<i128> partial(threads, 0);
        workers.reserve(threads);

        for (unsigned t = 0; t < threads; ++t) {
            const u64 step_begin = static_cast<u64>((static_cast<__int128>(step_count) * t) / threads);
            const u64 step_end_exclusive =
                static_cast<u64>((static_cast<__int128>(step_count) * (t + 1ULL)) / threads);

            if (step_begin >= step_end_exclusive) {
                continue;
            }

            const u64 begin_u = first_u + 2ULL * step_begin;
            const u64 end_u = first_u + 2ULL * (step_end_exclusive - 1ULL);
            workers.emplace_back([&, t, begin_u, end_u]() {
                partial[t] = count_circle_u_range(limit, parity, begin_u, end_u);
            });
        }

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

        for (const i128 value : partial) {
            total += value;
        }
    }

    return static_cast<u64>(total);
}

i128 solve_obtuse_count(u64 r,
                        bool allow_multithreading,
                        unsigned requested_threads) {
    const u64 a = r / 4ULL;

    // Angle at O is obtuse when x + y < 0.
    const i128 count_o = static_cast<i128>(r) * static_cast<i128>(r) +
                         static_cast<i128>(r / 2ULL);

    // Angle at C is obtuse when x + y > r/2.
    const i128 count_c =
        static_cast<i128>(r) * static_cast<i128>(2ULL * r + 1ULL) / static_cast<i128>(4);

    // Angle at B is obtuse when B is strictly inside the circle with diameter OC.
    const i128 count_b = static_cast<i128>(count_circle_points_with_parity(
        a,
        allow_multithreading,
        requested_threads));

    // Degenerate collinear points on y = x were included once above but must be excluded.
    const i128 degenerate = static_cast<i128>(r - 1ULL);

    return count_o + count_c + count_b - degenerate;
}

i128 brute_force_count(u64 r) {
    const i64 a = static_cast<i64>(r / 4ULL);

    i128 count = 0;
    for (i64 x = -static_cast<i64>(r); x <= static_cast<i64>(r); ++x) {
        for (i64 y = -static_cast<i64>(r); y <= static_cast<i64>(r); ++y) {
            if (std::llabs(x) + std::llabs(y) > static_cast<i64>(r)) {
                continue;
            }

            if (x == y) {
                // O, B, C are collinear -> largest angle is 180 degrees, excluded.
                continue;
            }

            const i64 sum = x + y;
            const i64 dot_o = a * sum;
            const i64 dot_b = x * x + y * y - a * sum;
            const i64 dot_c = 2LL * a * a - a * sum;

            if (dot_o < 0LL || dot_b < 0LL || dot_c < 0LL) {
                ++count;
            }
        }
    }

    return count;
}

std::string to_string_i128(i128 value) {
    if (value == 0) {
        return "0";
    }

    bool negative = value < 0;
    if (negative) {
        value = -value;
    }

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

    if (negative) {
        out.push_back('-');
    }
    std::reverse(out.begin(), out.end());
    return out;
}

bool run_checkpoints(const Options& options) {
    struct FixedCheckpoint {
        u64 r;
        u64 expected;
    };

    const std::vector<FixedCheckpoint> fixed = {
        {kCheckpointR1, kCheckpointExpected1},
        {kCheckpointR2, kCheckpointExpected2},
    };

    for (const FixedCheckpoint& cp : fixed) {
        const i128 got = solve_obtuse_count(cp.r, false, 1U);
        if (got != static_cast<i128>(cp.expected)) {
            std::cerr << "Checkpoint failed: N(" << cp.r << ") expected " << cp.expected
                      << ", got " << to_string_i128(got) << '\n';
            return false;
        }
        std::cout << "Checkpoint OK: N(" << cp.r << ") = " << cp.expected << '\n';
    }

    const std::vector<u64> brute_cases = {12ULL, 16ULL, 20ULL, 24ULL, 28ULL, 32ULL};
    for (const u64 r : brute_cases) {
        const i128 brute = brute_force_count(r);
        const i128 fast = solve_obtuse_count(r, false, 1U);
        if (brute != fast) {
            std::cerr << "Brute checkpoint failed: N(" << r << ") brute=" << to_string_i128(brute)
                      << ", fast=" << to_string_i128(fast) << '\n';
            return false;
        }
        std::cout << "Checkpoint OK: brute cross-check N(" << r << ") = "
                  << to_string_i128(fast) << '\n';
    }

    if (options.allow_multithreading) {
        const i128 single_thread = solve_obtuse_count(kThreadConsistencyR, false, 1U);
        const i128 multi_thread =
            solve_obtuse_count(kThreadConsistencyR, true, options.requested_threads);
        if (single_thread != multi_thread) {
            std::cerr << "Thread-consistency checkpoint failed at N(" << kThreadConsistencyR
                      << "): single=" << to_string_i128(single_thread)
                      << ", multi=" << to_string_i128(multi_thread) << '\n';
            return false;
        }
        std::cout << "Checkpoint OK: threaded consistency at N(" << kThreadConsistencyR
                  << ") = " << to_string_i128(single_thread) << '\n';
    }

    return true;
}

}  // namespace

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

    const auto start = std::chrono::steady_clock::now();

    if (options.run_checkpoints) {
        if (!run_checkpoints(options)) {
            return 1;
        }
    }

    const i128 answer =
        solve_obtuse_count(options.r, options.allow_multithreading, options.requested_threads);

    const auto finish = std::chrono::steady_clock::now();
    const std::chrono::duration<long double> elapsed = finish - start;

    std::cout << "Answer: " << to_string_i128(answer) << '\n';
    std::cout << "N(" << options.r << ") = " << to_string_i128(answer) << '\n';
    std::cout << "Elapsed: " << elapsed.count() << " s\n";

    return 0;
}

Python

import math

def solve():
    R = 1_000_000_000

    def isqrt_u64(n):
        r = math.isqrt(n)
        return r

    def count_circle_u_range(limit, parity, begin_u, end_u):
        if begin_u > end_u:
            return 0
        v = isqrt_u64(limit - begin_u * begin_u)
        if (v & 1) != parity:
            if v == 0:
                return 0
            v -= 1

        total = 0
        u = begin_u
        while u <= end_u:
            u2 = u * u
            while u2 + v * v > limit:
                if v <= 1:
                    v = -1
                    break
                v -= 2
            if v < 0:
                break
            count_v = v + 1
            multiplicity = 1 if u == 0 else 2
            total += multiplicity * count_v
            u += 2
        return total

    a = R // 4
    limit = 2 * a * a - 1
    parity = a & 1
    max_u = isqrt_u64(limit)
    if (max_u & 1) != parity:
        max_u -= 1
    first_u = parity

    count_b = count_circle_u_range(limit, parity, first_u, max_u)

    count_o = R * R + R // 2
    count_c = R * (2 * R + 1) // 4
    degenerate = R - 1

    answer = count_o + count_c + count_b - degenerate
    return str(answer)

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

Java

public class Euler210 {
    public static void main(String[] args) {
        long r = 1000000000L;
        long a = r / 4;
        long countO = r * r + r / 2;
        long countC = r * (2 * r + 1) / 4;
        // count_b: circle points with parity
        long limit = 2 * a * a - 1;
        int parity = (int) (a & 1);
        long maxU = isqrt(limit);
        if ((maxU & 1) != parity)
            maxU--;
        long firstU = parity;
        long totalB = 0;
        long v = isqrt(limit - firstU * firstU);
        if ((v & 1) != parity) {
            v--;
            if (v < 0)
                v = -1;
        }
        for (long u = firstU; u <= maxU; u += 2) {
            if (v < 0)
                break;
            long u2 = u * u;
            while (u2 + v * v > limit) {
                v -= 2;
                if (v < 0)
                    break;
            }
            if (v < 0)
                break;
            long countV = v + 1;
            long mult = (u == 0) ? 1 : 2;
            totalB += mult * countV;
        }
        long degenerate = r - 1;
        System.out.println(countO + countC + totalB - degenerate);
    }

    static long isqrt(long n) {
        if (n < 0)
            return 0;
        long s = (long) Math.sqrt((double) n);
        while (s * s > n)
            s--;
        while ((s + 1) * (s + 1) <= n)
            s++;
        return s;
    }
}