Problem 780: Toriangulations
View on Project EulerProject Euler Problem 780 Solution
EulerSolve provides an optimized solution for Project Euler Problem 780, Toriangulations, with C++, Python, Java, and a step-by-step mathematical explanation.
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....
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
- Problem page: https://projecteuler.net/problem=780
- Möbius inversion formula: Wikipedia — Möbius inversion formula
- Divisor summatory function: Wikipedia — Divisor summatory function
- Beatty sequence: Wikipedia — Beatty sequence
- Eisenstein integer: Wikipedia — Eisenstein integer
- Hexagonal lattice: Wikipedia — Hexagonal lattice
Problem 780 source code
C++
#include <algorithm>
#include <atomic>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <thread>
#include <unordered_map>
#include <utility>
#include <vector>
using i64 = std::int64_t;
using u64 = std::uint64_t;
using i128 = __int128_t;
using u128 = __uint128_t;
constexpr i64 MOD = 1'000'000'007LL;
static i128 abs_i128(i128 x) {
return x < 0 ? -x : x;
}
static i128 gcd_i128(i128 a, i128 b) {
a = abs_i128(a);
b = abs_i128(b);
while (b != 0) {
const i128 t = a % b;
a = b;
b = t;
}
return a;
}
static i128 floor_div_i128(const i128 a, const i128 b) {
i128 q = a / b;
i128 r = a % b;
if (r < 0) {
q -= 1;
}
return q;
}
static i128 isqrt_i128(const i128 x) {
if (x <= 0) {
return 0;
}
long double approx = std::sqrt(static_cast<long double>(x));
i128 r = static_cast<i128>(approx);
while ((r + 1) * (r + 1) <= x) {
++r;
}
while (r * r > x) {
--r;
}
return r;
}
static i128 floor_sqrt3_mul(const i128 n) {
return isqrt_i128(3 * n * n);
}
static i128 floor_div_sqrt3(const i128 n, const i128 m) {
const i128 nn = n * n;
const i128 den = 3 * m * m;
i128 t = isqrt_i128(nn / den);
while (den * (t + 1) * (t + 1) <= nn) {
++t;
}
while (den * t * t > nn) {
--t;
}
return t;
}
static i128 divisor_summatory(const i64 n) {
i128 res = 0;
i64 i = 1;
while (i <= n) {
const i64 q = n / i;
const i64 j = n / q;
res += static_cast<i128>(q) * static_cast<i128>(j - i + 1);
i = j + 1;
}
return res;
}
struct Alpha {
i128 a;
i128 b;
i128 c;
};
static Alpha normalize_alpha(i128 a, i128 b, i128 c) {
if (c < 0) {
a = -a;
b = -b;
c = -c;
}
const i128 g = gcd_i128(gcd_i128(abs_i128(a), abs_i128(b)), c);
if (g > 1) {
a /= g;
b /= g;
c /= g;
}
return {a, b, c};
}
static i128 floor_qsqrt3(const i128 a, const i128 b, const i128 c) {
if (b == 0) {
return floor_div_i128(a, c);
}
const i128 fb = floor_sqrt3_mul(b);
i128 x = floor_div_i128(a + fb, c);
const i128 bb3 = 3 * b * b;
while (true) {
const i128 y = (x + 1) * c - a;
if (y <= 0 || y * y <= bb3) {
++x;
} else {
break;
}
}
while (true) {
const i128 y = x * c - a;
if (y <= 0 || y * y <= bb3) {
break;
}
--x;
}
return x;
}
static i128 alpha_floor(const Alpha& alpha) {
return floor_qsqrt3(alpha.a, alpha.b, alpha.c);
}
static i128 alpha_mul_floor(const Alpha& alpha, const i128 n) {
return floor_qsqrt3(alpha.a * n, alpha.b * n, alpha.c);
}
static Alpha alpha_sub_int(const Alpha& alpha, const i128 k) {
return normalize_alpha(alpha.a - k * alpha.c, alpha.b, alpha.c);
}
static Alpha alpha_div_alpha_minus_one(const Alpha& alpha) {
const i128 ac = alpha.a - alpha.c;
i128 A = alpha.a * ac - 3 * alpha.b * alpha.b;
i128 B = -alpha.b * alpha.c;
i128 D = ac * ac - 3 * alpha.b * alpha.b;
if (D < 0) {
A = -A;
B = -B;
D = -D;
}
if (B < 0) {
A = -A;
B = -B;
}
return normalize_alpha(A, B, D);
}
static i128 tri(const i128 n) {
return n * (n + 1) / 2;
}
static i128 beatty_sum(Alpha alpha, i128 n) {
i128 res = 0;
i64 sign = 1;
while (n > 0) {
const i128 f = alpha_floor(alpha);
if (f > 1) {
res += static_cast<i128>(sign) * (f - 1) * tri(n);
alpha = alpha_sub_int(alpha, f - 1);
}
const i128 m = alpha_mul_floor(alpha, n);
res += static_cast<i128>(sign) * tri(m);
n = m - n;
if (n <= 0) {
break;
}
alpha = alpha_div_alpha_minus_one(alpha);
sign = -sign;
}
return res;
}
static i128 beatty_sqrt3(const i128 c, const i128 n) {
if (n <= 0) {
return 0;
}
return beatty_sum({0, c, 1}, n);
}
static void linear_sieve(const int n,
std::vector<int>& spf,
std::vector<int>& mu) {
spf.assign(n + 1, 0);
mu.assign(n + 1, 0);
std::vector<int> primes;
mu[1] = 1;
spf[1] = 1;
for (int i = 2; i <= n; ++i) {
if (spf[i] == 0) {
spf[i] = i;
primes.push_back(i);
mu[i] = -1;
}
for (const int p : primes) {
const i64 v = static_cast<i64>(i) * static_cast<i64>(p);
if (v > n) {
break;
}
spf[static_cast<std::size_t>(v)] = p;
if (i % p == 0) {
mu[static_cast<std::size_t>(v)] = 0;
break;
}
mu[static_cast<std::size_t>(v)] = -mu[i];
}
}
}
static std::vector<int> distinct_prime_factors(int x,
const std::vector<int>& spf) {
std::vector<int> ps;
while (x > 1) {
const int p = spf[static_cast<std::size_t>(x)];
ps.push_back(p);
while (x % p == 0) {
x /= p;
}
}
return ps;
}
static std::vector<std::vector<std::pair<int, int>>> build_squarefree_divs(
const int max_n,
const std::vector<int>& spf
) {
std::vector<std::vector<std::pair<int, int>>> sf(max_n + 1);
sf[1] = {{1, 1}};
for (int n = 2; n <= max_n; ++n) {
const std::vector<int> ps = distinct_prime_factors(n, spf);
std::vector<std::pair<int, int>> divs;
divs.push_back({1, 1});
for (const int p : ps) {
const std::size_t cur = divs.size();
for (std::size_t i = 0; i < cur; ++i) {
divs.push_back({divs[i].first * p, -divs[i].second});
}
}
sf[n] = std::move(divs);
}
return sf;
}
static std::vector<std::vector<int>> build_all_divisors(
const int max_n,
const std::vector<int>& spf
) {
std::vector<std::vector<int>> divs(max_n + 1);
divs[1] = {1};
for (int n = 2; n <= max_n; ++n) {
int x = n;
const int p = spf[static_cast<std::size_t>(x)];
int e = 0;
while (x % p == 0) {
x /= p;
++e;
}
const std::vector<int>& base = divs[static_cast<std::size_t>(x)];
std::vector<int> out;
out.reserve(base.size() * static_cast<std::size_t>(e + 1));
int pe = 1;
for (int i = 0; i <= e; ++i) {
for (const int d : base) {
out.push_back(d * pe);
}
pe *= p;
}
divs[static_cast<std::size_t>(n)] = std::move(out);
}
return divs;
}
static unsigned choose_threads(const i64 work_items) {
if (work_items < 2'000) {
return 1;
}
unsigned threads = std::thread::hardware_concurrency();
if (threads == 0) {
threads = 1;
}
if (static_cast<i64>(threads) > work_items) {
threads = static_cast<unsigned>(work_items);
}
return std::max(1U, threads);
}
struct StripPartial {
i128 s1 = 0;
i128 s2 = 0;
i128 diag1 = 0;
i128 diag2 = 0;
};
static std::pair<i128, i128> strip_hyperbola_sum(
const i64 N,
const i64 V,
const i64 L,
const std::vector<std::vector<std::pair<int, int>>>& sf_divs,
const std::vector<std::vector<int>>& all_divs
) {
const unsigned threads = choose_threads(V);
std::vector<StripPartial> partial(threads);
std::vector<std::thread> workers;
workers.reserve(threads);
std::atomic<i64> next_v{1};
for (unsigned tid = 0; tid < threads; ++tid) {
workers.emplace_back([&, tid]() {
StripPartial local;
while (true) {
const i64 v = next_v.fetch_add(1, std::memory_order_relaxed);
if (v > V) {
break;
}
const i64 Umax = L / v;
const std::vector<int>& dv = all_divs[static_cast<std::size_t>(v)];
for (const i64 d : dv) {
const i64 m = v / d;
const i64 hi = Umax / d;
if (hi < m) {
continue;
}
const i64 w = N / (2 * d);
const i64 lo_minus = m - 1;
i128 cnt = 0;
i128 sm = 0;
const auto& sfm = sf_divs[static_cast<std::size_t>(m)];
for (const auto& [q, muq] : sfm) {
cnt += static_cast<i128>(muq) *
(static_cast<i128>(hi / q) - static_cast<i128>(lo_minus / q));
const i64 hiq = hi / q;
const i64 loq = lo_minus / q;
if (hiq > 0) {
const i128 c = static_cast<i128>(v) * static_cast<i128>(q);
sm += static_cast<i128>(muq) *
(beatty_sqrt3(c, hiq) - beatty_sqrt3(c, loq));
}
}
local.s1 += static_cast<i128>(w) * cnt;
local.s2 += sm;
}
local.diag1 += static_cast<i128>(N / (2 * v));
local.diag2 += floor_sqrt3_mul(v);
}
partial[tid] = local;
});
}
for (std::thread& worker : workers) {
worker.join();
}
i128 s1 = 0;
i128 s2 = 0;
i128 diag1 = 0;
i128 diag2 = 0;
for (const StripPartial& p : partial) {
s1 += p.s1;
s2 += p.s2;
diag1 += p.diag1;
diag2 += p.diag2;
}
const i128 S1 = 2 * s1 - diag1;
const i128 S2 = 2 * s2 - diag2;
return {S1, S2};
}
static i64 count_mod3_res(const i64 lo, const i64 hi, const i64 r) {
if (hi < lo) {
return 0;
}
const i64 rem = ((lo % 3) + 3) % 3;
const i64 delta = (r - rem + 3) % 3;
const i64 first = lo + delta;
if (first > hi) {
return 0;
}
return (hi - first) / 3 + 1;
}
static i128 hex_Hsum(
const i64 X,
const std::vector<std::vector<std::pair<int, int>>>& sf_divs
) {
if (X <= 0) {
return 0;
}
std::unordered_map<i64, i128> d_cache;
d_cache.reserve(1 << 15);
auto D = [&](const i64 n) -> i128 {
if (n <= 0) {
return 0;
}
const auto it = d_cache.find(n);
if (it != d_cache.end()) {
return it->second;
}
const i128 val = divisor_summatory(n);
d_cache.emplace(n, val);
return val;
};
const i64 V = static_cast<i64>(isqrt_i128(X));
i128 extra = 0;
for (i64 v = 1; v <= V; ++v) {
const i128 disc0 = 4 * static_cast<i128>(X) - 3 * static_cast<i128>(v) * static_cast<i128>(v);
if (disc0 <= 0) {
break;
}
const i64 Umax = static_cast<i64>((-static_cast<i128>(v) + isqrt_i128(disc0)) / 2);
i64 u = v + 1;
const auto& sfv = sf_divs[static_cast<std::size_t>(v)];
const i64 vmod = v % 3;
while (u <= Umax) {
const i128 t = static_cast<i128>(u) * u + static_cast<i128>(u) * v + static_cast<i128>(v) * v;
const i64 q = static_cast<i64>(X / t);
if (q == 0) {
break;
}
const i64 T = X / q;
const i128 disc = 4 * static_cast<i128>(T) - 3 * static_cast<i128>(v) * static_cast<i128>(v);
i64 uhi = static_cast<i64>((-static_cast<i128>(v) + isqrt_i128(disc)) / 2);
if (uhi > Umax) {
uhi = Umax;
}
const i64 lo1 = u - 1;
i64 total = 0;
for (const auto& [d, mud] : sfv) {
total += mud * (uhi / d - lo1 / d);
}
if (vmod != 0) {
i64 bad = 0;
for (const auto& [d, mud] : sfv) {
const i64 dm3 = d % 3;
const i64 inv = (dm3 == 1 ? 1 : 2);
const i64 tlo = (u + d - 1) / d;
const i64 thi = uhi / d;
const i64 rr = (vmod * inv) % 3;
bad += mud * count_mod3_res(tlo, thi, rr);
}
total -= bad;
}
if (total != 0) {
extra += static_cast<i128>(total) * D(q);
}
u = uhi + 1;
}
}
return D(X) + 2 * extra;
}
static i64 solve_fast(const i64 N) {
if (N <= 0) {
return 0;
}
const i64 M = N / 2;
const i64 L = static_cast<i64>(floor_div_sqrt3(M, 1));
const i64 V_strip = static_cast<i64>(isqrt_i128(L));
const i64 X = N / 4;
const i64 V_hex = (X > 0 ? static_cast<i64>(isqrt_i128(X)) : 0);
const int max_pre = static_cast<int>(std::max<i64>({1, V_strip, V_hex}));
std::vector<int> spf;
std::vector<int> mu;
linear_sieve(max_pre, spf, mu);
const auto sf_divs = build_squarefree_divs(max_pre, spf);
const auto all_divs = build_all_divisors(max_pre, spf);
const i128 base = 2 * divisor_summatory(M);
const auto [S1, S2] = strip_hyperbola_sum(N, V_strip, L, sf_divs, all_divs);
const i128 strip_part = base + 4 * (S1 - S2);
const i128 H = hex_Hsum(X, sf_divs);
const i128 total = strip_part - 4 * H;
i128 modded = total % MOD;
if (modded < 0) {
modded += MOD;
}
return static_cast<i64>(modded);
}
int main() {
const i64 g5 = solve_fast(5);
const i64 g6 = solve_fast(6);
assert(((g6 - g5) % MOD + MOD) % MOD == 8);
assert(g6 == 14);
assert(solve_fast(100) == 8'090);
assert(solve_fast(100'000) == 645'124'048);
std::cout << solve_fast(1'000'000'000) << '\n';
return 0;
}
Python
import math
def solve():
MOD = 10**9+7; N = 10**9
def isqrt(x):
r=int(math.isqrt(x))
while (r+1)*(r+1)<=x: r+=1
while r*r>x: r-=1
return r
def isqrt_big(x):
if x<=0: return 0
r=int(math.isqrt(x))
while (r+1)*(r+1)<=x: r+=1
while r*r>x: r-=1
return r
def floor_sqrt3_mul(n):
return isqrt_big(3*n*n)
def floor_div_sqrt3(n,m):
nn=n*n; den=3*m*m
t=isqrt_big(nn//den)
while den*(t+1)*(t+1)<=nn: t+=1
while den*t*t>nn: t-=1
return t
def div_sum(n):
res=0; i=1
while i<=n:
q=n//i; j=n//q
res+=q*(j-i+1); i=j+1
return res
def tri(n): return n*(n+1)//2
def floor_qsqrt3(a,b,c):
if b==0: return a//c if (a>=0 or a%c==0) else -((-a+c-1)//c)
fb=isqrt_big(3*b*b)
x=(a+fb)//c
bb3=3*b*b
while True:
y=(x+1)*c-a
if y<=0 or y*y<=bb3: x+=1
else: break
while True:
y=x*c-a
if y<=0 or y*y<=bb3: break
x-=1
return x
def beatty_sum(a,b,c,n):
res=0; sign=1
while n>0:
f=floor_qsqrt3(a,b,c)
if f>1:
res+=sign*(f-1)*tri(n)
a-=(f-1)*c
m=floor_qsqrt3(a*n,b*n,c)
res+=sign*tri(m)
n=m-n
if n<=0: break
# alpha/(alpha-1): alpha = (a+b*sqrt3)/c
ac=a-c
A=a*ac-3*b*b; B=-b*c; D=ac*ac-3*b*b
if D<0: A,B,D=-A,-B,-D
if B<0: A,B=-A,-B
from math import gcd
g=gcd(gcd(abs(A),abs(B)),abs(D))
if g>1: A//=g; B//=g; D//=g
a,b,c=A,B,D; sign=-sign
return res
def beatty_sqrt3(cv,n):
if n<=0: return 0
return beatty_sum(0,cv,1,n)
def linear_sieve(n):
spf=[0]*(n+1); mu=[0]*(n+1); primes=[]; mu[1]=1; spf[1]=1
for i in range(2,n+1):
if spf[i]==0:
spf[i]=i; primes.append(i); mu[i]=-1
for p in primes:
if i*p>n: break
spf[i*p]=p
if i%p==0: mu[i*p]=0; break
mu[i*p]=-mu[i]
return spf,mu
def sf_divs(n,spf):
divs=[[] for _ in range(n+1)]
divs[1]=[(1,1)]
for x in range(2,n+1):
v=x; ps=[]
while v>1:
p=spf[v]; ps.append(p)
while v%p==0: v//=p
k=len(ps); dd=[]
for mask in range(1<<k):
d=1; bc=0
for i in range(k):
if mask&(1<<i): d*=ps[i]; bc+=1
dd.append((d,-1 if bc&1 else 1))
divs[x]=dd
return divs
def all_divs(n,spf):
divs=[[] for _ in range(n+1)]
divs[1]=[1]
for x in range(2,n+1):
v=x; p=spf[v]; e=0
while v%p==0: v//=p; e+=1
base=divs[v]; out=[]
pe=1
for i in range(e+1):
for d in base: out.append(d*pe)
pe*=p
divs[x]=out
return divs
M=N//2; L=floor_div_sqrt3(M,1); V_strip=isqrt_big(L)
X=N//4; V_hex=isqrt_big(X) if X>0 else 0
mx=max(1,V_strip,V_hex)
spf,mu=linear_sieve(mx)
sfd=sf_divs(mx,spf); ald=all_divs(mx,spf)
base=2*div_sum(M)
s1=0; s2=0; diag1=0; diag2=0
for v in range(1,V_strip+1):
Umax=L//v
for d in ald[v]:
m=v//d; hi=Umax//d
if hi<m: continue
w=N//(2*d); lo=m-1
cnt=0; sm=0
for q,muq in sfd[m]:
cnt+=muq*(hi//q-lo//q)
hiq=hi//q; loq=lo//q
if hiq>0:
c=v*q
sm+=muq*(beatty_sqrt3(c,hiq)-beatty_sqrt3(c,loq))
s1+=w*cnt; s2+=sm
diag1+=N//(2*v); diag2+=floor_sqrt3_mul(v)
S1=2*s1-diag1; S2=2*s2-diag2
strip=base+4*(S1-S2)
# Hex sum
d_cache={}
def D(n2):
if n2<=0: return 0
if n2 in d_cache: return d_cache[n2]
v=div_sum(n2); d_cache[n2]=v; return v
def cnt_mod3(lo,hi,r):
if hi<lo: return 0
rem=lo%3
if rem<0: rem+=3
delta=(r-rem)%3
first=lo+delta
if first>hi: return 0
return (hi-first)//3+1
extra=0
for v in range(1,V_hex+1):
disc0=4*X-3*v*v
if disc0<=0: break
Umax=(-v+isqrt_big(disc0))//2; u=v+1
sfv=sfd[v]; vm=v%3
while u<=Umax:
t=u*u+u*v+v*v; q=X//t
if q==0: break
T=X//q
disc=4*T-3*v*v
uhi=min((-v+isqrt_big(disc))//2,Umax)
lo1=u-1; total=0
for d,mud in sfv:
total+=mud*(uhi//d-lo1//d)
if vm!=0:
bad=0
for d,mud in sfv:
dm3=d%3; inv=1 if dm3==1 else 2
tlo=(u+d-1)//d; thi=uhi//d
rr=(vm*inv)%3
bad+=mud*cnt_mod3(tlo,thi,rr)
total-=bad
if total!=0: extra+=total*D(q)
u=uhi+1
H=D(X)+2*extra
total=strip-4*H
return str(total%MOD)
if __name__=='__main__':
print(solve())
Java
import java.util.ArrayList;
import java.util.HashMap;
public class Euler780 {
static final long MOD = 1000000007L;
static long isqrt(long x) {
if (x <= 0)
return 0;
long r = (long) Math.sqrt((double) x);
while ((r + 1) * (r + 1) <= x)
r++;
while (r * r > x && r > 0)
r--;
return r;
}
static long floorSqrt3Mul(long n) {
return isqrt(3L * n * n);
}
static long floorDivSqrt3(long n, long m) {
long nn = n * n;
long den = 3L * m * m;
long t = isqrt(nn / den);
while (den * (t + 1) * (t + 1) <= nn)
t++;
while (den * t * t > nn && t > 0)
t--;
return t;
}
static long divisorSummatory(long n) {
long res = 0;
long i = 1;
while (i <= n) {
long q = n / i;
long j = n / q;
res += q * (j - i + 1);
i = j + 1;
}
return res;
}
static long gcd(long a, long b) {
a = Math.abs(a);
b = Math.abs(b);
while (b != 0) {
long t = a % b;
a = b;
b = t;
}
return a;
}
static class Alpha {
long a, b, c;
Alpha(long a, long b, long c) {
this.a = a;
this.b = b;
this.c = c;
}
}
static Alpha normalizeAlpha(long a, long b, long c) {
if (c < 0) {
a = -a;
b = -b;
c = -c;
}
long g = gcd(gcd(a, b), c);
if (g > 1) {
a /= g;
b /= g;
c /= g;
}
return new Alpha(a, b, c);
}
static long floorDiv(long a, long b) {
long q = a / b;
long r = a % b;
if ((r < 0 && b > 0) || (r > 0 && b < 0))
q--;
return q;
}
static long floorQsqrt3(long a, long b, long c) {
if (b == 0)
return floorDiv(a, c);
long fb = floorSqrt3Mul(b);
long x = floorDiv(a + fb, c);
long bb3 = 3L * b * b;
while (true) {
long y = (x + 1) * c - a;
if (y <= 0 || y * y <= bb3)
x++;
else
break;
}
while (true) {
long y = x * c - a;
if (y <= 0 || y * y <= bb3)
break;
x--;
}
return x;
}
static long alphaFloor(Alpha alpha) {
return floorQsqrt3(alpha.a, alpha.b, alpha.c);
}
static long alphaMulFloor(Alpha alpha, long n) {
return floorQsqrt3(alpha.a * n, alpha.b * n, alpha.c);
}
static Alpha alphaSubInt(Alpha alpha, long k) {
return normalizeAlpha(alpha.a - k * alpha.c, alpha.b, alpha.c);
}
static Alpha alphaDivAlphaMinusOne(Alpha alpha) {
long ac = alpha.a - alpha.c;
long A = alpha.a * ac - 3L * alpha.b * alpha.b;
long B = -alpha.b * alpha.c;
long D = ac * ac - 3L * alpha.b * alpha.b;
if (D < 0) {
A = -A;
B = -B;
D = -D;
}
if (B < 0) {
A = -A;
B = -B;
}
return normalizeAlpha(A, B, D);
}
static long tri(long n) {
return n * (n + 1) / 2L;
}
static long beattySum(Alpha alpha, long n) {
long res = 0;
long sign = 1;
while (n > 0) {
long f = alphaFloor(alpha);
if (f > 1) {
res += sign * (f - 1) * tri(n);
alpha = alphaSubInt(alpha, f - 1);
}
long m = alphaMulFloor(alpha, n);
res += sign * tri(m);
n = m - n;
if (n <= 0)
break;
alpha = alphaDivAlphaMinusOne(alpha);
sign = -sign;
}
return res;
}
static long beattySqrt3(long c, long n) {
if (n <= 0)
return 0;
return beattySum(new Alpha(0, c, 1), n);
}
static void linearSieve(int n, int[] spf, int[] mu) {
ArrayList<Integer> primes = new ArrayList<>();
mu[1] = 1;
spf[1] = 1;
for (int i = 2; i <= n; i++) {
if (spf[i] == 0) {
spf[i] = i;
primes.add(i);
mu[i] = -1;
}
for (int p : primes) {
long v = (long) i * p;
if (v > n)
break;
spf[(int) v] = p;
if (i % p == 0) {
mu[(int) v] = 0;
break;
}
mu[(int) v] = -mu[i];
}
}
}
static class Pair {
int first;
int second;
Pair(int first, int second) {
this.first = first;
this.second = second;
}
}
static ArrayList<ArrayList<Pair>> buildSquarefreeDivs(int maxN, int[] spf) {
ArrayList<ArrayList<Pair>> sf = new ArrayList<>(maxN + 1);
for (int i = 0; i <= maxN; i++)
sf.add(null);
ArrayList<Pair> initial = new ArrayList<>();
initial.add(new Pair(1, 1));
sf.set(1, initial);
for (int n = 2; n <= maxN; n++) {
ArrayList<Integer> ps = new ArrayList<>();
int x = n;
while (x > 1) {
int p = spf[x];
ps.add(p);
while (x % p == 0)
x /= p;
}
ArrayList<Pair> divs = new ArrayList<>();
divs.add(new Pair(1, 1));
for (int p : ps) {
int cur = divs.size();
for (int i = 0; i < cur; i++) {
divs.add(new Pair(divs.get(i).first * p, -divs.get(i).second));
}
}
sf.set(n, divs);
}
return sf;
}
static ArrayList<ArrayList<Integer>> buildAllDivisors(int maxN, int[] spf) {
ArrayList<ArrayList<Integer>> divs = new ArrayList<>(maxN + 1);
for (int i = 0; i <= maxN; i++)
divs.add(null);
ArrayList<Integer> initial = new ArrayList<>();
initial.add(1);
divs.set(1, initial);
for (int n = 2; n <= maxN; n++) {
int x = n;
int p = spf[x];
int e = 0;
while (x % p == 0) {
x /= p;
e++;
}
ArrayList<Integer> base = divs.get(x);
ArrayList<Integer> out = new ArrayList<>();
int pe = 1;
for (int i = 0; i <= e; i++) {
for (int d : base) {
out.add(d * pe);
}
pe *= p;
}
divs.set(n, out);
}
return divs;
}
static long[] stripHyperbolaSum(long N, long V, long L, ArrayList<ArrayList<Pair>> sfDivs,
ArrayList<ArrayList<Integer>> allDivs) {
long s1 = 0, s2 = 0, diag1 = 0, diag2 = 0;
for (int v = 1; v <= V; v++) {
long Umax = L / v;
ArrayList<Integer> dv = allDivs.get(v);
for (long d : dv) {
long m = v / d;
long hi = Umax / d;
if (hi < m)
continue;
long w = N / (2L * d);
long loMinus = m - 1;
long cnt = 0;
long sm = 0;
ArrayList<Pair> sfm = sfDivs.get((int) m);
for (Pair p : sfm) {
long q = p.first;
long muq = p.second;
cnt += muq * ((hi / q) - (loMinus / q));
long hiq = hi / q;
long loq = loMinus / q;
if (hiq > 0) {
long c = (long) v * q;
sm += muq * (beattySqrt3(c, hiq) - beattySqrt3(c, loq));
}
}
s1 += w * cnt;
s2 += sm;
}
diag1 += N / (2L * v);
diag2 += floorSqrt3Mul(v);
}
long S1 = 2L * s1 - diag1;
long S2 = 2L * s2 - diag2;
return new long[] { S1, S2 };
}
static long countMod3Res(long lo, long hi, long r) {
if (hi < lo)
return 0;
long rem = ((lo % 3L) + 3L) % 3L;
long delta = (r - rem + 3L) % 3L;
long first = lo + delta;
if (first > hi)
return 0;
return (hi - first) / 3L + 1L;
}
static long hexHsum(long X, ArrayList<ArrayList<Pair>> sfDivs) {
if (X <= 0)
return 0;
HashMap<Long, Long> dCache = new HashMap<>();
long V = isqrt(X);
long extra = 0;
for (long v = 1; v <= V; v++) {
long disc0 = 4L * X - 3L * v * v;
if (disc0 <= 0)
break;
long Umax = (-v + isqrt(disc0)) / 2L;
long u = v + 1;
ArrayList<Pair> sfv = sfDivs.get((int) v);
long vmod = v % 3L;
while (u <= Umax) {
long t = u * u + u * v + v * v;
long q = X / t;
if (q == 0)
break;
long T = X / q;
long disc = 4L * T - 3L * v * v;
long uhi = (-v + isqrt(disc)) / 2L;
if (uhi > Umax)
uhi = Umax;
long lo1 = u - 1;
long total = 0;
for (Pair p : sfv) {
long d = p.first;
long mud = p.second;
total += mud * (uhi / d - lo1 / d);
}
if (vmod != 0) {
long bad = 0;
for (Pair p : sfv) {
long d = p.first;
long mud = p.second;
long dm3 = d % 3L;
long inv = (dm3 == 1 ? 1L : 2L);
long tlo = (u + d - 1L) / d;
long thi = uhi / d;
long rr = (vmod * inv) % 3L;
bad += mud * countMod3Res(tlo, thi, rr);
}
total -= bad;
}
if (total != 0) {
Long cacheVal = dCache.get(q);
long dq;
if (cacheVal != null) {
dq = cacheVal;
} else {
dq = divisorSummatory(q);
dCache.put(q, dq);
}
extra += total * dq;
}
u = uhi + 1;
}
}
Long cacheValX = dCache.get(X);
long dx = (cacheValX != null) ? cacheValX : divisorSummatory(X);
return dx + 2L * extra;
}
static long solveFast(long N) {
if (N <= 0)
return 0;
long M = N / 2L;
long L = floorDivSqrt3(M, 1L);
long VStrip = isqrt(L);
long X = N / 4L;
long VHex = (X > 0 ? isqrt(X) : 0);
int maxPre = (int) Math.max(1, Math.max(VStrip, VHex));
int[] spf = new int[maxPre + 1];
int[] mu = new int[maxPre + 1];
linearSieve(maxPre, spf, mu);
ArrayList<ArrayList<Pair>> sfDivs = buildSquarefreeDivs(maxPre, spf);
ArrayList<ArrayList<Integer>> allDivs = buildAllDivisors(maxPre, spf);
long base = 2L * divisorSummatory(M);
long[] S = stripHyperbolaSum(N, VStrip, L, sfDivs, allDivs);
long S1 = S[0], S2 = S[1];
long stripPart = base + 4L * (S1 - S2);
long H = hexHsum(X, sfDivs);
long total = stripPart - 4L * H;
long modded = total % MOD;
if (modded < 0)
modded += MOD;
return modded;
}
public static String solve() {
return Long.toString(solveFast(1000000000L));
}
public static void main(String[] args) {
System.out.println(solve());
}
}