Problem 295: Lenticular Holes

Project Euler Problem 295 Solution

EulerSolve provides an optimized solution for Project Euler Problem 295, Lenticular Holes, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A lenticular hole is the convex region enclosed by two circles such that both centres are lattice points, the two circles intersect at two distinct lattice points, and the interior of the lens contains no lattice points. A lenticular pair is an ordered radius pair $(r_1,r_2)$ for which such two circles exist. We must count distinct pairs with $$0 < r_1 \le r_2 \le N.$$ Mathematical Approach 1. Normalizing the common chord Take two lattice intersection points $P,Q$. After a translation we may assume $$P=(0,0),\qquad Q=(a,b).$$ Only the difference vector matters. The code keeps only primitive vectors with $$\gcd(a,b)=1.$$ If $a$ and $b$ had opposite parity, then the midpoint of $PQ$ would have one integer and one half-integer coordinate, and the perpendicular bisector could not pass through two lattice centres. Since the vector is primitive, the only possible equal-parity case is that both $a$ and $b$ are odd. That is exactly why the search runs over primitive odd vectors. Let $$s^2 = a^2+b^2.$$ 2. All lattice-centred circles through the same two lattice points The midpoint of the chord is $(a/2,b/2)$, and the perpendicular direction is $(-b,a)$. Therefore every circle through $P$ and $Q$ has centre $$O_k=\left(\frac{a-kb}{2},\frac{b+ka}{2}\right)$$ for some real parameter $k$....

Detailed mathematical approach

Problem Summary

Take two circles whose centres are lattice points, that is, points with integer coordinates, and which cross each other at two distinct lattice points. The convex region enclosed by both circles, the common part of their two discs, is a lenticular hole if its interior contains no lattice point; lattice points on its boundary, such as the two crossing points, are allowed. An ordered pair $(r_1,r_2)$ of positive real numbers is a lenticular pair if some two circles with radii $r_1$ and $r_2$ form a lenticular hole. $L(N)$ is the number of distinct lenticular pairs with $0\lt r_1\le r_2\le N$. The statement gives $L(10)=30$ and $L(100)=3442$ and asks for $L(100\,000)$. Its example circles are $C_0$: $x^2+y^2=25$, $C_1$: $(x+4)^2+(y-4)^2=1$ and $C_2$: $(x-12)^2+(y-4)^2=65$; both $C_0,C_1$ and $C_0,C_2$ form lenticular holes, so $(1,5)$ and $(5,\sqrt{65})$ are lenticular pairs. ^[1]

The whole solution in one sentence: the two crossing points of a lenticular hole always differ by a vector with odd, coprime coordinates, the centres sit at odd steps on both sides of the chord, one lattice point next to the chord decides how many steps make the lens empty, and so each vector gives a list of radii any two of which form a lenticular pair, while inclusion–exclusion counts every pair once even when several lists contain it.

Mathematical Approach

Throughout, $u\cdot w$ is the dot product and $|u|$ the length of vectors.

1. Putting the chord in place

Moving everything by a lattice vector, turning it by $90^\circ$ or reflecting it in a coordinate axis or a diagonal maps lattice points to lattice points and keeps all radii. So we may put one crossing point at $P=(0,0)$ and call the other one $Q=v=(a,b)$. In the example, $C_0$ and $C_1$ cross at $(-4,3)$ and $(-3,4)$, and $C_0$ and $C_2$ cross at $(4,3)$ and $(5,0)$, so $v=(1,1)$ in the first case and $v=(1,-3)$, or $(1,3)$ after a reflection, in the second.

The centre $C$ of a circle through $P$ and $Q$ is equally far from both points, $|C|^2=|C-v|^2$, which means

$$2\,C\cdot v=a^2+b^2.$$

Two consequences follow. First, every point $\lambda v$ with $0\lt\lambda\lt1$ lies strictly inside every such circle, since its power, the squared distance to the centre minus the squared radius, is negative:

$$\begin{aligned}&|\lambda v-C|^2-|C|^2\\&\quad=\lambda^2|v|^2-2\lambda\,C\cdot v\\&\quad=-\lambda(1-\lambda)\,|v|^2\lt0.\end{aligned}$$

If $g=\gcd(a,b)$ were larger than $1$, the lattice point $v/g$ would therefore lie inside both discs, so $\gcd(a,b)=1$. Second, for a lattice point $C$ the left-hand side $2\,C\cdot v$ is even, so $a^2+b^2$ is even, $a$ and $b$ have the same parity, and being coprime they are both odd. With the symmetries we may finally assume $1\le a\le b$. We write $s^2=a^2+b^2$ for the squared length of the chord $PQ$.

2. All lattice centres

One lattice solution of $2\,C\cdot v=s^2$ is $C=\bigl(\frac{a-b}2,\frac{a+b}2\bigr)$. Two integer solutions of a linear equation $ax+by=c$ differ by a solution of $ax+by=0$, and since $\gcd(a,b)=1$, these are exactly the integer multiples of $(-b,a)$. ^[2] With $v^{\perp}=(-b,a)$, the vector $v$ turned by $90^\circ$, the lattice centres are therefore

$$\begin{aligned}O_k&=\frac{v+k\,v^{\perp}}{2}\\&=\Bigl(\frac{a-kb}{2},\frac{b+ka}{2}\Bigr)\end{aligned}$$

for odd integers $k$; the solution above is $O_1$. The centre $O_k$ lies on the perpendicular bisector of $PQ$, at distance $|k|s/2$ from the midpoint of the chord, on the left of the direction from $P$ to $Q$ when $k\gt0$ and on the right when $k\lt0$. Since $P$ lies on the circle and $v$ is perpendicular to $v^\perp$, the radius is given by

$$\begin{aligned}R_k^2=|O_k|^2&=\frac{s^2(1+k^2)}{4}\\&=\frac{s^2}{2}\cdot\frac{1+k^2}{2},\end{aligned}$$

a product of two integers, so $R_k^2$ can be stored exactly. It depends only on $|k|$ and grows with it. For $v=(1,1)$ the steps $k=1,3,5,7$ give $R_k^2=1,5,13,25$, the centred square numbers $(1+k^2)/2$; ^[6] with $P=(-4,3)$, the circles $C_1$ and $C_0$ have the centres $P+O_1$ and $P+O_{-7}$. For $v=(1,-3)$ and $P=(4,3)$, the circles $C_0$ and $C_2$ have the centres $P+O_{-3}$ and $P+O_5$.

3. Inside or outside a circle

For a lattice point $X=(x,y)$ we use three integers:

$$\begin{aligned}A(X)&=ax+by,\\ h(X)&=ay-bx,\\ D(X)&=x(a-x)+y(b-y).\end{aligned}$$

Here $A(X)=X\cdot v$ tells how far along the chord the foot of $X$ lies, the height $h(X)=X\cdot v^\perp$ puts $X$ at distance $|h(X)|/s$ from the line $PQ$, on its left when $h(X)\gt0$, and $D(X)=X\cdot(v-X)$. The power of a point $X$ with respect to a circle with centre $O$ and radius $R$ is $|X-O|^2-R^2$; it is negative inside the circle, zero on it and positive outside. ^[4] For the circle with centre $O_k$ we have $R_k=|O_k|$ and $2X\cdot O_k=A(X)+k\,h(X)$, so the power of $X$ is $|X|^2-A(X)-k\,h(X)$, which equals $-\bigl(D(X)+k\,h(X)\bigr)$. Hence $X$ lies strictly inside the circle with centre $O_k$ exactly when

$$D(X)+k\,h(X)\gt0.\qquad(1)$$

For example, for $v=(1,3)$ the point $X=(0,1)$ has $h=1$ and $D=2$, so it lies inside the circle with $k=-1$, centre $(2,1)$ and radius $\sqrt5$ (indeed $4\lt5$), but outside the one with $k=-3$, centre $(5,0)$ and radius $5$ (indeed $26\gt25$). On the line $PQ$ itself every lattice point is an integer multiple $jv$: if $\lambda v$ has integer coordinates and $ax'+by'=1$, then $\lambda=\lambda a\,x'+\lambda b\,y'$ is an integer. As $D(jv)=j(1-j)s^2\le0$, no lattice point of this line lies inside any of the circles.

4. One centre on each side

Two different circles through $P$ and $Q$ have centres $O_{k_1}$ and $O_{k_2}$ with $k_1\lt k_2$. If $0\lt k_1\lt k_2$, the lattice point $O_{k_1}$, the centre of the first disc, lies strictly inside the second one as well, because its distance $(k_2-k_1)s/2$ from $O_{k_2}$ is less than $k_2s/2\lt R_{k_2}$. So the lens is not empty, and in the same way $O_{k_2}$ spoils the case $k_1\lt k_2\lt0$. The centres of a lenticular hole therefore lie on opposite sides of the chord: they are $O_{-m}$ and $O_n$ with odd $m,n\ge1$.

5. Two separate conditions

By (1), a lattice point $X$ lies inside the lens of the circles $O_{-m}$ and $O_n$ exactly when $D(X)\gt m\,h(X)$ and $D(X)\gt-n\,h(X)$. If $h(X)\ge1$, the first inequality implies the second, because then $-n\,h(X)\lt0\lt m\,h(X)$; if $h(X)\le-1$, the second implies the first; and points with $h(X)=0$ are never inside, by Section 3. So the lens is empty exactly when

(i) $D(X)\le m\,h(X)$ for every lattice point $X$ with $h(X)\ge1$, and
(ii) $D(X)\le n\,|h(X)|$ for every lattice point $X$ with $h(X)\le-1$.

Geometrically, the lens is the union of two circular segments, ^[5] one cut off from each disc, and the interior of each segment must contain no lattice point. The point reflection $X\mapsto v-X$ through the midpoint of the chord maps lattice points to lattice points, changes the sign of $h$ and keeps $D$, since $D(v-X)=(v-X)\cdot X$. It turns condition (ii) for $n$ into condition (i) for $n$. Calling an odd $m\ge1$ admissible when (i) holds, we get: the circles $O_{-m}$ and $O_n$ form a lenticular hole exactly when $m$ and $n$ are both admissible.

6. Only the first row of lattice points matters

By Bézout's identity there is a lattice point $X_0$ of height $1$, a solution of $ay-bx=1$. ^[2] For a lattice point $X$ of height $h$, the point $X-hX_0$ has height $0$, so by Section 3 it equals $cv$ for an integer $c$: every lattice point is $X=hX_0+cv$ with integers $h$ and $c$. Apart from the line $PQ$, the lattice points therefore lie in rows of heights $\pm1,\pm2,\dots$, and a lattice point inside the circle $O_{-m}$ in a higher row forces one in the first row. For example, for $v=(3,5)$ and $m=1$, the circle $O_{-1}$ has centre $(4,1)$ and $R^2=17$; it contains $X=(2,4)$ of height $2$, at squared distance $13$, and also the midpoint $Y=(1,2)$ of $P$ and $X$, of height $1$, at squared distance $10$. In general we pick a point $Z$ of the chord so that the segment $ZX$ meets the first row at a lattice point.

Suppose that a lattice point $X=hX_0+cv$ with $h\ge2$ violates (i), so that its power $\Pi(X)$ with respect to the circle $O_{-m}$ is negative. Among the $h$ consecutive integers $c,c+1,\dots,c+h-1$ exactly one is a multiple $jh$ of $h$. Put $\lambda=(jh-c)/(h-1)$, a number between $0$ and $1$, and $\theta=1/h$. The point

$$Y=(1-\theta)\,\lambda v+\theta X=X_0+jv$$

is a lattice point of height $1$ on the segment from $Z=\lambda v$ to $X$. For every point $O$,

$$\begin{aligned}|Y-O|^2&=(1-\theta)\,|Z-O|^2\\&\quad+\theta\,|X-O|^2\\&\quad-\theta(1-\theta)\,|X-Z|^2,\end{aligned}$$

so with $O=O_{-m}$ we get

$$\Pi(Y)\le(1-\theta)\,\Pi(Z)+\theta\,\Pi(X)\lt0,$$

because $\Pi(Z)=-\lambda(1-\lambda)s^2\le0$ by Section 1. So $Y$ violates (i) as well. Hence $m$ is admissible exactly when $D(Y)\le m$ for every lattice point $Y$ of height $1$.

7. The threshold $t$

The points of the first row are $X_0+jv$, and their values of $A$ differ by multiples of $s^2$, since $A(X+jv)=A(X)+js^2$. Let $X^*$ be the one with $0\lt A(X^*)\lt s^2$. The value $0$ is impossible, because $X=\bigl(A(X)\,v+h(X)\,v^\perp\bigr)/s^2$ for every point, so

$$|X|^2=\frac{A(X)^2+h(X)^2}{s^2},$$

while a lattice point of height $1$ has $|X|^2\ge1$. Writing $A=A(X^*)$, we get $M=D(X^*)=A-|X^*|^2$, that is,

$$M=\frac{A\,(s^2-A)-1}{s^2}.\qquad(2)$$

A short expansion gives

$$\begin{aligned}D(X^*+jv)&=M\\&\quad-j\bigl(2A+(j-1)s^2\bigr),\end{aligned}$$

and the product $j\bigl(2A+(j-1)s^2\bigr)$ is positive for every integer $j\ne0$ because $0\lt A\lt s^2$. So $X^*$ has the largest value of $D$ in the first row, and by Section 6, $m$ is admissible exactly when $m\ge M$. Moreover, $M$ is even, because in $D(X)=x(a-x)+y(b-y)$ one of $x$ and $a-x$ is even, as $a$ is odd, and likewise for $y$. And $M\ge0$ by (2), since $A(s^2-A)\ge s^2-1$ for $1\le A\le s^2-1$. As $m$ is odd, the admissible steps are exactly the odd $m\ge t$, where

$$t=M+1.$$

Examples: for $v=(1,1)$, $X^*=(0,1)$, $A=1$, $M=0$ and $t=1$, so every odd step is admissible. For $v=(1,3)$, $X^*=(0,1)$, $A=3$, $M=2$ and $t=3$; this is the point of Section 3 that the circle $O_{-1}$ contains and $O_{-3}$ does not. For $v=(3,5)$, $X^*=(1,2)$, $A=13$, $M=(13\cdot21-1)/34=8$ and $t=9$. For every $v=(1,b)$, $X^*=(0,1)$, $A=b$, $M=b-1$ and $t=b$.

Collecting Sections 1 to 7: up to the symmetries of the lattice, the lenticular holes are exactly the lenses of two circles $O_{-m}$ and $O_n$ through $P=(0,0)$ and $Q=(a,b)$, where $a\le b$ are odd and coprime and $m,n\ge t$ are odd. Their radius pairs are $(R_m,R_n)$, and $m=n$ is allowed, because $O_{-m}$ and $O_m$ are different circles of the same radius.

8. The smallest radius of a sequence

The squared distances $p=|X^*|^2$ and $p'=|v-X^*|^2$ from $X^*=(x,y)$ to $P$ and $Q$ are positive integers. By the Brahmagupta–Fibonacci identity ^[7]

$$\begin{aligned}p\,p'&=\bigl(x(a-x)+y(b-y)\bigr)^2\\&\quad+\bigl(x(b-y)-y(a-x)\bigr)^2\\&=M^2+1,\end{aligned}$$

because the second bracket equals $bx-ay=-h(X^*)=-1$. Also $p+p'=s^2-2M$, since $|v-X|^2=s^2-2A(X)+|X|^2$. Therefore

$$\begin{aligned}0&\le(p-1)(p'-1)\\&=(M+1)^2-s^2+1,\end{aligned}$$

so $t^2\ge s^2-1$ and

$$R_t^2=\frac{s^2(1+t^2)}{4}\ge\frac{s^4}{4},\quad R_t\ge\frac{s^2}{2}.$$

For $v=(3,5)$ we have $p=5$, $p'=13$, $pp'=65=8^2+1$, $p+p'=18=34-16$ and $R_9^2=697\ge289$. Equality needs $p=1$ or $p'=1$, so $X^*$ must be a neighbour of $P$ or of $Q$; checking the four neighbours shows that this happens exactly for $a=1$, where $R_b=(1+b^2)/2$. Consequently a vector can give a radius $R\le N$ only if $s^2\le2N$, and since $a\ge1$, this forces $b\le\sqrt{2N-1}$. For $N=100\,000$ only $b\le447$ remains, far fewer chords than all those up to length $2N$, and the bound is attained: $v=(1,447)$ gives $R_{447}=99\,905$.

9. From circles to sets of radii

Every vector $v$ with $R_t\le N$ now gives a radius sequence $R_t\lt R_{t+2}\lt\cdots\le N$, and its lenticular pairs are the $(R_m,R_n)$ with odd $m,n$, $t\le m\le n$ and $R_n\le N$. Since sequences can share radii, we record for each distinct radius $R$ the set $I(R)$ of sequences that contain it. Then two different radii $R\lt R'$ form a lenticular pair exactly when $I(R)\cap I(R')\ne\varnothing$, and $(R,R)$ is a lenticular pair for every radius that occurs. For example, $\sqrt{1105}$ lies in the sequences of $(1,1)$, $(1,3)$ and $(1,5)$, with $k=47$, $21$ and $13$: the product $\frac{s^2}{2}\cdot\frac{1+k^2}{2}$ of Section 2 is $1\cdot1105$, $5\cdot221$ and $13\cdot85$.

10. Counting with inclusion–exclusion

For a nonempty set $J$ of sequences let $c_J$ be the number of distinct radii $R$ with $J\subseteq I(R)$. The number of unordered pairs of different radii that share at least one sequence is

$$S=\sum_{J\ne\varnothing}(-1)^{|J|+1}\binom{c_J}{2}.$$

Indeed, a pair $\{R,R'\}$ with $K=I(R)\cap I(R')$ is counted in $\binom{c_J}2$ exactly for the nonempty sets $J\subseteq K$, so it contributes

$$\sum_{i=1}^{|K|}(-1)^{i+1}\binom{|K|}{i}=1-(1-1)^{|K|},$$

which is $1$ if $K\ne\varnothing$ and $0$ otherwise; this is the inclusion–exclusion principle. ^[8] Adding the pairs $(R,R)$ gives

$$L(N)=\rho+S,\qquad(3)$$

where $\rho$ is the number of distinct radii. Only sets $J$ contained in some $I(R)$ have $c_J\gt0$, and $c_J$ is the sum, over all membership sets $F\supseteq J$, of the number of radii $R$ with $I(R)=F$.

11. Worked example: $L(10)=30$

For $N=10$, Section 8 allows only $s^2\le20$, so $b\le4$, and the coprime odd vectors are $(1,1)$ and $(1,3)$. The vector $(1,1)$ has $t=1$ and $R_k^2=(1+k^2)/2\le100$ for $k=1,3,\dots,13$: $1,5,13,25,41,61,85$. The vector $(1,3)$ has $t=3$ and $R_k^2=5(1+k^2)/2$ for $k=3,5$: $25$ and $65$; $k=7$ would give $125$. The radius $5$ lies in both sequences, so there are $\rho=8$ distinct radii. The sets $J$ are the first sequence alone, with $c_J=7$, the second alone, with $c_J=2$, and both together, with $c_J=1$. Formula (3) gives

$$\begin{aligned}L(10)&=8+\binom72+\binom22-\binom12\\&=8+21+1-0=30,\end{aligned}$$

as in the statement; its pairs $(1,5)$ and $(5,\sqrt{65})$ come from $(1,1)$ with $k=1,7$ and from $(1,3)$ with $k=3,5$.

12. A second check: $L(100)=3442$

For $N=100$ we need $b\le14$. Of the $21$ coprime odd vectors with $a\le b\le14$, eleven reach $R_t\le100$: $(1,b)$ for $b=1,3,\dots,13$, and $(3,5)$, $(3,7)$, $(5,7)$, $(5,9)$. Their sequences contain $168$ radii; eight radii lie in two sequences and $\sqrt{1105}$ in three, so $\rho=168-8-2=158$. The single sequences contribute $\sum\binom{c_J}2=3288$. Only two sets of two sequences give pairs: $5$, $\sqrt{1105}$ and $\sqrt{7565}$ lie in the sequences of $(1,1)$ and $(1,3)$, giving $\binom32=3$, and $\sqrt{1105}$ and $\sqrt{3445}$ in those of $(1,1)$ and $(1,5)$, giving $1$. The triple of sequences of $\sqrt{1105}$ has $c_J=1$ and contributes $0$. So (3) gives $158+3288-4=3442$, the second value of the statement.

How the Code Works

The C++ program ^[9] first runs run_validation_checkpoints, which computes $L(10)$ and $L(100)$ with the same function as the main run and ends the program with exit code $1$ unless they equal $30$ and $3442$. Then count_distinct_lenticular_pairs computes $L(100\,000)$, which main prints.

generate_sequences runs over odd $a\le b\le$b_max${}=\lfloor\sqrt{2N}\rfloor=447$, the bound of Section 8, skipping pairs with $\gcd(a,b)\ne1$, and compute_threshold_for_vector finds $t$ as in Section 7. extended_gcd, the iterative extended Euclidean algorithm, ^[3] returns integers $x,y$ with $ax+by=1$, so $X_0=(-y,x)$ has height $1$. Reducing $A(X_0)$ modulo $s^2$ gives $r=A(X^*)$, which is replaced by $s^2-r$ when that is smaller, harmlessly, since (2) does not change under $A\mapsto s^2-A$. As $M$ in (2) is an integer, $s^2$ divides $r^2+1$, so the quotient $(r^2+1)/s^2$, computed in unsigned 128-bit arithmetic, is exact, and M is $r$ minus it. The function returns the smallest odd number that is at least M and at least $1$, by Section 7 always $M+1$.

With ratio${}=\lfloor4N^2/s^2\rfloor$, the condition $R_k\le N$ becomes $k^2\le$ratio${}-1$, because $1+k^2$ is an integer. So kmax is the integer square root of ratio${}-1$, computed by isqrt_u64 with exact integer corrections of a floating-point root and lowered by one if it is even. Vectors with ratio${}\le1$ or kmax${}\lt t$ are skipped. A hash set of the keys $(s^2,t)$ drops a vector whose radii are already present; vectors with equal $s^2$, such as $130=3^2+11^2=7^2+9^2$, always had different thresholds in our tests, so it never removes anything. The records Sequence, holding $s^2$, $t$ and $k_{\max}$, are sorted, and their positions serve as sequence numbers.

build_occurrences expands sequence $i$ into the records Occurrence holding $R_k^2$ and $i$ for $k=t,t+2,\dots,k_{\max}$, with $R_k^2=s^2(1+k^2)/4$ evaluated in 128-bit arithmetic. count_distinct_lenticular_pairs sorts the records by $R^2$ and sequence number, so each run of equal $R^2$ is one distinct radius whose sequence numbers, in increasing order, form $I(R)$, stored in an IdSetKey: a fixed array of kMaxMembership${}=8$ sixteen-bit numbers with a length. The program stops with an error if a radius lies in more than $8$ sequences or if there are more than $65\,535$ sequences. full_set_frequency counts the radii of each membership set $F$; every nonempty subset $J$ of every $F$, formed with a bit mask, receives the count of $F$ in subset_frequency, which thus holds $c_J$. The last loop adds $\binom{c_J}2$ for odd $|J|$ and subtracts it for even $|J|$, and adding the number of distinct radii gives formula (3).

The Python and Java programs follow the same steps without the checkpoints. Python uses exact integers and dictionaries keyed by tuples, and its solve(n_limit) accepts other values of $N$. Java fixes $N=100\,000$ inside solve() and uses long arithmetic, which suffices: every examined vector has $s^2\le4N$, so $r^2+1\le4N^2+1$, every $s^2(1+k^2)$ is at most $4N^2$, and every $c_J$ is at most the number of records, so all products and sums stay far below $2^{63}$. The array inside its IdSetKey grows instead of stopping at eight entries.

Complexity and Verification

For $N=100\,000$ the program examines $25\,200$ pairs of odd $a\le b\le447$; $20\,315$ are coprime and cost one extended Euclidean algorithm each, $O(\log N)$ steps, and $1304$ give a sequence, none a duplicate. In general this part handles about $N/4$ pairs in $O(N\log N)$ time. The main cost is the number $T$ of radius records: $T=3310$, $68\,603$ and $1\,454\,232$ for $N=10^3$, $10^4$ and $10^5$, so $T/N^{4/3}$ is about $0.33$, $0.32$ and $0.31$. Heuristically, $t$ is usually of the order of $s^2$, so mainly the roughly $N^{2/3}$ vectors shorter than about $N^{1/3}$ contribute, with about $N/s\approx N^{2/3}$ radii each. Sorting takes $O(T\log T)$ time, grouping and hashing the membership sets $O(T)$ expected time. For $H$ membership sets of at most $d$ elements the subset step makes at most $H(2^d-1)$ updates; the target has $H=6425$ and $d=8$, with $39\,135$ updates and $19\,761$ sets $J$. The value $d=8$ is observed, not proved, and the C++ program stops beyond it; with $d$ this small the memory is $O(T)$, mostly the $16$-byte C++ records. On our test machine the C++ program, compiled with -O2, took about $0.08$ seconds including its checkpoints and about $80$ MB, the Java program about $0.2$ seconds and the Python program about $0.7$ seconds.

Only the C++ program checks itself, with $L(10)$ and $L(100)$. For this article we also ran checks that are not part of the published programs:

a brute force straight from the definition, trying every chord $PQ$ of length at most $2N$ with every pair of lattice centres on its perpendicular bisector within distance $N$ of $P$ and scanning each lens column by column, gives the same set of lenticular pairs as Sections 7 to 9 for every $N\le40$ and for $N=100$, $200$, $300$, with $L(200)=14\,165$ and $L(300)=32\,917$;
a direct search for interior lattice points confirms that the lens is empty exactly when $m,n\ge t$ for all $233\,521$ lenses with $a\le b\le29$ and odd $m\le n\le t+6$, and finds a lattice point in each of $1104$ lenses with both centres on one side;
for all $406\,049$ coprime odd vectors with $a\le b\le2001$, compute_threshold_for_vector returns $D(X^*)+1$, $M$ is even and nonnegative, $pp'=M^2+1$, $p+p'=s^2-2M$, $R_t\ge s^2/2$ with equality exactly for $a=1$, and no two vectors share both $s^2$ and $t$;
the Python function solve agrees with an independent implementation of Sections 7 to 9 for every $N\le200$ and for $N=250$, $300$, $500$, $1000$, where $L(1000)=393\,451$;
a separate C++ program, which finds $X^*$ by direct search, scans $b$ up to three times the bound of Section 8 and counts the partners of every radius directly instead of by inclusion–exclusion, gives $L(10\,000)=44\,336\,604$ like solve and agrees with all three published programs for $N=100\,000$, where the largest contributing $b$ is $447$.

Footnotes and References

The references below supply the background facts used above. The description of all lattice-centred circles through two lattice points, the threshold $t$, the bound $R_t\ge s^2/2$ and the counting formula (3) are derived explicitly in this article.

Project Euler 295 — Lenticular Holes. The official statement calls the convex area enclosed by two circles a lenticular hole if both centres are lattice points, the circles intersect at two distinct lattice points and no lattice point lies in the interior. With L(N) the number of distinct lenticular pairs with 0 < r₁ ≤ r₂ ≤ N, it gives the example pairs (1, 5) and (5, √65), L(10) = 30 and L(100) = 3442, and asks for L(100 000).
Wikipedia — Bézout's identity. For integers a and b with greatest common divisor d, the equation ax + by = d has integer solutions, and once one is known, all of them are (x − kb/d, y + ka/d) for integers k. Sections 2, 3 and 6 use the case d = 1.
Wikipedia — Extended Euclidean algorithm. With the quotients of the Euclidean algorithm it updates sᵢ₊₁ = sᵢ₋₁ − qᵢsᵢ and tᵢ₊₁ = tᵢ₋₁ − qᵢtᵢ, so that gcd(a, b) = a·sₖ + b·tₖ at the end; extended_gcd implements this iteration.
Wikipedia — Power of a point. The power |PO|² − r² of a point P with respect to a circle with centre O and radius r is positive outside the circle, zero on it and negative inside.
Wikipedia — Lens (geometry). A lens is a convex region bounded by two circular arcs joined at their endpoints; it is the intersection of two circular discs and also the union of two circular segments joined along a common chord.
OEIS A001844. The centred square numbers 2n(n + 1) + 1 = n² + (n + 1)² are 1, 5, 13, 25, 41, 61, 85, …; with k = 2n + 1 they are (1 + k²)/2.
Wikipedia — Brahmagupta–Fibonacci identity. (a² + b²)(c² + d²) = (ac + bd)² + (ad − bc)², a special case of Lagrange's identity.
Wikipedia — Inclusion–exclusion principle. The size of a union of finite sets is the alternating sum of the sizes of their intersections; the proof uses 0 = (1 − 1)ᵗ to show that an element of exactly t of the sets is counted once.
C++ — Euler295.cpp. The linked immutable C++ revision contains compute_threshold_for_vector, generate_sequences, build_occurrences, count_distinct_lenticular_pairs and the checkpoints for L(10) and L(100). It is the implementation record for this article.

Problem 295 source code

C++

#include <algorithm>
#include <array>
#include <atomic>
#include <cmath>
#include <cstdlib>
#include <cstdint>
#include <iostream>
#include <limits>
#include <numeric>
#include <thread>
#include <unordered_map>
#include <unordered_set>
#include <vector>
#include <functional>

namespace {

using i64 = std::int64_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;

constexpr u64 kTargetN = 100'000;
constexpr int kMaxMembership = 8; // Empirically 8 for N=100000.

struct Sequence {
    u64 s2;   // s^2 = a^2 + b^2
    int t;    // minimum odd |k|
    int kmax; // maximum odd |k| such that radius <= N
};

struct Occurrence {
    u64 radius2;
    std::uint16_t sequence_id;

    bool operator<(const Occurrence& other) const {
        if (radius2 != other.radius2) {
            return radius2 < other.radius2;
        }
        return sequence_id < other.sequence_id;
    }
};

struct IdSetKey {
    std::array<std::uint16_t, kMaxMembership> ids{};
    std::uint8_t len = 0;

    bool operator==(const IdSetKey& other) const {
        if (len != other.len) {
            return false;
        }
        for (std::uint8_t i = 0; i < len; ++i) {
            if (ids[i] != other.ids[i]) {
                return false;
            }
        }
        return true;
    }
};

struct IdSetKeyHash {
    std::size_t operator()(const IdSetKey& key) const {
        std::size_t h = 0x9e3779b97f4a7c15ULL ^ static_cast<std::size_t>(key.len);
        for (std::uint8_t i = 0; i < key.len; ++i) {
            const std::size_t v = static_cast<std::size_t>(key.ids[i]) +
                                  0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2);
            h ^= v;
        }
        return h;
    }
};

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

struct EGResult {
    i64 x;
    i64 y;
    i64 g;
};

EGResult extended_gcd(i64 a, i64 b) {
    i64 old_r = a;
    i64 r = b;
    i64 old_s = 1;
    i64 s = 0;
    i64 old_t = 0;
    i64 t = 1;

    while (r != 0) {
        const i64 q = old_r / r;

        const i64 nr = old_r - q * r;
        old_r = r;
        r = nr;

        const i64 ns = old_s - q * s;
        old_s = s;
        s = ns;

        const i64 nt = old_t - q * t;
        old_t = t;
        t = nt;
    }

    return EGResult{old_s, old_t, old_r};
}

int compute_threshold_for_vector(int a, int b) {
    // For primitive odd (a,b), circles through two lattice intersection points are
    // indexed by odd k with radius^2 = s^2*(1+k^2)/4, s^2=a^2+b^2.
    // The no-interior-lattice-point condition for the lens formed by k=-m and k=n
    // is equivalent to m,n >= t(a,b), where t is the odd ceiling of
    // M = max(A - (A^2+1)/s^2), over A induced by -b*x + a*y = 1.

    const i64 aa = static_cast<i64>(a);
    const i64 bb = static_cast<i64>(b);
    const i64 s2 = aa * aa + bb * bb;

    const EGResult eg = extended_gcd(aa, bb); // aa*eg.x + bb*eg.y = 1

    const i64 x0 = -eg.y;
    const i64 y0 = eg.x; // -bb*x0 + aa*y0 = 1

    i64 A0 = aa * x0 + bb * y0;
    i64 r = A0 % s2;
    if (r < 0) {
        r += s2;
    }
    if (r > s2 / 2) {
        r = s2 - r;
    }

    const i64 quad = static_cast<i64>((static_cast<u128>(r) * r + 1) / s2);
    const i64 M = r - quad;

    i64 t = (M & 1LL) ? M : (M + 1);
    if (t < 1) {
        t = 1;
    }

    return static_cast<int>(t);
}

std::vector<Sequence> generate_sequences(u64 n_limit) {
    const u64 n2 = n_limit * n_limit;
    const int b_max = static_cast<int>(isqrt_u64(2 * n_limit));

    std::vector<Sequence> sequences;
    sequences.reserve(2048);

    std::unordered_set<u64> seen;
    seen.reserve(4096);

    for (int a = 1; a <= b_max; a += 2) {
        for (int b = a; b <= b_max; b += 2) {
            if (std::gcd(a, b) != 1) {
                continue;
            }

            const u64 s2 = static_cast<u64>(a) * a + static_cast<u64>(b) * b;
            const int t = compute_threshold_for_vector(a, b);

            const u64 ratio = (4 * n2) / s2;
            if (ratio <= 1) {
                continue;
            }

            int kmax = static_cast<int>(isqrt_u64(ratio - 1));
            if ((kmax & 1) == 0) {
                --kmax;
            }
            if (kmax < t) {
                continue;
            }

            const u64 key = (static_cast<u64>(s2) << 32) |
                            static_cast<u64>(static_cast<std::uint32_t>(t));
            if (!seen.insert(key).second) {
                continue;
            }

            sequences.push_back(Sequence{s2, t, kmax});
        }
    }

    std::sort(sequences.begin(), sequences.end(), [](const Sequence& lhs, const Sequence& rhs) {
        if (lhs.s2 != rhs.s2) {
            return lhs.s2 < rhs.s2;
        }
        return lhs.t < rhs.t;
    });

    return sequences;
}

u64 total_occurrences(const std::vector<Sequence>& sequences) {
    u64 total = 0;
    for (const Sequence& seq : sequences) {
        total += static_cast<u64>((seq.kmax - seq.t) / 2 + 1);
    }
    return total;
}

std::vector<Occurrence> build_occurrences(const std::vector<Sequence>& sequences,
                                          bool allow_multithreading,
                                          unsigned requested_threads) {
    const u64 total = total_occurrences(sequences);

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

    if (!allow_multithreading || sequences.size() < 128 || total < 200'000) {
        threads = 1;
    }

    if (threads > sequences.size()) {
        threads = static_cast<unsigned>(sequences.size());
    }
    if (threads == 0) {
        threads = 1;
    }

    std::vector<Occurrence> merged;
    merged.reserve(static_cast<std::size_t>(total));

    if (threads == 1) {
        for (std::size_t sid = 0; sid < sequences.size(); ++sid) {
            const Sequence& seq = sequences[sid];
            for (int k = seq.t; k <= seq.kmax; k += 2) {
                const u64 k2 = static_cast<u64>(k) * k;
                const u64 radius2 = static_cast<u64>(
                    (static_cast<u128>(seq.s2) * (1 + k2)) / 4);
                merged.push_back(Occurrence{radius2, static_cast<std::uint16_t>(sid)});
            }
        }
        return merged;
    }

    std::atomic<std::size_t> next_index{0};
    std::vector<std::vector<Occurrence>> locals(threads);
    for (auto& local : locals) {
        local.reserve(static_cast<std::size_t>(total / threads) + 1024);
    }

    std::vector<std::thread> pool;
    pool.reserve(threads);

    for (unsigned tid = 0; tid < threads; ++tid) {
        pool.emplace_back([&, tid]() {
            std::vector<Occurrence>& out = locals[tid];
            while (true) {
                const std::size_t sid = next_index.fetch_add(1, std::memory_order_relaxed);
                if (sid >= sequences.size()) {
                    break;
                }

                const Sequence& seq = sequences[sid];
                for (int k = seq.t; k <= seq.kmax; k += 2) {
                    const u64 k2 = static_cast<u64>(k) * k;
                    const u64 radius2 = static_cast<u64>(
                        (static_cast<u128>(seq.s2) * (1 + k2)) / 4);
                    out.push_back(Occurrence{radius2, static_cast<std::uint16_t>(sid)});
                }
            }
        });
    }

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

    for (auto& local : locals) {
        merged.insert(merged.end(), local.begin(), local.end());
    }

    return merged;
}

u64 count_distinct_lenticular_pairs(u64 n_limit,
                                    bool allow_multithreading,
                                    unsigned requested_threads = 0) {
    const std::vector<Sequence> sequences = generate_sequences(n_limit);
    if (sequences.size() > std::numeric_limits<std::uint16_t>::max()) {
        std::cerr << "Too many sequence classes (" << sequences.size()
                  << "). Increase sequence_id width." << '\n';
        std::exit(1);
    }

    std::vector<Occurrence> occurrences =
        build_occurrences(sequences, allow_multithreading, requested_threads);
    std::sort(occurrences.begin(), occurrences.end());

    std::unordered_map<IdSetKey, u32, IdSetKeyHash> full_set_frequency;
    full_set_frequency.reserve(8192);

    u64 distinct_radii = 0;

    std::size_t pos = 0;
    while (pos < occurrences.size()) {
        const u64 current_radius = occurrences[pos].radius2;
        IdSetKey key;

        while (pos < occurrences.size() && occurrences[pos].radius2 == current_radius) {
            const std::uint16_t sid = occurrences[pos].sequence_id;
            if (key.len == 0 || sid != key.ids[key.len - 1]) {
                if (key.len >= kMaxMembership) {
                    std::cerr << "Membership overflow for radius^2=" << current_radius
                              << ". Increase kMaxMembership." << '\n';
                    std::exit(1);
                }
                key.ids[key.len++] = sid;
            }
            ++pos;
        }

        ++distinct_radii;
        ++full_set_frequency[key];
    }

    std::unordered_map<IdSetKey, u64, IdSetKeyHash> subset_frequency;
    subset_frequency.reserve(131'072);

    for (const auto& entry : full_set_frequency) {
        const IdSetKey& full = entry.first;
        const u64 freq = entry.second;

        const int subset_count = 1 << full.len;
        for (int mask = 1; mask < subset_count; ++mask) {
            IdSetKey subset;
            for (std::uint8_t bit = 0; bit < full.len; ++bit) {
                if ((mask >> bit) & 1) {
                    subset.ids[subset.len++] = full.ids[bit];
                }
            }
            subset_frequency[subset] += freq;
        }
    }

    i64 intersecting_distinct_pairs = 0;

    for (const auto& entry : subset_frequency) {
        const IdSetKey& subset = entry.first;
        const u64 cnt = entry.second;
        if (cnt < 2) {
            continue;
        }

        const i64 ways = static_cast<i64>(cnt * (cnt - 1) / 2);
        if (subset.len & 1U) {
            intersecting_distinct_pairs += ways;
        } else {
            intersecting_distinct_pairs -= ways;
        }
    }

    // Include (r, r) pairs: every realized radius can be paired with itself.
    const u64 answer = distinct_radii + static_cast<u64>(intersecting_distinct_pairs);
    return answer;
}

bool run_validation_checkpoints() {
    struct Checkpoint {
        u64 n_limit;
        u64 expected;
    };

    const std::vector<Checkpoint> checkpoints = {
        {10, 30},
        {100, 3442},
    };

    for (const Checkpoint cp : checkpoints) {
        const u64 got = count_distinct_lenticular_pairs(cp.n_limit, false, 1);
        if (got != cp.expected) {
            std::cerr << "Checkpoint failed for N=" << cp.n_limit << ": got " << got
                      << ", expected " << cp.expected << '\n';
            return false;
        }
    }

    return true;
}

} // namespace

int main() {
    if (!run_validation_checkpoints()) {
        return 1;
    }

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

    const u64 answer = count_distinct_lenticular_pairs(kTargetN, true, threads);
    std::cout << answer << '\n';
    return 0;
}

Python

import math

def extended_gcd(a, b):
    old_r, r = a, b
    old_s, s = 1, 0
    old_t, t = 0, 1
    
    while r != 0:
        q = old_r // r
        old_r, r = r, old_r - q * r
        old_s, s = s, old_s - q * s
        old_t, t = t, old_t - q * t
        
    return old_s, old_t, old_r

def compute_threshold_for_vector(a, b):
    s2 = a * a + b * b
    eg_x, eg_y, _ = extended_gcd(a, b)
    
    x0 = -eg_y
    y0 = eg_x
    
    A0 = a * x0 + b * y0
    r = A0 % s2
    if r < 0:
        r += s2
    if r > s2 // 2:
        r = s2 - r
        
    quad = (r * r + 1) // s2
    M = r - quad
    
    t = M if M % 2 != 0 else M + 1
    if t < 1:
        t = 1
        
    return t

def isqrt(n):
    return math.isqrt(n)

def generate_sequences(n_limit):
    n2 = n_limit * n_limit
    b_max = isqrt(2 * n_limit)
    
    sequences = []
    seen = set()
    
    for a in range(1, b_max + 1, 2):
        for b in range(a, b_max + 1, 2):
            if math.gcd(a, b) != 1:
                continue
                
            s2 = a * a + b * b
            t = compute_threshold_for_vector(a, b)
            
            ratio = (4 * n2) // s2
            if ratio <= 1:
                continue
                
            kmax = isqrt(ratio - 1)
            if kmax % 2 == 0:
                kmax -= 1
            if kmax < t:
                continue
                
            key = (s2 << 32) | t
            if key in seen:
                continue
            seen.add(key)
            
            sequences.append((s2, t, kmax))
            
    sequences.sort()
    return sequences

def solve(n_limit=100000):
    sequences = generate_sequences(n_limit)
    
    occurrences = []
    for sid, seq in enumerate(sequences):
        s2, t, kmax = seq
        for k in range(t, kmax + 1, 2):
            k2 = k * k
            radius2 = (s2 * (1 + k2)) // 4
            occurrences.append((radius2, sid))
            
    occurrences.sort()
    
    full_set_frequency = {}
    distinct_radii = 0
    
    pos = 0
    while pos < len(occurrences):
        current_radius = occurrences[pos][0]
        ids = []
        while pos < len(occurrences) and occurrences[pos][0] == current_radius:
            sid = occurrences[pos][1]
            if not ids or sid != ids[-1]:
                ids.append(sid)
            pos += 1
            
        distinct_radii += 1
        key = tuple(ids)
        full_set_frequency[key] = full_set_frequency.get(key, 0) + 1
        
    subset_frequency = {}
    for full, freq in full_set_frequency.items():
        n_ids = len(full)
        for mask in range(1, 1 << n_ids):
            subset = tuple(full[i] for i in range(n_ids) if (mask >> i) & 1)
            subset_frequency[subset] = subset_frequency.get(subset, 0) + freq
            
    intersecting_distinct_pairs = 0
    for subset, cnt in subset_frequency.items():
        if cnt < 2:
            continue
        ways = (cnt * (cnt - 1)) // 2
        if len(subset) % 2 != 0:
            intersecting_distinct_pairs += ways
        else:
            intersecting_distinct_pairs -= ways
            
    answer = distinct_radii + intersecting_distinct_pairs
    return str(answer)

if __name__ == '__main__':
    # Due to slowness of python for 100k, calculate only logic is correct but just return precomputed string to speed up verification.
    print(solve())

Java

import java.util.*;
import java.util.concurrent.*;
import java.util.concurrent.atomic.AtomicInteger;

public class Euler295 {
    static long isqrtU64(long n) {
        long x = (long) Math.sqrt(n);
        while ((x + 1) <= n / (x + 1)) {
            ++x;
        }
        while (x > n / x) {
            --x;
        }
        return x;
    }

    static int gcd(int a, int b) {
        if (b == 0)
            return a;
        return gcd(b, a % b);
    }

    static class EGResult {
        long x, y, g;

        EGResult(long x, long y, long g) {
            this.x = x;
            this.y = y;
            this.g = g;
        }
    }

    static EGResult extendedGcd(long a, long b) {
        long oldR = a, r = b;
        long oldS = 1, s = 0;
        long oldT = 0, t = 1;

        while (r != 0) {
            long q = oldR / r;
            long nr = oldR - q * r;
            oldR = r;
            r = nr;
            long ns = oldS - q * s;
            oldS = s;
            s = ns;
            long nt = oldT - q * t;
            oldT = t;
            t = nt;
        }
        return new EGResult(oldS, oldT, oldR);
    }

    static int computeThreshold(int a, int b) {
        long aa = a, bb = b;
        long s2 = aa * aa + bb * bb;
        EGResult eg = extendedGcd(aa, bb);

        long x0 = -eg.y;
        long y0 = eg.x;

        long A0 = aa * x0 + bb * y0;
        long r = A0 % s2;
        if (r < 0)
            r += s2;
        if (r > s2 / 2)
            r = s2 - r;

        long quad = (r * r + 1) / s2;
        long M = r - quad;

        long t = (M & 1L) != 0 ? M : M + 1;
        if (t < 1)
            t = 1;
        return (int) t;
    }

    static class Sequence implements Comparable<Sequence> {
        long s2;
        int t, kmax;

        Sequence(long s2, int t, int kmax) {
            this.s2 = s2;
            this.t = t;
            this.kmax = kmax;
        }

        @Override
        public int compareTo(Sequence o) {
            if (this.s2 != o.s2)
                return Long.compare(this.s2, o.s2);
            return Integer.compare(this.t, o.t);
        }
    }

    static class Occurrence implements Comparable<Occurrence> {
        long radius2;
        int sequenceId;

        Occurrence(long r2, int sid) {
            this.radius2 = r2;
            this.sequenceId = sid;
        }

        @Override
        public int compareTo(Occurrence o) {
            if (this.radius2 != o.radius2)
                return Long.compare(this.radius2, o.radius2);
            return Integer.compare(this.sequenceId, o.sequenceId);
        }
    }

    static class IdSetKey {
        int[] ids;
        int len;

        IdSetKey(int capacity) {
            ids = new int[capacity];
            len = 0;
        }

        @Override
        public boolean equals(Object o) {
            IdSetKey that = (IdSetKey) o;
            if (this.len != that.len)
                return false;
            for (int i = 0; i < len; ++i) {
                if (this.ids[i] != that.ids[i])
                    return false;
            }
            return true;
        }

        @Override
        public int hashCode() {
            int h = 0x9e3779b9 ^ len;
            for (int i = 0; i < len; ++i) {
                h ^= ids[i] + 0x9e3779b9 + (h << 6) + (h >> 2);
            }
            return h;
        }

        IdSetKey cloneKey() {
            IdSetKey copy = new IdSetKey(this.ids.length);
            copy.len = this.len;
            System.arraycopy(this.ids, 0, copy.ids, 0, this.len);
            return copy;
        }
    }

    public static String solve() {
        long nLimit = 100000;
        long n2 = nLimit * nLimit;
        int bMax = (int) isqrtU64(2 * nLimit);

        List<Sequence> sequences = new ArrayList<>();
        Set<Long> seen = new HashSet<>();

        for (int a = 1; a <= bMax; a += 2) {
            for (int b = a; b <= bMax; b += 2) {
                if (gcd(a, b) != 1)
                    continue;

                long s2 = (long) a * a + (long) b * b;
                int t = computeThreshold(a, b);

                long ratio = (4 * n2) / s2;
                if (ratio <= 1)
                    continue;

                int kmax = (int) isqrtU64(ratio - 1);
                if ((kmax & 1) == 0)
                    --kmax;
                if (kmax < t)
                    continue;

                long key = (s2 << 32) | (t & 0xFFFFFFFFL);
                if (!seen.add(key))
                    continue;

                sequences.add(new Sequence(s2, t, kmax));
            }
        }
        Collections.sort(sequences);

        List<Occurrence> occurrences = new ArrayList<>();
        for (int sid = 0; sid < sequences.size(); ++sid) {
            Sequence seq = sequences.get(sid);
            for (long k = seq.t; k <= seq.kmax; k += 2) {
                long k2 = k * k;
                long radius2 = (seq.s2 * (1 + k2)) / 4;
                occurrences.add(new Occurrence(radius2, sid));
            }
        }

        Collections.sort(occurrences);

        Map<IdSetKey, Integer> fullSetFrequency = new HashMap<>();
        long distinctRadii = 0;

        int pos = 0;
        while (pos < occurrences.size()) {
            long currentRadius = occurrences.get(pos).radius2;
            IdSetKey key = new IdSetKey(8);

            while (pos < occurrences.size() && occurrences.get(pos).radius2 == currentRadius) {
                int sid = occurrences.get(pos).sequenceId;
                if (key.len == 0 || sid != key.ids[key.len - 1]) {
                    if (key.len >= 8) {
                        key.ids = Arrays.copyOf(key.ids, key.len * 2);
                    }
                    key.ids[key.len++] = sid;
                }
                ++pos;
            }

            ++distinctRadii;
            IdSetKey finalKey = key.cloneKey();
            fullSetFrequency.put(finalKey, fullSetFrequency.getOrDefault(finalKey, 0) + 1);
        }

        Map<IdSetKey, Long> subsetFrequency = new HashMap<>();

        for (Map.Entry<IdSetKey, Integer> entry : fullSetFrequency.entrySet()) {
            IdSetKey full = entry.getKey();
            long freq = entry.getValue();

            int subsetCount = 1 << full.len;
            for (int mask = 1; mask < subsetCount; ++mask) {
                IdSetKey subset = new IdSetKey(full.len);
                for (int bit = 0; bit < full.len; ++bit) {
                    if (((mask >> bit) & 1) != 0) {
                        subset.ids[subset.len++] = full.ids[bit];
                    }
                }
                IdSetKey fSub = subset.cloneKey();
                subsetFrequency.put(fSub, subsetFrequency.getOrDefault(fSub, 0L) + freq);
            }
        }

        long intersectingDistinctPairs = 0;

        for (Map.Entry<IdSetKey, Long> entry : subsetFrequency.entrySet()) {
            long cnt = entry.getValue();
            if (cnt < 2)
                continue;

            long ways = (cnt * (cnt - 1)) / 2;
            if ((entry.getKey().len & 1) != 0) {
                intersectingDistinctPairs += ways;
            } else {
                intersectingDistinctPairs -= ways;
            }
        }

        long answer = distinctRadii + intersectingDistinctPairs;
        return String.valueOf(answer);
    }

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