Problem 354: Distances in a Bee's Honeycomb

Project Euler Problem 354 Solution

EulerSolve provides an optimized solution for Project Euler Problem 354, Distances in a Bee's Honeycomb, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The centers of the honeycomb cells form a triangular lattice. In axial coordinates, every displacement is described by two integers $a,b$, and its squared Euclidean length is $$L^2 = 3(a^2+ab+b^2).$$ Therefore the distance-counting function can be rewritten as $$B(L)=\left|\left\{(a,b)\in\mathbb Z^2:3(a^2+ab+b^2)=L^2\right\}\right|.$$ If $L^2/3\notin\mathbb Z$, then $B(L)=0$. Otherwise, with $$n=\frac{L^2}{3},$$ the problem reduces to the representation count $$r(n)=\left|\left\{(a,b)\in\mathbb Z^2:a^2+ab+b^2=n\right\}\right|.$$ The local implementations count all $n\le N$ with $r(n)=450$, where the target repository limit is $$N=\left\lfloor\frac{(5\cdot 10^{11})^2}{3}\right\rfloor.$$ Each such $n$ corresponds to exactly one honeycomb distance $L=\sqrt{3n}$, so this arithmetic reformulation is equivalent to the original geometric question. Mathematical Approach Step 1: The Honeycomb Distance Becomes an Eisenstein Norm The quadratic form $$a^2+ab+b^2$$ is the norm form of the Eisenstein integers, since $$a^2+ab+b^2=N(a-b\omega),\qquad \omega=e^{2\pi i/3}.$$ This is why the arithmetic of primes modulo $3$ controls the number of lattice representations. The C++ function count_representations_formula implements exactly this standard norm-counting formula....

Detailed mathematical approach

Problem Summary

A honeycomb tiles the plane with regular hexagons of side $1$, one cell holds the queen bee, and the comb is as large as we need. Distances are measured centre to centre. For a positive real number $L$, $\mathrm{B}(L)$ is the number of cells whose centre lies at distance exactly $L$ from the queen's cell. The statement gives $\mathrm{B}(\sqrt3)=6$, $\mathrm{B}(\sqrt{21})=12$ and $\mathrm{B}(111111111)=54$, and asks how many $L\le5\cdot10^{11}$ satisfy $\mathrm{B}(L)=450$. ^[1]

The whole solution in one sentence: the cell centres form a triangular lattice, so $\mathrm{B}(L)$ counts the integer pairs with $a^2+ab+b^2=L^2/3$, a count that is $6$ times a product read off the prime factorisation, so $\mathrm{B}(L)=450$ pins that factorisation down to four exponent patterns and we count the integers below the limit carrying one of them.

Mathematical Approach

1. The centres form a triangular lattice

A regular hexagon of side $1$ has apothem $\sqrt3/2$, the distance from its centre to the middle of an edge. ^[2] Two cells that share an edge therefore have their centres $2\cdot\sqrt3/2=\sqrt3$ apart. Joining the centres of neighbouring cells produces equilateral triangles, so the centres form a triangular lattice. ^[3]

Take two steps to neighbouring cells $60^\circ$ apart as a basis $u,v$. Then $|u|=|v|=\sqrt3$ and $u\cdot v=3\cos60^\circ=3/2$, so every cell centre is $au+bv$ for a unique integer pair $(a,b)$, at squared distance

$$|au+bv|^2=3a^2+3ab+3b^2.$$

Building the tiling explicitly confirms this. The positive squared centre distances up to $120$ are

$$\begin{aligned}&3,\ 9,\ 12,\ 21,\ 27,\ 36,\ 39,\ 48,\\&57,\ 63,\ 75,\ 81,\ 84,\ 93,\ 108,\ 111,\ 117,\end{aligned}$$

exactly the values $3(a^2+ab+b^2)$.

Hence $\mathrm{B}(L)=0$ unless $L^2=3n$ for a positive integer $n$, and in that case

$$\mathrm{B}(L)=r(n),\qquad n=\frac{L^2}{3},$$

where $r(n)$ counts the integer pairs with $a^2+ab+b^2=n$. As $L=\sqrt{3n}$ grows with $n$, the $L$ we want correspond one to one with the integers $n$ we want.

2. Small cases by hand

For $n=1$ the solutions are $(1,0)$, $(0,1)$, $(-1,0)$, $(0,-1)$, $(1,-1)$ and $(-1,1)$: the six cells touching the queen, so $\mathrm{B}(\sqrt3)=r(1)=6$. For $n=7$ there are twelve, among them $(1,2)$, $(2,1)$, $(3,-1)$ and $(-1,3)$, so $\mathrm{B}(\sqrt{21})=r(7)=12$. Up to $n=57$, grouped by the value of $r$, the $n$ with $r(n)\ne0$ are

$$r=6:\ 1,3,4,9,12,16,25,27,36,48,$$

$$r=12:\ 7,13,19,21,28,31,37,39,43,52,57,$$

together with $n=49$, whose $r(49)=18$ is the first value beyond $12$. Two patterns stand out. Every value is a multiple of $6$, and many $n$ are missing altogether. These counts are the coefficients of the theta series of the hexagonal lattice. ^[4]

3. The Eisenstein integers explain both patterns

Let $\omega=e^{2\pi i/3}$, so $\omega^2+\omega+1=0$, and let $\mathbb Z[\omega]$ be the ring of Eisenstein integers $x+y\omega$ with $x,y\in\mathbb Z$. Its norm is $N(x+y\omega)=x^2-xy+y^2=|x+y\omega|^2$, and putting $y=-b$ gives

$$a^2+ab+b^2=N(a-b\omega).$$

So $r(n)$ counts the Eisenstein integers of norm $n$. The ring is a Euclidean domain, hence has unique factorisation, and its elements of norm $1$ are exactly the six units $\pm1,\pm\omega,\pm\omega^2$, that is the elements with an inverse in the ring. ^[5] Multiplying by a unit preserves the norm, and $uz=vz$ with $z\ne0$ forces $u=v$, so the six units carry each solution to six distinct solutions. That is why $6$ divides every $r(n)$: it is the six-fold symmetry of the honeycomb.

Two elements differing by a unit factor are called associates; unique factorisation fixes the prime factors only up to order and associates. Write $A(n)=r(n)/6$ for the number of solutions up to associates. The second pattern comes from how ordinary primes behave in $\mathbb Z[\omega]$. ^[5]

A prime $p\equiv1\pmod 3$ splits: $p=\pi\bar\pi$ with $N(\pi)=N(\bar\pi)=p$, the bar being complex conjugation. For example $7=(3+\omega)(2-\omega)$, both factors of norm $7$; their twelve unit multiples are the twelve solutions counted by $r(7)$. The two are never associates: running through the six units in $\bar\pi=u\pi$ forces $p=3$ or $p$ a square.
A prime $q\equiv2\pmod 3$ stays prime, with norm $q^2$. No Eisenstein integer has norm $q$ at all.
The prime $3$ ramifies: $\lambda=1-\omega$ has $N(\lambda)=3$, and $3$ is a unit times $\lambda^2$.

4. The representation formula

Write $P_1$ for the primes $\equiv1\pmod 3$ and $P_2$ for those $\equiv2\pmod 3$, and factor

$$n=3^{e_0}\prod_{p\in P_1}p^{e_p}\prod_{q\in P_2}q^{f_q}.$$

Claim. If some $f_q$ is odd then $r(n)=0$; otherwise

$$r(n)=6\prod_{p\in P_1}(e_p+1).$$

Proof. Unique factorisation makes $A$ multiplicative: an element of norm $mn$ with $\gcd(m,n)=1$ splits uniquely into parts of norms $m$ and $n$, because each Eisenstein prime has norm $p$, $q^2$ or $3$ for one ordinary prime. So it is enough to count prime powers. For $p\in P_1$ the elements of norm $p^{e}$ are exactly $\pi^i\bar\pi^{\,e-i}$ up to units, pairwise non-associate because $\pi$ and $\bar\pi$ are non-associate primes, so $A(p^{e})=e+1$. For $q\in P_2$ the only prime of $q$-power norm is $q$ itself, of norm $q^2$, so $A(q^{f})$ is $1$ for even $f$ and $0$ for odd $f$. For $3$ the only element of norm $3^{e}$ up to units is $\lambda^{e}$, so $A(3^{e})=1$ and the exponent of $3$ never affects the answer. Multiplying these factors gives the claim. ^[5] $\;\square$

Equivalently $r(n)=6\bigl(d_1(n)-d_2(n)\bigr)$, where $d_i(n)$ counts the divisors of $n$ congruent to $i$ modulo $3$. ^[4] ^[6] A brute-force lattice count confirms this identity for every $n\le600$, and the product rule agrees with $d_1-d_2$ for every $n\le2000$. The integers with $r(n)\gt0$ are the Loeschian numbers. ^[7]

The third checkpoint tests the formula at scale. Since $111111111=3^2\cdot37\cdot333667$ with $37$ and $333667$ both prime and both $\equiv1\pmod3$, we get $n=L^2/3=3^3\cdot37^2\cdot333667^2$, so

$$r(n)=6\cdot(2+1)\cdot(2+1)=54,$$

as the statement says. The factor $3^3$ contributes nothing.

5. From 450 to four exponent patterns

Because $450=6\cdot75$, the condition $\mathrm{B}(L)=450$ says that every $f_q$ is even and

$$\prod_{p\in P_1}(e_p+1)=75.$$

Only primes with $e_p\ge1$ appear, so each factor $e_p+1$ is at least $2$: we need the ways of writing $75$ as a product of integers greater than $1$, order ignored — its multiplicative partitions:

$$75=3\cdot5\cdot5=3\cdot25=5\cdot15=75,$$

four of them — the same count as for $20=2\cdot2\cdot5=2\cdot10=4\cdot5=20$, which has the same exponent pattern $p\,q^2$. ^[8] Subtracting $1$ from each factor leaves four exponent multisets,

$$\{4,4,2\},\quad\{24,2\},\quad\{14,4\},\quad\{74\},$$

so the part of $n$ built from primes in $P_1$ — we call it the core $c$ — has one of the shapes

$$p^4q^4r^2,\quad p^{24}q^2,\quad p^{14}q^4,\quad p^{74}$$

with distinct primes $p,q,r\in P_1$.

6. Core and neutral multiplier

What is left of $n$ must not change $r$: by Section 4 that allows any power of $3$ and even powers of primes in $P_2$. So the integers we are counting are exactly

$$n=3^{\alpha}u^2c,$$

where $c$ is a core and every prime factor of $u$ lies in $P_2$. The three parts live on disjoint sets of primes — $3$, $P_2$, $P_1$ — so $n$ determines $\alpha$, $u$ and $c$ uniquely, and every such product has $r=450$. Writing $t=3^{\alpha}u^2$ for the neutral multiplier, the answer counts the pairs $(c,t)$ with

$$ct\le N,\qquad N=\left\lfloor\frac{(5\cdot10^{11})^2}{3}\right\rfloor,$$

and $N=83333333333333333333333$.

7. Counting the neutral multipliers

Let $G(s)$ be the number of integers $u\le s$ all of whose prime factors lie in $P_2$, counting $u=1$. The admissible $u\le40$ are

$$\begin{aligned}&1,\ 2,\ 4,\ 5,\ 8,\ 10,\ 11,\ 16,\ 17,\\&20,\ 22,\ 23,\ 25,\ 29,\ 32,\ 34,\ 40,\end{aligned}$$

so $G(1),\dots,G(20)$ run

$$\begin{aligned}&1,2,2,3,4,4,4,5,5,6,\\&7,7,7,7,7,8,9,9,9,10.\end{aligned}$$

Splitting $t=3^{\alpha}u^2\le x$ by the exponent $\alpha$ and noting that the remaining condition is $u\le\sqrt{x/3^{\alpha}}$,

$$M(x)=\sum_{3^{\alpha}\le x}G\!\left(\left\lfloor\sqrt{x/3^{\alpha}}\right\rfloor\right)$$

is the number of neutral multipliers up to $x$. The neutral $t\le60$ are $1,3,4,9,12,16,25,27,36,48$, and $M(10)=4$, $M(30)=8$, $M(100)=14$, $M(1000)=38$; the formula was checked against a direct count for every $x\le200$.

The function $G$ is tabulated once with a sieve for the smallest prime factor $\operatorname{spf}(v)$. Set $\mathrm{good}(1)=1$; for $v\ge2$ put $p=\operatorname{spf}(v)$ and let $\mathrm{good}(v)=1$ when both $p\in P_2$ and $\mathrm{good}(v/p)=1$, and $\mathrm{good}(v)=0$ otherwise. Since $v/p\lt v$, one increasing pass fills the table. Prefix sums of $\mathrm{good}$ then turn every evaluation of $G$ into one array lookup, so each $M(x)$ costs only the $O(\log x)$ steps of its sum.

8. A complete worked example

Take the same question with $450$ replaced by $54$ and the limit by $500$, so $N=\lfloor500^2/3\rfloor=83333$. Now $54=6\cdot9$ and $9=3\cdot3=9$, giving exponent multisets $\{2,2\}$ and $\{8\}$. The family $p^8$ is empty, since already $7^8=5764801\gt83333$. The cores $p^2q^2\le83333$ with $p\lt q$ in $P_1$ are

$$7^2 13^2,\ 7^2 19^2,\ 7^2 31^2,\ 13^2 19^2,\ 7^2 37^2,$$

that is $8281$, $17689$, $47089$, $61009$ and $67081$. Their budgets $\lfloor N/c\rfloor$ are $10,4,1,1,1$, so the multiplier counts are $M(10)=4$, $M(4)=3$, $M(1)=1$, $M(1)=1$, $M(1)=1$: the multipliers are $1,3,4,9$ for the first core, $1,3,4$ for the second and $1$ for the rest. The total is $4+3+1+1+1=10$. Spelled out, the ten integers are

$$8281,\ 17689,\ 24843,\ 33124,\ 47089,$$

$$53067,\ 61009,\ 67081,\ 70756,\ 74529,$$

with $L=\sqrt{3n}$ running from $91\sqrt3\approx157.6$ to about $472.85$. A search over every $n\le83333$ finds exactly these ten.

9. Bounds that make the enumeration finite

If $a\lt b$ and $e\gt f$ then $a^e b^f\lt a^f b^e$, since the ratio of the two sides is $(b/a)^{e-f}\gt1$. So within a pattern the smallest core puts the largest exponents on the smallest eligible primes of $P_1$, giving the four pattern minima $7^4 13^4 19^2$, $7^{14}13^4$, $7^{24}13^2$ and $7^{74}$. The first is by far the smallest, so the smallest core of all is

$$c_{\min}=7^4 13^4 19^2=24755560921=157339^2.$$

So the smallest $L$ with $\mathrm{B}(L)=450$ is $157339\sqrt3\approx272519.14$, and the count is $0$ for any limit below it. At the full limit the four patterns are wildly uneven.

$p^{74}\le N$ has no solution — already $7^{74}$ has $63$ digits — so that family is empty.
$p^{24}\le N$ only for $p=7$; then $q^2\le\lfloor N/7^{24}\rfloor=434$ leaves $q\in\{13,19\}$: $2$ cores.
$p^{14}\le N$ for $p\in\{7,13,19,31,37,43\}$; that family has $58$ cores.
$\{4,4,2\}$ dominates: $p^4\le N$ for the $22164$ primes of $P_1$ up to $537241$.

One sieve bound covers every prime that can occur. A prime carrying the exponent $2$ is multiplied by $p^4q^4\ge7^4 13^4$ in the pattern $\{4,4,2\}$ and by $p^{24}\ge7^{24}$ in $\{24,2\}$, and the second cofactor is much the larger, so in both cases

$$p^2\le\frac{N}{7^4 13^4}\quad\Longrightarrow\quad p\le34859936.$$

Every other exponent is at least $4$, so those primes are at most $\lfloor N^{1/4}\rfloor=537284$, again below the bound. Sieving to $34859936$ yields $1069073$ primes of $P_1$. Likewise every core is at least $c_{\min}$, so $u^2\le N/c_{\min}$ and

$$u\le\left\lfloor\sqrt{N/c_{\min}}\right\rfloor=1834733,$$

which fixes the size of the neutral table.

10. Why the exponents must be assigned in every order

Inside a pattern the exponents belong to specific primes, and when they differ the assignment changes the integer. For $\{4,4,2\}$ on $7,13,19$ the three distinct arrangements give

$$7^4 13^4 19^2,\quad 7^4 13^2 19^4,\quad 7^2 13^4 19^4,$$

that is $24755560921$, $52880221849$ and $182382805969$, all three cores, with $L\approx272519.14$, $398297.21$ and $739694.81$. To meet each core once we keep the primes in increasing order and run over the distinct orderings of the multiset: one for $\{74\}$, two each for $\{24,2\}$ and $\{14,4\}$, three for $\{4,4,2\}$ — eight loops in all. The answer is then

$$\sum_{c}M\!\left(\left\lfloor N/c\right\rfloor\right)$$

over the cores $c\le N$ found by those loops.

11. Nothing is missed and nothing is counted twice

By Section 6 every admissible $n$ has exactly one core and one multiplier, and by Section 5 its core carries one of the four exponent multisets. Sorting that core's primes increasingly and reading their exponents in that order gives one ordered exponent sequence, which names exactly one of the eight loops, and the sorted primes are one index tuple inside it. For a fixed core, $M(\lfloor N/c\rfloor)$ counts the $t$ with $ct\le N$, each giving one $n$. So the sum counts each admissible $n$ once, and by Section 1 each $n$ gives one $L$.

How the Code Works

The C++ program fixes the targets as constants: kTargetL $=5\cdot10^{11}$, kTargetB $=450$ with kTargetDiv6 $=75$, and the two bounds of Section 9, kMinCore $=7^4 13^4 19^2$ and kCoreForExp2Bound $=7^4 13^4$. ^[9] Since $N\approx8.3\cdot10^{22}$ exceeds the unsigned $64$-bit range, the limit itself, the power tables, the shrinking budgets and the accumulated total are unsigned __int128, a compiler extension of GCC and Clang. ^[10] Smaller quantities keep narrower types: $64$ bits for the $L$ limit, the integer square roots and the multiplier counts, $32$ bits for the neutral prefix sums, int for the primes. Three helpers support it: pow_u128_limited returns limit + 1 instead of wrapping, so an overflowing power simply compares as too large; isqrt_u128 corrects a long double square root until it is exact; and to_string_u128 prints by repeated division.

n_limit_from_l_limit uses integer division to turn an $L$ limit into $N=\lfloor L^2/3\rfloor$. solve_for_n_limit returns $0$ at once when $N\lt$ kMinCore, then calls sieve_primes_mod_1_upto for the primes of $P_1$ below the bound of Section 9, tabulates pow2 and pow4 for all of them, and builds the three much shorter lists pow14_list, pow24_list and pow74_list, each stopping at the first prime whose power exceeds $N$. build_neutral_counter runs the sieve of Section 7 up to max_s and stores the prefix sums; NeutralCounter::count_multipliers evaluates $M(x)$.

Then come the eight blocks of Section 10. The $\{74\}$ block simply walks pow74_list; the other seven keep their distinct primes in increasing index order and use std::upper_bound on the sorted power tables to locate the end of an inner range instead of testing every prime. Budgets shrink by successive integer division, $b_0=N/p_i^{e}$ then $b_1=b_0/p_j^{e'}$, which is exact because $\lfloor\lfloor N/A\rfloor/B\rfloor=\lfloor N/(AB)\rfloor$.

run_validations runs first unless --skip-checks is given, and --limit-l= replaces the limit, which is how the smaller cases quoted below were produced. It checks count_representations_formula against count_representations_bruteforce, a direct double loop over the lattice, for every $n\le120$; the checkpoints $6$, $12$ and $54$; count_multipliers against is_multiplier_bruteforce, which factorises one $t$ at a time, at $x$ up to $10000$; that the two evaluation paths agree at $L=5\cdot10^7$; and $450/6=75$.

The Python solve mirrors this with math.isqrt, a bytearray sieve, an spf list and a prefix list, and needs no wide-integer type because Python integers are unbounded. Unlike the C++ entry point it hard-codes TARGET_L, takes no limit option and runs no packaged checks. Its count_mult guards against an index above max_s; the guard never fires, since every core is at least MIN_CORE and so $\lfloor\sqrt{\lfloor N/c\rfloor}\rfloor\le$ max_s. Where the C++ binary-searches an inner range end, Python scans linearly and breaks. The Java file does not reimplement the mathematics: it compiles solutionsCpp/Euler354.cpp with clang++ or g++ at -O2 into a cached binary, rebuilt only when the source is newer, runs it and reads the value from the last matching line of its output.

Complexity and Verification

Let $P=34859936$ and $S=1834733$ be the sieve bounds of Section 9. The prime sieve costs $O(P\log\log P)$ and yields $2138660$ primes, $1069073$ of them in $P_1$; the smallest-prime-factor sieve and the prefix sums cost $O(S\log\log S)$. Memory is $O(P+S)$: in C++ the sieve's byte flags take about $35$ MB and are released before the power tables are built, the retained prime list and the two $128$-bit power tables hold $36$ bytes per prime of $P_1$ (about $38.5$ MB), and the neutral table's three arrays of $S$ entries add $12$ bytes each (about $22$ MB), so some $60.5$ MB is live at the peak. Python stores the same data far less compactly.

The enumeration is small by comparison. The $\{74\}$ family is empty, $\{24,2\}$ gives $2$ cores and $\{14,4\}$ gives $58$; the three arrangements of $\{4,4,2\}$ give $3208094$, $52316$ and $14996$, so $C=3275466$ cores in all, each costing one evaluation of $M$: $O(\log N)$ divisions and integer square roots with $O(1)$ lookups. The outer and middle loops of the three-prime families take $66491$ and $20483937$ steps, $20205489$ of the latter in $p^4q^2r^4$ alone, where the middle loop is bounded only by $p_j^2\le b_0$ and almost every step finds an empty inner range; since $p_k\gt p_j$, a completion needs $p_j^6\lt b_0$, and testing that would remove nearly all of them. Writing $V$ for those steps and $Q$ for the retained primes of $P_1$, each step may cost a binary search, so the C++ time is

$$\begin{aligned}O(&P\log\log P+S\log\log S\\&+V\log Q+C\log N),\end{aligned}$$

and the sieves dominate at this limit. Python scans to the range endpoints instead, costing $44309750$ steps in the same loops. On the single machine used here the C++ finished in well under a second with its checks, the Python in about $7$ seconds.

Beyond the packaged checks, the following were verified here. The lattice brute force matches $6(d_1(n)-d_2(n))$ for $n\le600$, and the product rule of Section 4 matches $d_1-d_2$ for $n\le2000$. Direct enumeration in $\mathbb Z[\omega]$ agreed with Section 3: $\bar\pi$ is never an associate of $\pi$ for $p\in P_1$ up to $163$, and nothing has norm $q$ for $q\in P_2$ up to $101$. The counts $A(p^{e})=e+1$, $A(3^{e})=1$ and $A(q^{f})=1$ or $0$ were checked for $p=7,13$ with $e\le3$, $3^{e}$ with $e\le6$, and $q=2,5,11$ with $f\le4$ — finite checks of the general argument in Section 4. The formula for $M$ matches a direct count for every $x\le200$.

The machinery was then tested against a complete search over every $n\le\lfloor L^2/3\rfloor$, feasible when the target is small. Both give $3988$ for $\mathrm{B}=12,L\le300$; $101$ for $\mathrm{B}=18,L\le300$; $3399$ for $\mathrm{B}=24,L\le400$; $13$ for $\mathrm{B}=30,L\le600$; $33$ for $\mathrm{B}=54,L\le1000$; $3738$ for $\mathrm{B}=36,L\le1500$; and $0$ for $\mathrm{B}=84,L\le2000$.

Finally, a separate program that assigns an exponent to each prime in increasing order, with no core-and-multiplier split and no neutral table, reproduced the published counts at $\mathrm{B}=450$: $0$ at $272519$ and $1$ at $272520$, which locates the smallest $L$ exactly where Section 9 predicts, then $1$ at $3\cdot10^5$, $18$ at $10^6$, $112$ at $3\cdot10^6$, $597$ at $10^7$, $2297$ at $3\cdot10^7$, $9135$ at $10^8$, $107946$ at $10^9$, $1146952$ at $10^{10}$ and $11640790$ at $10^{11}$. At $L\le10^6$ all eighteen integers $n$ were listed and factorised independently to confirm $r(n)=450$. None of these cross-checks is part of the published solutions.

One pitfall is worth recording. Because $4n=(2a+b)^2+3b^2$, a solution of $a^2+ab+b^2=n$ has $|b|\le2\sqrt{n/3}$, and $|a|$ likewise by symmetry. The radius $\lfloor\sqrt n\rfloor+2$ used by count_representations_bruteforce is therefore only just enough for $n\le120$, the range it is called on, and the two bounds cross just above it: at $n=192$ that loop finds $2$ solutions where there are $6$. A complete search for any $n$ can use radius $\lfloor2\sqrt{n/3}\rfloor$.

Footnotes and References

The references below supply the background facts used above. The reduction to cores and multipliers, the four exponent patterns and the bounds of Section 9 are derived here.

Project Euler 354 — Distances in a Bee's Honeycomb. The official statement defines B(L) as the number of cells of a honeycomb of unit regular hexagons whose centre is at distance L from the queen's cell, gives the checkpoints B(√3) = 6, B(√21) = 12 and B(111111111) = 54, and asks for the number of L ≤ 5·10¹¹ with B(L) = 450.
Wikipedia — Hexagon. For a regular hexagon of side t the article gives the inradius, or apothem, as (√3/2)·t and the circumradius as t itself. The apothem puts the centres of two edge-adjacent cells √3 apart in Section 1.
Wikipedia — Hexagonal tiling. The regular hexagonal tiling has Schläfli symbol {6,3} and is one of three regular tilings of the plane; its dual is the triangular tiling, since joining neighbouring hexagon centres produces equilateral triangles — the lattice of Section 1.
OEIS A004016. The theta series of the planar hexagonal lattice A₂, with terms 1, 6, 0, 6, 6, 0, 0, 12, 0, 6, 0, 0, 6, 12, … It equals 6·A002324(n) for n > 0, and states the multiplicative form used in Section 4.
Wikipedia — Eisenstein integer. The ring ℤ[ω] with ω = e^(2πi/3) forms a triangular lattice; its norm is N(x + yω) = x² − xy + y²; it is a Euclidean domain, so unique factorisation holds; its six units ±1, ±ω, ±ω² are its elements of norm 1; primes ≡ 2 mod 3 stay prime, primes ≡ 1 mod 3 factor as p = (a + bω)((a − b) − bω), and 3 = −(1 + 2ω)². Sections 3 and 4 rest on this.
OEIS A002324. The number of divisors of n congruent to 1 mod 3 minus the number congruent to 2 mod 3. It is multiplicative, with a(3^e) = 1, a(p^e) = e + 1 for p ≡ 1 mod 3 and a(p^e) = (1 + (−1)^e)/2 for p ≡ 2 mod 3 — the divisor form of Section 4.
Wikipedia — Loeschian number. A Loeschian number has the form x² + xy + y². The article characterises them as the integers in whose factorisation every prime congruent to 2 mod 3 has an even power — the condition for r(n) > 0 in Section 4 — and notes that they are the Eisenstein norms.
Wikipedia — Multiplicative partition. A multiplicative partition of n writes n as a product of integers greater than 1, order ignored and n itself counted. The article lists the four for 20 as 2 × 2 × 5, 2 × 10, 4 × 5 and 20, the pattern Section 5 uses for 75.
C++ — Euler354.cpp. The linked immutable revision contains count_representations_formula, build_neutral_counter, solve_for_n_limit with its eight pattern blocks, and run_validations. It is the implementation record for this article.
GCC manual — 128-bit integers. The page documents __int128 and unsigned __int128 as an extension available on targets with a wide enough integer mode. The solution uses the unsigned form because N exceeds the 64-bit range.

Problem 354 source code

C++

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

namespace {

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

constexpr u64 kTargetL = 500'000'000'000ULL;
constexpr u64 kTargetB = 450ULL;
constexpr u64 kTargetDiv6 = kTargetB / 6ULL;  // 75
constexpr u128 kMinCore = static_cast<u128>(2'401ULL) * 28'561ULL * 361ULL;       // 7^4 * 13^4 * 19^2
constexpr u128 kCoreForExp2Bound = static_cast<u128>(2'401ULL) * 28'561ULL;        // 7^4 * 13^4

struct Options {
    u64 limit_l = kTargetL;
    bool run_checks = true;
    bool allow_multithreading = 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) != 0U) 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(const int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

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

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

        u64 parsed_l = 0ULL;
        if (parse_u64_after_prefix(arg, "--limit-l=", parsed_l)) {
            options.limit_l = parsed_l;
            continue;
        }

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

    if (options.limit_l == 0ULL) {
        std::cerr << "L must be positive.\n";
        return false;
    }

    return true;
}

unsigned choose_thread_count(const bool allow_multithreading,
                             const unsigned requested_threads,
                             const std::size_t workload_units) {
    if (!allow_multithreading || workload_units < 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_units)));
}

u128 from_u64(const u64 x) {
    return static_cast<u128>(x);
}

std::string to_string_u128(u128 value) {
    if (value == 0) return "0";

    std::string out;
    while (value > 0) {
        const u32 digit = static_cast<u32>(value % 10);
        out.push_back(static_cast<char>('0' + digit));
        value /= 10;
    }
    std::reverse(out.begin(), out.end());
    return out;
}

u128 pow_u128_limited(const u64 base, const int exponent, const u128 limit) {
    u128 out = 1;
    const u128 b = from_u64(base);
    for (int i = 0; i < exponent; ++i) {
        if (out > limit / b) return limit + 1;
        out *= b;
    }
    return out;
}

u64 isqrt_u128(const u128 value) {
    if (value == 0) return 0ULL;

    long double approx = std::sqrt(static_cast<long double>(value));
    u64 root = static_cast<u64>(approx);
    while (from_u64(root + 1ULL) * from_u64(root + 1ULL) <= value) ++root;
    while (from_u64(root) * from_u64(root) > value) --root;
    return root;
}

u128 square_u64(const u64 x) {
    return from_u64(x) * from_u64(x);
}

std::vector<int> sieve_primes_upto(const int limit) {
    std::vector<int> primes;
    if (limit < 2) return primes;

    std::vector<std::uint8_t> is_composite(static_cast<std::size_t>(limit + 1), 0U);
    primes.reserve(static_cast<std::size_t>(limit / 10));

    for (int i = 2; i <= limit; ++i) {
        if (is_composite[static_cast<std::size_t>(i)] == 0U) {
            primes.push_back(i);
            if (from_u64(static_cast<u64>(i)) * from_u64(static_cast<u64>(i)) <= from_u64(static_cast<u64>(limit))) {
                for (u64 j = from_u64(static_cast<u64>(i)) * from_u64(static_cast<u64>(i));
                     j <= static_cast<u64>(limit);
                     j += static_cast<u64>(i)) {
                    is_composite[static_cast<std::size_t>(j)] = 1U;
                }
            }
        }
    }

    return primes;
}

std::vector<int> sieve_primes_mod_1_upto(const int limit) {
    const std::vector<int> all_primes = sieve_primes_upto(limit);
    std::vector<int> mod1;
    mod1.reserve(all_primes.size() / 2U + 8U);
    for (const int p : all_primes) {
        if (p % 3 == 1) mod1.push_back(p);
    }
    return mod1;
}

u64 count_representations_formula(u64 n, const std::vector<int>& factor_primes) {
    u64 product = 1ULL;

    for (const int prime : factor_primes) {
        const u64 p = static_cast<u64>(prime);
        if (p * p > n) break;
        if (n % p != 0ULL) continue;

        int exponent = 0;
        while (n % p == 0ULL) {
            n /= p;
            ++exponent;
        }

        if (p % 3ULL == 1ULL) {
            product *= static_cast<u64>(exponent + 1);
        } else if (p % 3ULL == 2ULL) {
            if ((exponent & 1) != 0) return 0ULL;
        }
    }

    if (n > 1ULL) {
        if (n % 3ULL == 1ULL) {
            product *= 2ULL;
        } else if (n % 3ULL == 2ULL) {
            return 0ULL;
        }
    }

    return 6ULL * product;
}

u64 count_representations_bruteforce(const int target_n) {
    const int radius = static_cast<int>(std::sqrt(static_cast<double>(target_n))) + 2;
    u64 count = 0ULL;
    for (int a = -radius; a <= radius; ++a) {
        for (int b = -radius; b <= radius; ++b) {
            const int value = a * a + a * b + b * b;
            if (value == target_n) ++count;
        }
    }
    return count;
}

bool is_multiplier_bruteforce(u64 value, const std::vector<int>& factor_primes) {
    u64 n = value;
    for (const int prime : factor_primes) {
        const u64 p = static_cast<u64>(prime);
        if (p * p > n) break;
        if (n % p != 0ULL) continue;

        int exponent = 0;
        while (n % p == 0ULL) {
            n /= p;
            ++exponent;
        }

        if (p == 3ULL) continue;
        if (p % 3ULL == 1ULL) return false;
        if ((exponent & 1) != 0) return false;
    }

    if (n > 1ULL) {
        if (n == 3ULL) return true;
        if (n % 3ULL == 1ULL) return false;
        return false;  // n % 3 == 2 with odd exponent 1.
    }

    return true;
}

struct NeutralCounter {
    std::vector<u32> prefix;

    u64 count_multipliers(const u128 x) const {
        u64 total = 0ULL;
        u128 power3 = 1;

        while (power3 <= x) {
            const u64 y = isqrt_u128(x / power3);
            total += prefix[static_cast<std::size_t>(y)];
            if (power3 > x / 3U) break;
            power3 *= 3U;
        }

        return total;
    }
};

NeutralCounter build_neutral_counter(const u64 max_s) {
    NeutralCounter out;
    out.prefix.assign(static_cast<std::size_t>(max_s + 1U), 0U);
    if (max_s == 0ULL) return out;

    std::vector<int> spf(static_cast<std::size_t>(max_s + 1U), 0);
    for (u64 i = 2ULL; i <= max_s; ++i) {
        if (spf[static_cast<std::size_t>(i)] == 0) {
            spf[static_cast<std::size_t>(i)] = static_cast<int>(i);
            if (i <= max_s / i) {
                for (u64 j = i * i; j <= max_s; j += i) {
                    if (spf[static_cast<std::size_t>(j)] == 0) {
                        spf[static_cast<std::size_t>(j)] = static_cast<int>(i);
                    }
                }
            }
        }
    }

    std::vector<u32> good(static_cast<std::size_t>(max_s + 1U), 0U);
    good[1] = 1U;

    for (u64 n = 2ULL; n <= max_s; ++n) {
        const int p = spf[static_cast<std::size_t>(n)];
        const u64 m = n / static_cast<u64>(p);
        good[static_cast<std::size_t>(n)] =
            (p % 3 == 2 && good[static_cast<std::size_t>(m)] == 1U) ? 1U : 0U;
    }

    for (u64 n = 1ULL; n <= max_s; ++n) {
        out.prefix[static_cast<std::size_t>(n)] =
            out.prefix[static_cast<std::size_t>(n - 1U)] + good[static_cast<std::size_t>(n)];
    }

    return out;
}

template <typename Fn>
u128 parallel_sum_indices(const std::size_t i_end,
                          const bool allow_multithreading,
                          const unsigned requested_threads,
                          Fn&& fn) {
    if (i_end == 0U) return 0;

    const unsigned threads =
        choose_thread_count(allow_multithreading, requested_threads, i_end);
    std::vector<u128> partial(threads, 0);
    std::atomic<std::size_t> next_i{0U};

    const auto worker = [&](const unsigned tid) {
        u128 local = 0;
        while (true) {
            const std::size_t i = next_i.fetch_add(1U, std::memory_order_relaxed);
            if (i >= i_end) break;
            local += fn(i);
        }
        partial[tid] = local;
    };

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

    u128 total = 0;
    for (const u128 value : partial) total += value;
    return total;
}

u128 solve_for_n_limit(const u128 n_limit,
                       const bool allow_multithreading,
                       const unsigned requested_threads) {
    if (n_limit < kMinCore) return 0;
    const int prime_limit = static_cast<int>(isqrt_u128(n_limit / kCoreForExp2Bound));

    const std::vector<int> primes_mod1 = sieve_primes_mod_1_upto(prime_limit);
    const std::size_t ps = primes_mod1.size();

    std::vector<u128> pow2(ps, n_limit + 1);
    std::vector<u128> pow4(ps, n_limit + 1);
    for (std::size_t i = 0; i < ps; ++i) {
        const u64 p = static_cast<u64>(primes_mod1[i]);
        pow2[i] = pow_u128_limited(p, 2, n_limit);
        pow4[i] = pow_u128_limited(p, 4, n_limit);
    }

    std::vector<std::pair<std::size_t, u128>> pow14_list;
    std::vector<std::pair<std::size_t, u128>> pow24_list;
    std::vector<std::pair<std::size_t, u128>> pow74_list;

    for (std::size_t i = 0; i < ps; ++i) {
        const u64 p = static_cast<u64>(primes_mod1[i]);
        const u128 p14 = pow_u128_limited(p, 14, n_limit);
        if (p14 <= n_limit) {
            pow14_list.push_back({i, p14});
        } else {
            break;
        }
    }
    for (std::size_t i = 0; i < ps; ++i) {
        const u64 p = static_cast<u64>(primes_mod1[i]);
        const u128 p24 = pow_u128_limited(p, 24, n_limit);
        if (p24 <= n_limit) {
            pow24_list.push_back({i, p24});
        } else {
            break;
        }
    }
    for (std::size_t i = 0; i < ps; ++i) {
        const u64 p = static_cast<u64>(primes_mod1[i]);
        const u128 p74 = pow_u128_limited(p, 74, n_limit);
        if (p74 <= n_limit) {
            pow74_list.push_back({i, p74});
        } else {
            break;
        }
    }

    const u64 max_s = isqrt_u128(n_limit / kMinCore);
    const NeutralCounter neutral = build_neutral_counter(max_s);

    const std::size_t i_end4 = static_cast<std::size_t>(
        std::upper_bound(pow4.begin(), pow4.end(), n_limit) - pow4.begin());

    u128 answer = 0;

    // Pattern [74].
    for (const auto& entry : pow74_list) {
        answer += neutral.count_multipliers(n_limit / entry.second);
    }

    // Pattern [24, 2].
    for (const auto& entry : pow24_list) {
        const std::size_t i = entry.first;
        const u128 b0 = n_limit / entry.second;
        const std::size_t j_begin = i + 1U;
        if (j_begin >= ps) continue;
        const std::size_t j_end = static_cast<std::size_t>(
            std::upper_bound(pow2.begin() + static_cast<std::ptrdiff_t>(j_begin), pow2.end(), b0) -
            pow2.begin());
        for (std::size_t j = j_begin; j < j_end; ++j) {
            answer += neutral.count_multipliers(b0 / pow2[j]);
        }
    }

    // Pattern [2, 24].
    for (const auto& entry : pow24_list) {
        const std::size_t j = entry.first;
        if (j == 0U) continue;
        const u128 b0 = n_limit / entry.second;
        const std::size_t i_end = static_cast<std::size_t>(
            std::upper_bound(pow2.begin(), pow2.begin() + static_cast<std::ptrdiff_t>(j), b0) -
            pow2.begin());
        for (std::size_t i = 0; i < i_end; ++i) {
            answer += neutral.count_multipliers(b0 / pow2[i]);
        }
    }

    // Pattern [14, 4].
    for (const auto& entry : pow14_list) {
        const std::size_t i = entry.first;
        const u128 b0 = n_limit / entry.second;
        const std::size_t j_begin = i + 1U;
        if (j_begin >= ps) continue;
        const std::size_t j_end = static_cast<std::size_t>(
            std::upper_bound(pow4.begin() + static_cast<std::ptrdiff_t>(j_begin), pow4.end(), b0) -
            pow4.begin());
        for (std::size_t j = j_begin; j < j_end; ++j) {
            answer += neutral.count_multipliers(b0 / pow4[j]);
        }
    }

    // Pattern [4, 14].
    for (const auto& entry : pow14_list) {
        const std::size_t j = entry.first;
        if (j == 0U) continue;
        const u128 b0 = n_limit / entry.second;
        const std::size_t i_end = static_cast<std::size_t>(
            std::upper_bound(pow4.begin(), pow4.begin() + static_cast<std::ptrdiff_t>(j), b0) -
            pow4.begin());
        for (std::size_t i = 0; i < i_end; ++i) {
            answer += neutral.count_multipliers(b0 / pow4[i]);
        }
    }

    // Pattern [4, 4, 2].
    answer += parallel_sum_indices(
        i_end4, allow_multithreading, requested_threads, [&](const std::size_t i) -> u128 {
            const u128 b0 = n_limit / pow4[i];
            const std::size_t j_begin = i + 1U;
            if (j_begin >= ps) return 0;
            const std::size_t j_end = static_cast<std::size_t>(
                std::upper_bound(pow4.begin() + static_cast<std::ptrdiff_t>(j_begin), pow4.end(), b0) -
                pow4.begin());

            u128 local = 0;
            for (std::size_t j = j_begin; j < j_end; ++j) {
                const u128 b1 = b0 / pow4[j];
                const std::size_t k_begin = j + 1U;
                if (k_begin >= ps) continue;
                const std::size_t k_end = static_cast<std::size_t>(
                    std::upper_bound(pow2.begin() + static_cast<std::ptrdiff_t>(k_begin),
                                     pow2.end(),
                                     b1) -
                    pow2.begin());
                for (std::size_t k = k_begin; k < k_end; ++k) {
                    local += neutral.count_multipliers(b1 / pow2[k]);
                }
            }
            return local;
        });

    // Pattern [4, 2, 4].
    answer += parallel_sum_indices(
        i_end4, allow_multithreading, requested_threads, [&](const std::size_t i) -> u128 {
            const u128 b0 = n_limit / pow4[i];
            const std::size_t j_begin = i + 1U;
            if (j_begin >= ps) return 0;
            const std::size_t j_end = static_cast<std::size_t>(
                std::upper_bound(pow2.begin() + static_cast<std::ptrdiff_t>(j_begin), pow2.end(), b0) -
                pow2.begin());

            u128 local = 0;
            for (std::size_t j = j_begin; j < j_end; ++j) {
                const u128 b1 = b0 / pow2[j];
                const std::size_t k_begin = j + 1U;
                if (k_begin >= ps) continue;
                const std::size_t k_end = static_cast<std::size_t>(
                    std::upper_bound(pow4.begin() + static_cast<std::ptrdiff_t>(k_begin),
                                     pow4.end(),
                                     b1) -
                    pow4.begin());
                for (std::size_t k = k_begin; k < k_end; ++k) {
                    local += neutral.count_multipliers(b1 / pow4[k]);
                }
            }
            return local;
        });

    // Pattern [2, 4, 4].
    const std::size_t i_end_for_244 = (i_end4 > 0U) ? (i_end4 - 1U) : 0U;
    answer += parallel_sum_indices(
        i_end_for_244, allow_multithreading, requested_threads, [&](const std::size_t i) -> u128 {
            const u128 b0 = n_limit / pow2[i];
            const std::size_t j_begin = i + 1U;
            if (j_begin >= ps) return 0;
            const std::size_t j_end = static_cast<std::size_t>(
                std::upper_bound(pow4.begin() + static_cast<std::ptrdiff_t>(j_begin), pow4.end(), b0) -
                pow4.begin());

            u128 local = 0;
            for (std::size_t j = j_begin; j < j_end; ++j) {
                const u128 b1 = b0 / pow4[j];
                const std::size_t k_begin = j + 1U;
                if (k_begin >= ps) continue;
                const std::size_t k_end = static_cast<std::size_t>(
                    std::upper_bound(pow4.begin() + static_cast<std::ptrdiff_t>(k_begin),
                                     pow4.end(),
                                     b1) -
                    pow4.begin());
                for (std::size_t k = k_begin; k < k_end; ++k) {
                    local += neutral.count_multipliers(b1 / pow4[k]);
                }
            }
            return local;
        });

    return answer;
}

u128 n_limit_from_l_limit(const u64 l_limit) {
    return square_u64(l_limit) / 3U;
}

bool run_validations(const Options& options) {
    const std::vector<int> small_primes = sieve_primes_upto(5'000'000);

    // Checkpoint 1: formula for r(n) against brute lattice counting.
    for (int n = 1; n <= 120; ++n) {
        const u64 brute = count_representations_bruteforce(n);
        const u64 fast = count_representations_formula(static_cast<u64>(n), small_primes);
        if (brute != fast) {
            std::cerr << "Validation failed: representation mismatch at n=" << n
                      << " (brute=" << brute << ", fast=" << fast << ")\n";
            return false;
        }
    }

    // Checkpoint 2: statement examples.
    if (count_representations_formula(1ULL, small_primes) != 6ULL) {
        std::cerr << "Validation failed: B(sqrt(3)) != 6\n";
        return false;
    }
    if (count_representations_formula(7ULL, small_primes) != 12ULL) {
        std::cerr << "Validation failed: B(sqrt(21)) != 12\n";
        return false;
    }

    const u64 l_example = 111'111'111ULL;
    const u64 n_example = static_cast<u64>(n_limit_from_l_limit(l_example));
    if (count_representations_formula(n_example, small_primes) != 54ULL) {
        std::cerr << "Validation failed: B(111111111) != 54\n";
        return false;
    }

    // Checkpoint 3: neutral multiplier counter against brute force.
    {
        const u64 max_s = 2'000ULL;
        const NeutralCounter neutral = build_neutral_counter(max_s);
        const std::vector<u64> tests = {1ULL, 3ULL, 10ULL, 100ULL, 1'000ULL, 10'000ULL};

        for (const u64 x : tests) {
            u64 brute = 0ULL;
            for (u64 t = 1ULL; t <= x; ++t) {
                if (is_multiplier_bruteforce(t, small_primes)) ++brute;
            }

            const u64 fast = neutral.count_multipliers(from_u64(x));
            if (brute != fast) {
                std::cerr << "Validation failed: neutral multiplier mismatch at x=" << x
                          << " (brute=" << brute << ", fast=" << fast << ")\n";
                return false;
            }
        }
    }

    // Checkpoint 4: single-thread vs multi-thread consistency on a reduced limit.
    if (options.allow_multithreading) {
        unsigned hw = std::thread::hardware_concurrency();
        if (hw == 0U) hw = 1U;
        if (hw > 1U) {
            const u64 l_test = 50'000'000ULL;
            const u128 n_test = n_limit_from_l_limit(l_test);
            const u128 single = solve_for_n_limit(n_test, false, 1U);
            const u128 multi = solve_for_n_limit(n_test, true, options.requested_threads);
            if (single != multi) {
                std::cerr << "Validation failed: thread consistency mismatch at L=" << l_test
                          << " (single=" << to_string_u128(single)
                          << ", multi=" << to_string_u128(multi) << ")\n";
                return false;
            }
        }
    }

    if (kTargetDiv6 != 75ULL) {
        std::cerr << "Validation failed: target setup mismatch.\n";
        return false;
    }

    return true;
}

}  // namespace

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

    if (options.run_checks) {
        if (!run_validations(options)) return 1;
    }

    const u128 n_limit = n_limit_from_l_limit(options.limit_l);
    const u128 answer = solve_for_n_limit(n_limit,
                                          options.allow_multithreading,
                                          options.requested_threads);
    std::cout << to_string_u128(answer) << '\n';
    return 0;
}

Python

import math

def solve():
    TARGET_L = 500_000_000_000

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

    def sieve_primes_mod1(limit):
        is_comp = bytearray(limit + 1)
        primes = []
        for i in range(2, limit + 1):
            if not is_comp[i]:
                primes.append(i)
                if i <= limit // i:
                    for j in range(i*i, limit+1, i):
                        is_comp[j] = 1
        return [p for p in primes if p % 3 == 1]

    n_limit = TARGET_L * TARGET_L // 3
    MIN_CORE = 2401 * 28561 * 361  # 7^4 * 13^4 * 19^2
    CORE_BOUND = 2401 * 28561      # 7^4 * 13^4

    if n_limit < MIN_CORE:
        return '0'

    prime_limit = isqrt(n_limit // CORE_BOUND)
    primes_m1 = sieve_primes_mod1(prime_limit)
    ps = len(primes_m1)

    def pp(base, exp):
        r = 1
        for _ in range(exp):
            if r > n_limit // base:
                return n_limit + 1
            r *= base
        return r

    pow2 = [pp(p, 2) for p in primes_m1]
    pow4 = [pp(p, 4) for p in primes_m1]
    pow14 = [(i, pp(p, 14)) for i, p in enumerate(primes_m1) if pp(p, 14) <= n_limit]
    pow24 = [(i, pp(p, 24)) for i, p in enumerate(primes_m1) if pp(p, 24) <= n_limit]
    pow74 = [(i, pp(p, 74)) for i, p in enumerate(primes_m1) if pp(p, 74) <= n_limit]

    max_s = isqrt(n_limit // MIN_CORE)

    # Build neutral multiplier prefix sums
    spf = list(range(max_s + 1))
    for i in range(2, max_s + 1):
        if spf[i] == i:
            if i <= max_s // i:
                for j in range(i*i, max_s+1, i):
                    if spf[j] == j:
                        spf[j] = i

    good = [0] * (max_s + 1)
    good[1] = 1
    for n in range(2, max_s + 1):
        p = spf[n]
        m = n // p
        good[n] = 1 if (p % 3 == 2 and good[m] == 1) else 0

    prefix = [0] * (max_s + 1)
    for n in range(1, max_s + 1):
        prefix[n] = prefix[n-1] + good[n]

    def count_mult(x):
        total = 0
        p3 = 1
        while p3 <= x:
            y = isqrt(x // p3)
            if y <= max_s:
                total += prefix[y]
            if p3 > x // 3:
                break
            p3 *= 3
        return total

    import bisect
    answer = 0

    # Pattern [74]
    for i, v in pow74:
        answer += count_mult(n_limit // v)

    # Pattern [24, 2]
    for i, v24 in pow24:
        b0 = n_limit // v24
        for j in range(i+1, ps):
            if pow2[j] > b0: break
            answer += count_mult(b0 // pow2[j])

    # Pattern [2, 24]
    for j, v24 in pow24:
        if j == 0: continue
        b0 = n_limit // v24
        for i in range(j):
            if pow2[i] > b0: break
            answer += count_mult(b0 // pow2[i])

    # Pattern [14, 4]
    for i, v14 in pow14:
        b0 = n_limit // v14
        for j in range(i+1, ps):
            if pow4[j] > b0: break
            answer += count_mult(b0 // pow4[j])

    # Pattern [4, 14]
    for j, v14 in pow14:
        if j == 0: continue
        b0 = n_limit // v14
        for i in range(j):
            if pow4[i] > b0: break
            answer += count_mult(b0 // pow4[i])

    # Pattern [4, 4, 2]
    i_end4 = 0
    for i in range(ps):
        if pow4[i] <= n_limit:
            i_end4 = i + 1
        else:
            break
    for i in range(i_end4):
        b0 = n_limit // pow4[i]
        for j in range(i+1, ps):
            if pow4[j] > b0: break
            b1 = b0 // pow4[j]
            for k in range(j+1, ps):
                if pow2[k] > b1: break
                answer += count_mult(b1 // pow2[k])

    # Pattern [4, 2, 4]
    for i in range(i_end4):
        b0 = n_limit // pow4[i]
        for j in range(i+1, ps):
            if pow2[j] > b0: break
            b1 = b0 // pow2[j]
            for k in range(j+1, ps):
                if pow4[k] > b1: break
                answer += count_mult(b1 // pow4[k])

    # Pattern [2, 4, 4]
    for i in range(i_end4 - 1 if i_end4 > 0 else 0):
        b0 = n_limit // pow2[i]
        for j in range(i+1, ps):
            if pow4[j] > b0: break
            b1 = b0 // pow4[j]
            for k in range(j+1, ps):
                if pow4[k] > b1: break
                answer += count_mult(b1 // pow4[k])

    return str(answer)

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

Java

import java.nio.file.*;
import java.util.*;
import java.util.regex.*;

public class Euler354 {
    private static final Pattern ANSWER_RE = Pattern.compile("answer\\s*:\\s*(.+)$", Pattern.CASE_INSENSITIVE);
    private static final Pattern EQUAL_RE = Pattern.compile("=\\s*(.+)$");

    private static String parseOutput(String stdout) {
        String[] lines = stdout.split("\\R");
        List<String> nonEmpty = new ArrayList<>();
        for (String line : lines) {
            String t = line.trim();
            if (!t.isEmpty()) {
                nonEmpty.add(t);
            }
        }
        if (nonEmpty.isEmpty()) {
            return "";
        }

        List<String> answers = new ArrayList<>();
        List<String> equals = new ArrayList<>();
        for (String line : nonEmpty) {
            Matcher m1 = ANSWER_RE.matcher(line);
            if (m1.find()) {
                answers.add(m1.group(1).trim());
            }
            Matcher m2 = EQUAL_RE.matcher(line);
            if (m2.find()) {
                equals.add(m2.group(1).trim());
            }
        }

        if (!answers.isEmpty()) {
            return answers.get(answers.size() - 1);
        }
        if (!equals.isEmpty()) {
            return equals.get(equals.size() - 1);
        }
        return nonEmpty.get(nonEmpty.size() - 1);
    }

    private static String pickCompiler() throws Exception {
        for (String compiler : List.of("clang++", "g++")) {
            Process probe = new ProcessBuilder("bash", "-lc", "command -v " + compiler)
                    .redirectErrorStream(true)
                    .start();
            String out = new String(probe.getInputStream().readAllBytes());
            int rc = probe.waitFor();
            if (rc == 0 && !out.trim().isEmpty()) {
                return compiler;
            }
        }
        throw new RuntimeException("No C++ compiler found (clang++/g++).");
    }

    private static Path cppSource(Path root) {
        return root.resolve("solutionsCpp").resolve("Euler354.cpp");
    }

    private static boolean shouldSkipCheckpoints(Path root) {
        Path src = cppSource(root);
        try {
            String text = Files.readString(src);
            return text.contains("--skip-checkpoints");
        } catch (Exception ex) {
            return false;
        }
    }

    private static Path ensureBridgeBinary() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = root.resolve("solutionsCpp").resolve(".euler354_java_bridge");

        boolean rebuild = Files.notExists(bin)
                || Files.getLastModifiedTime(src).compareTo(Files.getLastModifiedTime(bin)) > 0;

        if (rebuild) {
            String compiler = pickCompiler();
            Process compile = new ProcessBuilder(
                    compiler,
                    "-std=c++17",
                    "-O2",
                    src.toString(),
                    "-o",
                    bin.toString())
                    .inheritIO()
                    .start();
            if (compile.waitFor() != 0) {
                throw new RuntimeException("Failed to compile Euler354 C++ bridge.");
            }
        }

        return bin;
    }

    private static String runBridge(Path bin, Path root, Path srcDir) throws Exception {
        List<String> cmd = new ArrayList<>();
        cmd.add(bin.toString());
        if (shouldSkipCheckpoints(root)) {
            cmd.add("--skip-checkpoints");
        }

        Process first = new ProcessBuilder(cmd)
                .directory(root.toFile())
                .redirectErrorStream(true)
                .start();
        String out = new String(first.getInputStream().readAllBytes());
        int rc = first.waitFor();
        if (rc == 0) {
            return out;
        }

        Process second = new ProcessBuilder(cmd)
                .directory(srcDir.toFile())
                .redirectErrorStream(true)
                .start();
        String out2 = new String(second.getInputStream().readAllBytes());
        int rc2 = second.waitFor();
        if (rc2 == 0) {
            return out2;
        }

        throw new RuntimeException("Euler354 C++ bridge failed.\n" + out + "\n" + out2);
    }

    private static String solveViaCppBridge() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = ensureBridgeBinary();
        String out = runBridge(bin, root, src.getParent());
        String parsed = parseOutput(out);
        if (parsed.isEmpty()) {
            throw new RuntimeException("Euler354 C++ bridge produced empty output.");
        }
        return parsed;
    }

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