Problem 311: Biclinic Integral Quadrilaterals

Project Euler Problem 311 Solution

EulerSolve provides an optimized solution for Project Euler Problem 311, Biclinic Integral Quadrilaterals, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Write $$r=AO=CO,\qquad s=BO=DO=\frac{BD}{2},\qquad 0<r\le s.$$ Then every biclinic quadrilateral lives on one fixed pair of concentric circles centered at $O$: the vertices $A,C$ lie on radius $r$, while $B,D$ are the opposite endpoints of a diameter of radius $s$. The key invariant is $$AB^2+BC^2+CD^2+AD^2=4(r^2+s^2).$$ So if we define $$n=r^2+s^2,$$ then the condition in the problem is simply $n\le N/4$. The whole problem becomes: for each $n$, how many biclinic quadrilaterals are produced by the representations of $n$ as a sum of two squares? Mathematical Approach 1) Every side pair gives a strict representation of the same $n$. Look at vertex $A$. Let $$u=AB,\qquad v=AD,\qquad u<v.$$ Define $$x=\frac{u+v}{2},\qquad y=\frac{v-u}{2},\qquad u=x-y,\quad v=x+y.$$ In triangle $BAD$, the point $O$ is the midpoint of $BD$. Apollonius gives $$u^2+v^2=2(AO^2+BO^2)=2(r^2+s^2)=2n.$$ But also $$u^2+v^2=(x-y)^2+(x+y)^2=2(x^2+y^2).$$ Therefore $$x^2+y^2=n.$$ Because $u$ and $v$ are distinct positive integers, we get a strict representation $x>y>0$. The same argument at vertex $C$ gives another strict representation of the same $n$. 2) The reverse construction. Suppose we already know one representation $$n=s^2+r^2,\qquad s\ge r>0,$$ which will be used for the diagonal data $BD=2s$ and $AO=CO=r$....

Detailed mathematical approach

Problem Summary

Let $ABCD$ be a convex quadrilateral with integer sides $1\le AB\lt BC\lt CD\lt AD$. The diagonal $BD$ has integer length, $O$ is its midpoint, and $AO$ is an integer as well. Such a quadrilateral is a biclinic integral quadrilateral when $AO=CO\le BO=DO$. The statement exhibits one, with $AB=19$, $BC=29$, $CD=37$, $AD=43$, $BD=48$ and $AO=CO=23$. Writing $B(N)$ for the number of distinct biclinic integral quadrilaterals with $AB^2+BC^2+CD^2+AD^2\le N$, it records $B(10\,000)=49$ and $B(1\,000\,000)=38239$ and asks for $B(10^{10})$. ^[1]

The whole solution in one sentence: the pair of radii and the two pairs of sides are three different ways of writing one integer $t\le N/4$ as a sum of two squares, so $B(N)$ sums binomial coefficients in the number $m(t)$ of those ways — $\binom{m(t)}{3}$, plus $\binom{m(t)}{2}$ when the radii may be equal — and one sweep over lattice points produces every $m(t)$.

Mathematical Approach

1. Two concentric circles

Put $O$ at the origin. As $O$ is the midpoint of $BD$, the vertices $B$ and $D$ are the ends of a diameter of the circle of radius $R=BD/2$ about $O$. The condition $AO=CO$ puts $A$ and $C$ on one circle of radius $r=AO$ about the same centre, and the biclinic condition is exactly $r\le R$. In a convex quadrilateral each diagonal separates the other two vertices, so $A$ and $C$ lie strictly on opposite sides of the line $BD$. The two radii and the positions of $A$ and $C$ fix the whole figure.

2. One invariant: $t=r^2+R^2$

In the triangle $ABD$ the segment $AO$ is the median to the side $BD$. Apollonius' theorem says that a median $m$ to a side of length $a$ satisfies $b^2+c^2=2m^2+a^2/2$; equivalently, here, ^[2]

$$AB^2+AD^2=2AO^2+2BO^2=2r^2+2R^2 .$$

The same theorem in the triangle $CBD$, whose median to $BD$ is $CO$, gives $BC^2+CD^2=2r^2+2R^2$. Two things follow.

First, $BD$ must be even: $AB^2+AD^2=2r^2+BD^2/2$, and the left side and $2r^2$ are integers, so $BD^2/2$ is an integer and $BD$ is even. Both radii are therefore integers with $1\le r\le R$.

Second, adding the two identities,

$$AB^2+BC^2+CD^2+AD^2=4(r^2+R^2).$$

Write $t=r^2+R^2$. The size condition becomes $4t\le N$, that is $t\le T$ with $T=\lfloor N/4\rfloor$. In the example $r=23$ and $R=24$, so $t=529+576=1105$ and the four squares add up to $4420=4\cdot1105$.

3. A vertex is a representation of $t$

For the sides $AB=19$ and $AD=43$ of the example, the half-sum is $31$ and the half-difference is $12$, and $31^2+12^2=1105$, the $t$ just found. That is the whole trick. In general put $u=AB$ and $v=AD$, so $u\lt v$ and $u^2+v^2=2t$ by Section 2. This sum is even, so $u$ and $v$ have the same parity and

$$x=\frac{u+v}{2},\qquad y=\frac{v-u}{2}$$

are integers, with $u=x-y$, $v=x+y$ and

$$x^2+y^2=\frac{u^2+v^2}{2}=t .$$

Here $x\gt y$ because $u\ge1$, and $y\ge1$ because $u\ne v$. Conversely any integers $x\gt y\ge1$ with $x^2+y^2=t$ give $u=x-y\ge1$ and $v=x+y$ with $u^2+v^2=2t$, and the two constructions invert each other. So the candidate side pairs at $A$, before the triangle inequality of Section 5 is imposed, are exactly the strict representations

$$t=x^2+y^2,\qquad x\gt y\ge1 .$$

Vertex $C$, with its sides $BC$ and $CD$, gives a second strict representation of the same $t$. Let $\mathcal R(t)$ be the set of strict representations and $m(t)=|\mathcal R(t)|$. Two members of $\mathcal R(t)$ have different small parts, because $y$ determines $x=\sqrt{t-y^2}$. For the example $m(1105)=4$; Section 10 lists the four.

4. The diagonal is a representation too

The radii satisfy $t=R^2+r^2$ with $1\le r\le R$. If $r\lt R$ this is the strict representation $(R,r)\in\mathcal R(t)$. If $r=R$ then $t=2r^2$, so $t$ is twice a square and $(r,r)$ is a representation that is not strict.

5. Which representation may be the diagonal

Let $(x,y)\in\mathcal R(t)$ belong to $A$, so the triangle $ABD$ has sides $AB=x-y$, $AD=x+y$ and $BD=2R$. Vertex $A$ is strictly off the line $BD$ by Section 1, so this triangle is non-degenerate, and three lengths form a non-degenerate triangle exactly when each is strictly less than the sum of the other two. ^[3] Since $AD-AB=2y$ and $AB+AD=2x$, the requirement is

$$2y\lt 2R\lt 2x .$$

The left inequality is free: $y\lt x$ gives $2y^2\lt t$, and $r\le R$ gives $t\le 2R^2$, so $y\lt R$ always. The right one is a genuine restriction with a clean form, since $x^2=t-y^2$ and $R^2=t-r^2$:

$$x\gt R\iff y^2\lt r^2\iff y\lt r .$$

So the rule is: the representation supplying the diagonal must have a strictly larger small part than one supplying a vertex.

Conversely, when $y\lt r$ the vertex exists and lands on the inner circle by itself: if $PB=u$, $PD=v$ and $u^2+v^2=2t$, Apollonius in $PBD$ gives $2PO^2=u^2+v^2-2R^2=2r^2$, so $PO=r$. Section 7 gives its coordinates.

6. A three-element subset gives exactly one quadrilateral

Order $\mathcal R(t)$ by decreasing small part $y$. For $0\le a\lt b\lt t/2$,

$$b(t-b)-a(t-a)=(b-a)(t-a-b)\gt0,$$

so $s(t-s)$ grows with $s$ on $0\le s\lt t/2$. Strict representations have $y^2\lt t/2$, so taking $s=y^2$ shows that $(xy)^2=y^2(t-y^2)$, and hence $xy$, decreases along the order. From $(x\pm y)^2=t\pm2xy$ the sum $x+y$ then decreases and the difference $x-y$ increases along it.

Now take any three strict representations, labelled $P_1,P_2,P_3$ in that order, with small parts $y_1\gt y_2\gt y_3$. The diagonal must have the largest small part of the three, so it must be $P_1$: either other choice would leave $P_1$ as a vertex with the larger small part, which Section 5 forbids. Setting $r=y_1$, $R=x_1$ and reading the sides off $P_2$ and $P_3$,

$$AB=x_2-y_2,\qquad AD=x_2+y_2,$$

$$BC=x_3-y_3,\qquad CD=x_3+y_3,$$

the monotonicity above gives $AB\lt BC\lt CD\lt AD$ at once, and $AB\ge1$. So every three-element subset of $\mathcal R(t)$ yields exactly one biclinic integral quadrilateral: which of $P_2,P_3$ belongs to $A$ is forced by the order of the sides.

It runs the other way for a quadrilateral with $r\lt R$: that hands back $(R,r)$ and the representations at $A$ and at $C$, three different members of $\mathcal R(t)$ — the first by Section 5, the other two because $AB\lt BC$. Three-element subsets therefore match the quadrilaterals with unequal radii exactly, so such a $t$ contributes $\binom{m(t)}{3}$. The case $r=R$ is Section 8.

7. Convexity comes for free

Place $B=(-R,0)$ and $D=(R,0)$. A point $P=(p,q)$ with $PB=u$, $PD=v$ satisfies $u^2-v^2=PB^2-PD^2=4pR$, so

$$p=\frac{u^2-v^2}{4R},\qquad p^2+q^2=t-R^2=r^2 ,$$

the second identity again by Apollonius. So $q^2=r^2-p^2$, and a point with $q\ne0$ exists exactly when $|p|\lt r$ — the strict triangle inequality of Section 5. Take $A$ with $q\gt0$ and $C$ with $q\lt0$, putting them on opposite sides of $BD$.

The other diagonal then behaves as it must. The segment $AC$ lies inside the closed disc of radius $r$ about $O$, so it crosses the line $BD$ at distance at most $r$ from $O$ — not exactly $r$, since that would make the crossing point $A$ or $C$, which are off the line. It is thus nearer to $O$ than $r\le R$, hence strictly between $B$ and $D$. Each diagonal separates the other two vertices, so $ABCD$ is convex. Flipping the sign of $q$ at both vertices only reflects the figure, so each subset gives one quadrilateral up to congruence.

8. The exceptional case $t=2k^2$

When $t=2k^2$ the symmetric representation $(k,k)$ serves as the diagonal, giving $BD=2k$ and $AO=CO=BO=DO=k$; the problem allows this, demanding only $AO=CO\le BO=DO$. Its small part $k=\sqrt{t/2}$ exceeds that of every strict representation, so Section 5 excludes nothing: any two members of $\mathcal R(t)$ may be the vertices, and Section 6 fixes which is which. Such a $t$ contributes a further $\binom{m(t)}{2}$ quadrilaterals, with no double counting, since a diagonal representation is either strict or symmetric, never both.

9. The counting formula

Collecting Sections 6 and 8, with $\varepsilon(t)=1$ when $t$ is twice a square and $\varepsilon(t)=0$ otherwise,

$$\begin{aligned}B(N)&=\sum_{t\le T}\binom{m(t)}{3}\\&\quad+\sum_{t\le T}\varepsilon(t)\binom{m(t)}{2},\end{aligned}$$

where $T=\lfloor N/4\rfloor$ and $m(t)$ counts the partitions of $t$ into two distinct nonzero squares — the sequence A025441 of the OEIS. ^[6]

10. Worked examples

The example of the statement. Here $t=1105=5\cdot13\cdot17$ and

$$\mathcal R(1105)=\{(24,23),(31,12),(32,9),(33,4)\},$$

already in the order of Section 6. In the subset $\{(24,23),(31,12),(33,4)\}$ the diagonal is $(24,23)$, so $BD=48$ and $AO=CO=23$, while $(31,12)$ gives $AB=19$, $AD=43$ and $(33,4)$ gives $BC=29$, $CD=37$ — exactly the quadrilateral of the statement. Since $m(1105)=4$, this one value of $t$ contributes $\binom43=4$ quadrilaterals, with sides

$(19,23,41,43)$, $(19,29,37,43)$ and $(23,29,37,41)$, each with $BD=48$ and $AO=23$;
$(23,29,37,41)$ once more, with $BD=62$ and $AO=12$.

The last two share their four sides but differ in their diagonals, and both count: "distinct" in the statement means distinct quadrilaterals, not distinct quadruples of sides.

The smallest one. The least $t$ that contributes anything is $t=325=5^2\cdot13$, where $m(325)=3$ and

$$\mathcal R(325)=\{(15,10),(17,6),(18,1)\}.$$

Its single subset makes $(15,10)$ the diagonal, so $BD=30$ and $AO=CO=10$, while $(17,6)$ and $(18,1)$ give $AB=11$, $AD=23$ and $BC=17$, $CD=19$. The squares add up to $121+289+361+529=1300$, so $B(N)=0$ for $N\le1299$ and $B(1300)=1$.

Why the second sum is needed. Take $t=1250=2\cdot25^2$, whose only strict representations are $(31,17)$ and $(35,5)$, so $\binom23=0$. But $(25,25)$ is admissible as the diagonal, and the two give $AB=14$, $AD=48$, $BC=30$, $CD=40$ with $BD=50$ and $AO=CO=BO=DO=25$; here $196+900+1600+2304=5000$, four times $1250$. Without the $\varepsilon$ term this quadrilateral would be lost.

11. How large $m(t)$ gets, and how to compute it

No integer has to be factorised: for every $y=1,2,\dots$ and every $x\gt y$ with $x^2+y^2\le T$, add one to a counter at index $x^2+y^2$; at the end the counter at $t$ holds $m(t)$.

How wide must a counter be? Let $r_2(t)$ count the pairs $(a,b)\in\mathbb Z^2$ with $a^2+b^2=t$ and put $D(t)=r_2(t)/4$. If $t=2^{g}p_1^{f_1}\cdots q_1^{h_1}\cdots$ with primes $p_i\equiv1$ and $q_i\equiv3 \pmod 4$, then $D(t)=(f_1+1)(f_2+1)\cdots$ when every $h_i$ is even and $D(t)=0$ otherwise; among the primes, only $2$ and those $\equiv1\pmod4$ are themselves sums of two squares. ^[4] ^[5] Each strict representation accounts for eight of the pairs counted by $r_2$, a representation with $b=0$ for four, and one with $a=b$ for four, so

$$m(t)=\tfrac12\bigl(D(t)-\sigma(t)-\varepsilon(t)\bigr),$$

where $\sigma(t)=1$ exactly when $t$ is a square. To bound this on $t\le2.5\cdot10^9$: a power of $2$ or an even power of a prime $\equiv3\pmod4$ enlarges $t$ without changing $D(t)$, and an odd such power makes $D(t)=0$, so only primes $\equiv1\pmod4$ matter. For a fixed multiset of exponents the smallest $t$ puts the largest exponents on the smallest of them, $5,13,17,29,37,41,\dots$, and seven distinct ones are already too many:

$$\begin{aligned}&5\cdot13\cdot17\cdot29\cdot37\cdot41\cdot53\\&\qquad=2\,576\,450\,045 .\end{aligned}$$

Since $5^{14}$ is out of range too, it is enough to walk the $108$ non-increasing exponent vectors on those six primes that stay inside $2.5\cdot10^9$; their largest $D$ values, for one up to six primes, are $14$, $40$, $72$, $96$, $120$ and $128$. So $D(t)\le128$, attained at

$$t=5^3\cdot13\cdot17\cdot29\cdot37\cdot41=1\,215\,306\,625,$$

which is neither a square nor twice a square, so $m(t)=64$ there. Hence $m(t)\le64$ over the whole range, and a 16-bit counter has room to spare.

How the Code Works

The C++ program ^[9] reads its options in parse_arguments: the numeric --limit= and --block= go through parse_u64_after_prefix, which rejects 64-bit overflow, while --skip-checkpoints is handled directly, and a zero block width is refused. The defaults are $N=10^{10}$ and a block width $W=4\cdot10^6$.

solve_biclinic sets max_t to $\lfloor N/4\rfloor$ and cuts $1,\dots,T$ into $\lceil T/W\rceil$ blocks of width $W$ — $625$ for the target. The blocks exist for memory: a counter array for the whole range would need gigabytes, one block needs $W$ entries.

process_block(low, high, counts, touched) is the sweep of Section 11 restricted to one block; write $L$ and $H$ for its bounds. The outer loop over $y$ stops once $y^2+(y+1)^2$ exceeds $H$, the least value a strict representation with this $y$ can reach. The window of $x$ is $x\ge y+1$, raised to $\lceil\sqrt{L-y^2}\,\rceil$ by isqrt_ceil_u64 when the block starts higher, and capped at $\lfloor\sqrt{H-y^2}\rfloor$ by isqrt_u64. Both roots start from a long double estimate and step until the inequality holds exactly, so no rounding error survives. The inner loop carries $x^2$ along, advancing it by $2x+1$, so it never recomputes x*x. counts is a vector of 16-bit counters indexed by $n-L$; the first time one leaves zero its index joins touched, and a counter that would pass $65535$ raises an exception.

After the sweep the function walks touched alone — the values of $n$ that occurred — adding, for every $m\ge2$, n_choose_3(m) and, when is_twice_square(n) holds, n_choose_2(m): the two terms of Section 9. It resets each counter to zero, so one array serves every block. n_choose_3 forms $m(m-1)(m-2)$ in an unsigned 128-bit integer before dividing by six ^[8] — generous, since with $m\le64$ that product is at most $249\,984$.

Unless --skip-checkpoints is given, main first calls run_checkpoints, which recomputes the statement's $B(10\,000)=49$ and $B(1\,000\,000)=38239$ with solve_biclinic itself and stops on a mismatch. A third check recomputes $B(4\,000\,000)$ twice under different internal settings and demands the same value; at the default block width that range fits into one block, so the two runs coincide and the check only bites when --block= is much smaller. main prints $B(N)$ and the elapsed time, checkpoints included.

Euler311.java computes the same formula with isqrt, isqrtCeil, isTwiceSquare, nChoose2, nChoose3 and BlockTask.call, the counterpart of process_block. Several differences are real. Its solve fixes $N=10^{10}$ and $W=4\cdot10^6$ in the source and has no checkpoint routine. Each block allocates a fresh byte[] of length high-low+1, so there is no reset and no index list: the finalisation scans the whole array, not only the entries that occurred. A signed Java byte holds at most $127$, so the bound $m\le64$ keeps the counters safe, but nothing there would notice an overflow, whereas the C++ throws. And isqrt corrects a double estimate from Math.sqrt where the C++ starts from a long double one. All its arithmetic, nChoose3 included, is signed long.

Euler311.py does not reimplement the method. It is a bridge: it takes the problem number from its own file name, looks for solutionsCpp/Euler311.cpp under the parent of its own directory, and caches the executable beside that source. It compiles with clang++, or g++ if that is absent, at -std=c++17 -O2 when the executable is missing or older than the source, runs it without arguments, and returns what parse_output finds: a line labelled answer: if present, else the last line holding an equals sign — here the $B(N)$ line, since the closing timing line has none.

Complexity and Verification

Each step of the sweep handles one lattice point $(x,y)$ with $1\le y\lt x$ and $x^2+y^2\le T$, once, in the single block holding $x^2+y^2$. Those points fill one eighth of the disc of radius $\sqrt T$ without its boundary rays, so by the Gauss circle problem their number is $\pi T/8+O(\sqrt T)$. ^[7] At $T=2\,500\,000\,000$ the exact count is $981\,704\,953$ against $\pi T/8\approx981\,747\,704$, and the outer loop reaches $y=35\,354$.

Two smaller terms sit beside it. Every block restarts the outer loop at $y=1$, so a block ending at $H$ sets up $\Theta(\sqrt H)$ rows whether they contribute or not: over $K=\lceil T/W\rceil$ blocks, $\Theta(T^{3/2}/W)$ row set-ups, $14\,748\,772$ at the default width, about $1.5\%$ of the increments. The Java finalisation adds $\Theta(KW)=\Theta(T)$. The total $O(W+T+T^{3/2}/W)$ is linear in $N$ while the block count stays bounded, though not for a fixed $W$. Processing all $625$ blocks one after another takes about $5.6$ seconds on a current desktop machine, checkpoints skipped.

One block's workspace is $W$ counters — $8$ MB for the 16-bit C++ array, $4$ MB for the Java byte[] — plus the C++ index list, whose reserve of $W/4+1024$ indices of 32 bits (about $4$ MB) is an initial capacity that could grow to $W$; at the default width fewer than $0.9$ million distinct $n$ occur per block, so it does not. Java also keeps one Future per block, $O(K)$ bookkeeping. Storage is thus $O(W)$ in C++ and $O(W+K)$ in Java: it grows with the block width, not with $N$.

The arithmetic is exact. No counter overflows, by the bound $m\le64$. One block sum is at most

$$W\left(\binom{64}{3}+\binom{64}{2}\right)=4\cdot10^6\cdot43\,680,$$

under $1.8\cdot10^{11}$, so $625$ blocks stay below $1.1\cdot10^{14}$, far inside both 64-bit accumulators. Addition of non-negative integers is exact and associative, so the order in which block sums are combined cannot matter.

Beyond the packaged checkpoints, the following were computed for this article.

A brute force from the geometry — run over all $r\le R$, test $AD-AB\lt BD\lt AB+AD$ per candidate vertex, count admissible pairs — agrees with Section 9 for every $N\le3000$ and at $N=5000$, $10\,000$, $20\,000$, $50\,000$, $100\,000$ and $200\,000$.
All $3203$ three-element subsets of $\mathcal R(t)$ for $t\le40\,000$ were tried with each of the three choices of diagonal: exactly one always satisfies every condition, and it is always the representation with the largest small part, as Sections 5 and 6 claim.
$m(t)$ agrees with A025441 up to $t=20$ and with $\tfrac12(D-\sigma-\varepsilon)$ for every $t\lt60\,000$; an exhaustive search over the exponent vectors confirms $D\le128$.
The coordinates of Section 7 were formed for all four quadrilaterals of $t=1105$ and for those of $t=325$ and $t=1250$: each time the side lengths are the integers claimed and the four cross products of consecutive edges share one sign, so the figure is convex.
The C++ returns $B(1\,000\,000)=38239$ for block widths $1$, $7$, $999$, $123\,457$ and $4\,000\,000$, and $B(100\,000)=1652$. The Java port, whose limit is a compiled-in constant, reproduces $49$, $38239$ and $1652$ when that constant is set to the matching $N$.

Footnotes and References

The references below supply the background facts used above. The reduction to three representations of $t$, the rule that picks the diagonal, the bijection with three-element subsets and the bound $m(t)\le64$ are derived explicitly in this article.

Project Euler 311 — Biclinic Integral Quadrilaterals. The official statement takes a convex integer-sided quadrilateral ABCD with 1 ≤ AB < BC < CD < AD, integer BD with midpoint O and integer AO, calls it biclinic integral when AO = CO ≤ BO = DO, gives the example with sides 19, 29, 37, 43, diagonal 48 and AO = 23, and asks for B(10 000 000 000) after the checkpoints B(10 000) = 49 and B(1 000 000) = 38239.
Wikipedia — Apollonius's theorem. For the median AD to the side BC of a triangle, AB² + AC² = 2(BD² + AD²), equivalent, the article notes, to the parallelogram law. Section 2 applies it in the triangles ABD and CBD; Sections 5 and 7 apply it in reverse, so a vertex at the right distances from B and D lies on the inner circle.
Wikipedia — Triangle inequality. The article states the converse too: if each of three lengths is less than the sum of the other two a triangle of positive area exists, while equality gives a degenerate one of zero area. Section 5 needs the strict form, since the vertex lies off the line BD.
Wikipedia — Fermat's theorem on sums of two squares. An odd prime p is a sum of two squares if and only if p ≡ 1 (mod 4): 5, 13, 17 and 29 are, 3, 7, 11 and 19 are not. Hence the values of t with many representations — 325, 1105, the record 1 215 306 625 — are built from those primes.
Wikipedia — Sum of squares function. For n = 2^g p₁^f₁ ⋯ q₁^h₁ ⋯ with pᵢ ≡ 1 and qᵢ ≡ 3 (mod 4), the number of ordered pairs of integers whose squares sum to n is r₂(n) = 4(f₁ + 1)(f₂ + 1) ⋯ when every hᵢ is even and 0 otherwise. Section 11 turns this into the bound m(t) ≤ 64.
OEIS A025441. "Number of partitions of n into 2 distinct nonzero squares" — the quantity m(t) of Section 3. With offset 0 the terms begin 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, so a(5) = 1 from 4 + 1 and a(13) = 1 from 9 + 4.
Wikipedia — Gauss circle problem. The number of integer points with m² + n² ≤ r² is πr² + E(r), and Gauss proved |E(r) − 2| ≤ 2√2 πr for r ≥ √2, so the error is O(r). The complexity section uses this with r = √T.
GCC manual — 128-bit integers. The unsigned __int128 extension exists on targets with an integer mode wide enough for 128 bits, though constants of that type cannot be written where long long is narrower. The C++ uses it for the product inside n_choose_3.
C++ — Euler311.cpp. The linked immutable revision contains process_block, solve_biclinic, n_choose_3 and run_checkpoints, and matches the source described above byte for byte. It is the implementation record for this article; the brute-force comparisons above are extra checks, not packaged tests.

Problem 311 source code

C++

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

namespace {

using u64 = std::uint64_t;
using u32 = std::uint32_t;
using u16 = std::uint16_t;
using u128 = unsigned __int128;

constexpr u64 kDefaultLimit = 10'000'000'000ULL;
constexpr u64 kDefaultBlockSpan = 4'000'000ULL;
constexpr u64 kCheckpointLimit1 = 10'000ULL;
constexpr u64 kCheckpointExpected1 = 49ULL;
constexpr u64 kCheckpointLimit2 = 1'000'000ULL;
constexpr u64 kCheckpointExpected2 = 38'239ULL;
constexpr u64 kThreadConsistencyLimit = 4'000'000ULL;

struct Options {
    u64 limit = kDefaultLimit;
    u64 block_span = kDefaultBlockSpan;
    bool allow_multithreading = true;
    bool run_checkpoints = true;
    unsigned requested_threads = 0U;
};

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 = 0ULL;
    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 = 0ULL;
    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 parsed_u64 = 0ULL;
        if (parse_u64_after_prefix(arg, "--limit=", parsed_u64)) {
            options.limit = parsed_u64;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--block=", parsed_u64)) {
            options.block_span = parsed_u64;
            continue;
        }

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

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

    if (options.block_span == 0ULL) {
        std::cerr << "--block must be at least 1.\n";
        return false;
    }

    return true;
}

u64 isqrt_u64(u64 n) {
    if (n == 0ULL) {
        return 0ULL;
    }
    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;
}

u64 isqrt_ceil_u64(u64 n) {
    const u64 r = isqrt_u64(n);
    return (r * r == n) ? r : (r + 1ULL);
}

bool is_twice_square(u64 n) {
    if ((n & 1ULL) != 0ULL) {
        return false;
    }
    const u64 m = n >> 1ULL;
    const u64 r = isqrt_u64(m);
    return r * r == m;
}

u64 n_choose_2(u64 n) {
    return (n < 2ULL) ? 0ULL : (n * (n - 1ULL) / 2ULL);
}

u64 n_choose_3(u64 n) {
    if (n < 3ULL) {
        return 0ULL;
    }
    const u128 numerator = static_cast<u128>(n) * (n - 1ULL) * (n - 2ULL);
    return static_cast<u64>(numerator / 6ULL);
}

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

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

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

u64 process_block(u64 low,
                  u64 high,
                  std::vector<u16>& counts,
                  std::vector<u32>& touched) {
    touched.clear();

    for (u64 y = 1ULL;; ++y) {
        const u64 y2 = y * y;
        const u64 min_n = y2 + (y + 1ULL) * (y + 1ULL);
        if (min_n > high) {
            break;
        }

        u64 x_min = y + 1ULL;
        const u64 x_min_n = y2 + x_min * x_min;
        if (x_min_n < low) {
            x_min = isqrt_ceil_u64(low - y2);
            if (x_min <= y) {
                x_min = y + 1ULL;
            }
        }

        const u64 x_max = isqrt_u64(high - y2);
        if (x_max < x_min) {
            continue;
        }

        u64 x2 = x_min * x_min;
        for (u64 x = x_min; x <= x_max; ++x) {
            const u64 n = y2 + x2;
            const u32 idx = static_cast<u32>(n - low);
            u16& count = counts[static_cast<std::size_t>(idx)];
            if (count == 0U) {
                touched.push_back(idx);
            }
            if (count == std::numeric_limits<u16>::max()) {
                throw std::runtime_error("Representation count overflowed uint16_t.");
            }
            ++count;
            x2 += (2ULL * x + 1ULL);
        }
    }

    u64 block_sum = 0ULL;
    for (const u32 idx : touched) {
        u16& count = counts[static_cast<std::size_t>(idx)];
        const u64 m = static_cast<u64>(count);
        if (m >= 2ULL) {
            const u64 n = low + static_cast<u64>(idx);
            block_sum += n_choose_3(m);
            if (is_twice_square(n)) {
                block_sum += n_choose_2(m);
            }
        }
        count = 0U;
    }

    return block_sum;
}

u64 solve_biclinic(u64 limit,
                   bool allow_multithreading,
                   unsigned requested_threads,
                   u64 block_span) {
    if (limit < 4ULL) {
        return 0ULL;
    }

    const u64 max_t = limit / 4ULL;
    const u64 block_count_u64 = (max_t + block_span - 1ULL) / block_span;
    if (block_count_u64 > static_cast<u64>(std::numeric_limits<std::size_t>::max())) {
        throw std::runtime_error("Too many blocks for this platform.");
    }
    const std::size_t block_count = static_cast<std::size_t>(block_count_u64);
    const unsigned threads =
        choose_thread_count(allow_multithreading, requested_threads, block_count);

    std::vector<u64> partial(threads, 0ULL);
    auto worker = [&](unsigned tid) {
        std::vector<u16> counts(static_cast<std::size_t>(block_span), 0U);
        std::vector<u32> touched;
        touched.reserve(static_cast<std::size_t>(block_span / 4ULL + 1024ULL));

        u64 local_sum = 0ULL;
        for (std::size_t bi = tid; bi < block_count; bi += threads) {
            const u64 low = 1ULL + static_cast<u64>(bi) * block_span;
            const u64 high = std::min<u64>(max_t, low + block_span - 1ULL);
            local_sum += process_block(low, high, counts, touched);
        }
        partial[tid] = local_sum;
    };

    std::vector<std::thread> pool;
    pool.reserve(threads > 0U ? threads - 1U : 0U);
    for (unsigned t = 1U; t < threads; ++t) {
        pool.emplace_back(worker, t);
    }
    worker(0U);
    for (std::thread& thread : pool) {
        thread.join();
    }

    u64 total = 0ULL;
    for (const u64 x : partial) {
        total += x;
    }
    return total;
}

bool run_checkpoints(const Options& options) {
    const u64 sample1 = solve_biclinic(
        kCheckpointLimit1, options.allow_multithreading, options.requested_threads, options.block_span);
    if (sample1 != kCheckpointExpected1) {
        std::cerr << "Checkpoint failed for N=" << kCheckpointLimit1
                  << ": expected " << kCheckpointExpected1
                  << ", got " << sample1 << '\n';
        return false;
    }
    std::cout << "Checkpoint passed: B(" << kCheckpointLimit1 << ") = " << sample1 << '\n';

    const u64 sample2 = solve_biclinic(
        kCheckpointLimit2, options.allow_multithreading, options.requested_threads, options.block_span);
    if (sample2 != kCheckpointExpected2) {
        std::cerr << "Checkpoint failed for N=" << kCheckpointLimit2
                  << ": expected " << kCheckpointExpected2
                  << ", got " << sample2 << '\n';
        return false;
    }
    std::cout << "Checkpoint passed: B(" << kCheckpointLimit2 << ") = " << sample2 << '\n';

    if (options.allow_multithreading) {
        const u64 multi = solve_biclinic(
            kThreadConsistencyLimit, true, options.requested_threads, options.block_span);
        const u64 single = solve_biclinic(
            kThreadConsistencyLimit, false, options.requested_threads, options.block_span);
        if (multi != single) {
            std::cerr << "Thread consistency failed at N=" << kThreadConsistencyLimit
                      << ": multi-thread=" << multi
                      << ", single-thread=" << single << '\n';
            return false;
        }
        std::cout << "Checkpoint passed: thread consistency at N="
                  << kThreadConsistencyLimit << '\n';
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    try {
        const auto start_time = std::chrono::steady_clock::now();

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

        const u64 result = solve_biclinic(options.limit,
                                          options.allow_multithreading,
                                          options.requested_threads,
                                          options.block_span);

        const auto end_time = std::chrono::steady_clock::now();
        const std::chrono::duration<double> elapsed = end_time - start_time;

        std::cout << "B(" << options.limit << ") = " << result << '\n';
        std::cout << "Elapsed: " << elapsed.count() << " seconds\n";
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }

    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = subprocess.check_output([str(binary)], text=True)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

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

public class Euler311 {
    static long isqrt(long n) {
        long r = (long) Math.sqrt(n);
        while ((r + 1) * (r + 1) <= n)
            r++;
        while (r * r > n)
            r--;
        return r;
    }

    static long isqrtCeil(long n) {
        long r = isqrt(n);
        return (r * r == n) ? r : r + 1;
    }

    static boolean isTwiceSquare(long n) {
        if ((n & 1) != 0)
            return false;
        long m = n >> 1;
        long r = isqrt(m);
        return r * r == m;
    }

    static long nChoose2(long n) {
        return (n < 2) ? 0 : (n * (n - 1)) / 2;
    }

    static long nChoose3(long n) {
        return (n < 3) ? 0 : (n * (n - 1) * (n - 2)) / 6;
    }

    static class BlockTask implements Callable<Long> {
        long low, high;

        BlockTask(long low, long high) {
            this.low = low;
            this.high = high;
        }

        @Override
        public Long call() {
            int len = (int) (high - low + 1);
            byte[] counts = new byte[len];

            for (long y = 1;; ++y) {
                long y2 = y * y;
                long min_n = y2 + (y + 1) * (y + 1);
                if (min_n > high)
                    break;

                long xMin = y + 1;
                long xMinN = y2 + xMin * xMin;
                if (xMinN < low) {
                    xMin = isqrtCeil(low - y2);
                    if (xMin <= y)
                        xMin = y + 1;
                }

                long xMax = isqrt(high - y2);
                if (xMax < xMin)
                    continue;

                long x2 = xMin * xMin;
                for (long x = xMin; x <= xMax; ++x) {
                    long n = y2 + x2;
                    int idx = (int) (n - low);
                    counts[idx]++;
                    x2 += 2 * x + 1;
                }
            }

            long blockSum = 0;
            for (int idx = 0; idx < len; ++idx) {
                int count = counts[idx];
                if (count >= 2) {
                    long n = low + idx;
                    blockSum += nChoose3(count);
                    if (isTwiceSquare(n)) {
                        blockSum += nChoose2(count);
                    }
                }
            }

            return blockSum;
        }
    }

    public static String solve() {
        long limit = 10000000000L;
        long maxT = limit / 4;
        long blockSpan = 4000000;

        int numThreads = Runtime.getRuntime().availableProcessors();
        ExecutorService executor = Executors.newFixedThreadPool(numThreads);
        List<Future<Long>> futures = new ArrayList<>();

        long low = 1;
        while (low <= maxT) {
            long high = Math.min(maxT, low + blockSpan - 1);
            futures.add(executor.submit(new BlockTask(low, high)));
            low = high + 1;
        }

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

        return String.valueOf(total);
    }

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