Detailed mathematical approach

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.

Step 2: Representation Formula for \(a^2+ab+b^2\)

Write the prime factorization of \(n\) as

$n=3^{e_0}\prod_{p\equiv 1 \pmod{3}}p^{e_p}\prod_{q\equiv 2 \pmod{3}}q^{f_q}.$

Then the representation count is determined by the congruence class of each prime:

if some exponent \(f_q\) is odd, then

$r(n)=0;$

otherwise,

$r(n)=6\prod_{p\equiv 1 \pmod{3}}(e_p+1).$

The reason is structural. Primes \(p\equiv 1 \pmod{3}\) split in the Eisenstein integers and contribute a factor \(e_p+1\). Primes \(q\equiv 2 \pmod{3}\) are inert, so they can only appear with even exponent if the norm is to remain representable. The prime \(3\) is ramified, but it does not change the multiplicity factor beyond the universal leading constant \(6\).

Step 3: Turn \(r(n)=450\) into Exponent Patterns

Since

$450=6\cdot 75,$

we must have

$\prod_{p\equiv 1 \pmod{3}}(e_p+1)=75.$

Now factor \(75\) into integers greater than \(1\):

$75,\qquad 25\cdot 3,\qquad 15\cdot 5,\qquad 5\cdot 5\cdot 3.$

Subtracting \(1\) from each factor gives the only possible exponent multisets on primes \(p\equiv 1 \pmod{3}\):

$[74],\qquad [24,2],\qquad [14,4],\qquad [4,4,2].$

So every admissible core is of one of the forms

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

with distinct primes \(p,q,r\equiv 1 \pmod{3}\). Whenever the exponents are not all equal, their assignment to distinct primes matters. This explains the separate code paths for \([24,2]\) and \([2,24]\), for \([14,4]\) and \([4,14]\), and for the three ordered versions of \([4,4,2]\).

Step 4: Core and Neutral Multiplier Decomposition

Every valid integer \(n\) can be written uniquely as

$n=c\cdot t,$

where \(c\) is a core carrying one of the exponent patterns above on primes \(p\equiv 1 \pmod{3}\), and \(t\) is a multiplier that does not change the value \(r(n)=450\).

The only neutral factors are powers of \(3\) and even powers of primes \(q\equiv 2 \pmod{3}\). Therefore

$t=3^a u^2,$

where every prime divisor of \(u\) satisfies

$q\equiv 2 \pmod{3}.$

No prime \(p\equiv 1 \pmod{3}\) may appear inside \(u\), even squared, because that would increase some exponent \(e_p\) and change the product \(\prod (e_p+1)\).

Step 5: Fast Counting of Neutral Multipliers

For \(s\ge 1\), define

$G(s)=\left|\left\{u\le s:\forall q\mid u,\ q\equiv 2 \pmod{3}\right\}\right|.$

If \(x\) is the remaining budget after fixing the core \(c\), then the number of admissible multipliers \(t\le x\) is

$M(x)=\sum_{\substack{a\ge 0\\3^a\le x}} G\!\left(\left\lfloor\sqrt{\frac{x}{3^a}}\right\rfloor\right).$

This formula is exactly what the C++ method NeutralCounter::count_multipliers and the Python helper count_mult compute. The prefix table is built from a smallest-prime-factor sieve. The boolean recurrence used by the code is simple:

$\texttt{good}(1)=1,\qquad \texttt{good}(n)=1 \iff \left(\operatorname{spf}(n)\equiv 2 \pmod{3}\right)\land \texttt{good}\!\left(\frac{n}{\operatorname{spf}(n)}\right)=1.$

After prefix sums, each evaluation of \(G(s)\) is \(O(1)\), so each call to \(M(x)\) only iterates over powers of \(3\).

Step 6: Enumerating the Core Families

The smallest possible core is obtained from the smallest three primes congruent to \(1 \pmod{3}\):

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

This is why the solver returns \(0\) immediately when \(N<c_{\min}\). It also explains the constants

$7^4\cdot 13^4\cdot 19^2,\qquad 7^4\cdot 13^4$

stored in the C++ source as lower structural bounds.

The program sieves primes \(p\equiv 1 \pmod{3}\), precomputes \(p^2\) and \(p^4\) for all such primes, and stores the much shorter lists of primes with

$p^{14}\le N,\qquad p^{24}\le N,\qquad p^{74}\le N.$

Then it runs the ordered families

$[74],\ [24,2],\ [2,24],\ [14,4],\ [4,14],\ [4,4,2],\ [4,2,4],\ [2,4,4],$

and for each admissible core \(c\) it adds

$M\!\left(\left\lfloor\frac{N}{c}\right\rfloor\right).$

This is the whole algorithm: enumerate the rare core patterns, then attach every neutral multiplier that still fits below the limit.

Worked Checks from the Source Files

The C++ implementation validates the representation formula against brute-force lattice counting for every \(1\le n\le 120\). It also includes explicit problem checkpoints:

$r(1)=6 \quad\Rightarrow\quad B(\sqrt{3})=6,$

$r(7)=12 \quad\Rightarrow\quad B(\sqrt{21})=12.$

For the larger integer distance used in the validation block,

$L=111111111,\qquad n=\frac{L^2}{3}=3^3\cdot 37^2\cdot 333667^2,$

so the formula gives

$B(L)=r(n)=6(2+1)(2+1)=54,$

exactly matching the checkpoint in the repository.

How the Code Works

The C++ file is the full reference implementation. It contains the number-theoretic formula, brute-force validators, the neutral-multiplier sieve, precomputed power tables, and the final core enumeration. The expensive three-prime families are parallelized with parallel_sum_indices, but the mathematics is unchanged.

The Python file is a direct arithmetic transcription of the same method: the same exponent families, the same neutral counter, and the same final sum. The Java file does not reimplement the mathematics independently; instead, it compiles and runs the C++ solver and returns the parsed final output. So all three local solution files encode the same argument, just with different engineering choices.

Complexity Analysis

Let

$N=\left\lfloor\frac{L_{\max}^2}{3}\right\rfloor,\qquad P=\left\lfloor\sqrt{\frac{N}{7^4\cdot 13^4}}\right\rfloor,\qquad S=\left\lfloor\sqrt{\frac{N}{7^4\cdot 13^4\cdot 19^2}}\right\rfloor.$

Sieving primes \(p\equiv 1 \pmod{3}\) up to \(P\) costs near \(O(P\log\log P)\) time. Building the smallest-prime-factor table and the neutral prefix counts up to \(S\) also costs near \(O(S\log\log S)\) time. Both structures use linear memory in their respective limits.

After that preprocessing, each multiplier query \(M(x)\) uses only \(O(\log x)\) iterations over powers of \(3\), with \(O(1)\) prefix lookups inside. The family loops are heavily pruned because exponents \(74\), \(24\), \(14\), \(4\), and \(2\) make the core values grow very quickly. In practice the algorithm is dominated by the sieves and this pruned enumeration, not by any geometric search.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=354
2. Eisenstein integers and the norm \(a^2+ab+b^2\): https://en.wikipedia.org/wiki/Eisenstein_integer
3. Triangular lattice: https://en.wikipedia.org/wiki/Triangular_lattice
4. Hexagonal tiling / honeycomb geometry: https://en.wikipedia.org/wiki/Hexagonal_tiling
5. Binary quadratic forms: https://en.wikipedia.org/wiki/Binary_quadratic_form