Detailed mathematical approach

Problem Summary

Problem 780 asks for the arithmetic counting function \(G(N)\) associated with toriangulations, evaluated at \(N=10^9\) modulo \(10^9+7\). The C++, Python, and Java implementations show that the geometric statement can be rewritten as an exact count of primitive directions in the triangular lattice, so the task becomes a problem about gcd filters, divisor sums, and the quadratic form \(u^2+uv+v^2\).

Instead of enumerating candidate configurations directly, the solution splits the total into an axial contribution, a strip-hyperbola contribution, and a hexagonal correction. Each of those pieces can be evaluated with arithmetic tools that scale roughly like \(O(\sqrt N\,\mathrm{polylog}\,N)\) rather than linearly in \(N\).

Mathematical Approach

The implementations reduce the final answer to an exact arithmetic decomposition. Define

$m_0=\left\lfloor\frac N2\right\rfloor,\qquad \lambda=\left\lfloor\frac{m_0}{\sqrt3}\right\rfloor,\qquad x_0=\left\lfloor\frac N4\right\rfloor,$

and let

$D(t)=\sum_{k=1}^{t}\left\lfloor\frac{t}{k}\right\rfloor.$

Then the count is organized as

$G(N)=2D(m_0)+4\bigl(A(\lambda)-B(\lambda)\bigr)-4C(x_0).$

Step 1: Separate the Three Arithmetic Pieces

The term \(2D(m_0)\) accounts for the two axial families that are already visible in one-dimensional divisor counting. The remaining work comes from genuinely two-dimensional directions in the triangular lattice.

The strip part is captured by two sums over positive integer pairs with \(uv\le \lambda\):

$A(\lambda)=\sum_{uv\le \lambda}\left\lfloor\frac{N}{2\gcd(u,v)}\right\rfloor,$

$B(\lambda)=\sum_{uv\le \lambda}\left\lfloor\frac{\sqrt3\,uv}{\gcd(u,v)}\right\rfloor.$

The first counts admissible widths after separating out the common divisor of \((u,v)\), while the second subtracts the exact \(\sqrt3\)-boundary coming from the same primitive data.

Step 2: Rewrite Coprimality with Möbius Inversion

To evaluate the strip sums efficiently, write

$u=d\,a,\qquad v=d\,b,\qquad \gcd(a,b)=1.$

For fixed \(v\) and fixed divisor \(d\mid v\), the reduced parameter \(b=v/d\) is fixed, and the admissible values of \(a\) lie in an interval determined by \(uv\le \lambda\). The coprimality condition is converted into an inclusion-exclusion sum

$\mathbf{1}_{\gcd(a,b)=1}=\sum_{q\mid b}\mu(q)\,\mathbf{1}_{q\mid a},$

where \(\mu\) is the Möbius function. This turns a scan over individual \(a\) values into a sum over the squarefree divisors of \(b\), which is exactly what the implementations precompute.

Step 3: Evaluate the \(\sqrt3\)-Weighted Boundary Exactly

After the Möbius step, the weighted part of the strip term becomes sums of the form

$\sum_{k=1}^{n}\left\lfloor c\,k\sqrt3\right\rfloor,$

for integer \(c\). These are Beatty-type floor sums. The implementations do not approximate them numerically; instead they apply an exact recursive floor-sum transformation for quadratic irrationals.

That matters because the final answer is sensitive to unit-size errors. The recursive reduction keeps everything in integer arithmetic, using only corrected integer square roots and algebraic transformations, so the strip contribution remains exact.

Step 4: Add the Hexagonal Norm Correction

The second geometric regime is controlled by the Eisenstein norm

$Q(u,v)=u^2+uv+v^2.$

The correction term can be written as

$C(x_0)=D(x_0)+2\sum_{\substack{u>v\ge 1\\Q(u,v)\le x_0\\\gcd(u,v)=1\\u\not\equiv v\pmod 3}} D\!\left(\left\lfloor\frac{x_0}{Q(u,v)}\right\rfloor\right).$

The condition \(u\not\equiv v\pmod 3\) removes the extra common factor coming from the Eisenstein prime above \(3\). In the implementations, that residue-class exclusion is handled arithmetically: first count all coprime values, then subtract the bad congruence class modulo \(3\).

Step 5: Compress the Hexagonal Sum by Quotient Plateaus

For fixed \(v\), the quotient

$\left\lfloor\frac{x_0}{Q(u,v)}\right\rfloor$

stays constant on intervals of \(u\). The implementations therefore move through the admissible \(u\)-range in blocks rather than one value at a time. Each block contributes a multiplicity times one cached divisor-summatory value.

This is the same philosophy as the Dirichlet hyperbola method: group equal floor quotients, evaluate each group once, and reuse the result.

Worked Example: \(N=10\)

For a small sanity check, take \(N=10\). Then

$m_0=\left\lfloor\frac{10}{2}\right\rfloor=5,\qquad \lambda=\left\lfloor\frac{5}{\sqrt3}\right\rfloor=2,\qquad x_0=\left\lfloor\frac{10}{4}\right\rfloor=2.$

Also

$D(5)=\left\lfloor\frac51\right\rfloor+\left\lfloor\frac52\right\rfloor+\left\lfloor\frac53\right\rfloor+\left\lfloor\frac54\right\rfloor+\left\lfloor\frac55\right\rfloor=5+2+1+1+1=10.$

The ordered pairs with \(uv\le 2\) are \((1,1)\), \((1,2)\), and \((2,1)\). Since all three pairs have gcd \(1\),

$A(2)=3\left\lfloor\frac{10}{2}\right\rfloor=15.$

For the weighted term,

$B(2)=\lfloor\sqrt3\rfloor+\lfloor2\sqrt3\rfloor+\lfloor2\sqrt3\rfloor=1+3+3=7.$

Finally, \(Q(2,1)=7>2\), so the inner hexagonal sum is empty and

$C(2)=D(2)=\left\lfloor\frac21\right\rfloor+\left\lfloor\frac22\right\rfloor=3.$

Therefore

$G(10)=2\cdot 10+4(15-7)-4\cdot 3=20+32-12=40.$

This example mirrors the exact decomposition used at the full input size.

How the Code Works

The C++, Python, and Java implementations first compute the three cutoffs \(m_0\), \(\lambda\), and \(x_0\). They then precompute smallest prime factors up to \(\max(\lfloor\sqrt{\lambda}\rfloor,\lfloor\sqrt{x_0}\rfloor)\), together with two divisor tables: all divisors of each integer, and squarefree divisors carrying Möbius signs.

For the strip contribution, the implementation loops over the second coordinate \(v\), splits by each divisor \(d\mid v\), and uses Möbius inversion to count reduced partners \(a\) with \(\gcd(a,v/d)=1\). The corresponding \(\sqrt3\)-weighted floor sums are evaluated by the exact Beatty recursion instead of direct iteration.

For the hexagonal correction, the implementation scans pairs ordered by \(v\) and advances \(u\) in maximal ranges where \(\left\lfloor x_0/Q(u,v)\right\rfloor\) is constant. Each distinct divisor-summatory query is cached, so repeated values are computed once and reused. The final arithmetic combination is then reduced modulo \(10^9+7\).

The three language versions follow the same mathematics. The C++ implementation additionally spreads the strip accumulation across multiple worker threads, while the Python and Java versions keep the same logic in serial form.

Complexity Analysis

Let \(P=\max(\lfloor\sqrt{\lambda}\rfloor,\lfloor\sqrt{x_0}\rfloor)\). Building the factor tables costs \(O(P\log\log P)\) time and \(O(P)\) memory. The strip and hexagonal phases both exploit divisor compression, quotient plateaus, and cached divisor-summatory values, so they avoid scanning all candidate pairs explicitly.

In practice the total running time behaves like \(O(\sqrt N\,\mathrm{polylog}\,N)\), with modest polylogarithmic factors coming from divisor enumeration and recursive floor-sum evaluation. Memory usage remains \(O(\sqrt N)\).

Footnotes and References

1. Problem page: https://projecteuler.net/problem=780
2. Möbius inversion formula: Wikipedia — Möbius inversion formula
3. Divisor summatory function: Wikipedia — Divisor summatory function
4. Beatty sequence: Wikipedia — Beatty sequence
5. Eisenstein integer: Wikipedia — Eisenstein integer
6. Hexagonal lattice: Wikipedia — Hexagonal lattice