Problem 880: Nested Radicals
View on Project EulerProject Euler Problem 880 Solution
EulerSolve provides an optimized solution for Project Euler Problem 880, Nested Radicals, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary The implementations compute a modular quantity \(H(N)\) for \(N=10^{15}\), with modulus $$M=1031^3+2.$$ Instead of scanning all integer pairs \((x,y)\) with \(\max(|x|,|y|)\le N\), the solution starts from the algebraic classification used for the nested-radical identities behind the problem. Once that classification is available, the task becomes: $$H(N)=\sum_{\substack{(x,y)\ \text{admissible}\\ \max(|x|,|y|)\le N}} (|x|+|y|)\pmod{M}.$$ The key point is that every admissible pair comes from a primitive parameter pair together with a square scaling factor, so the search can be reorganized arithmetically rather than geometrically. Mathematical Approach The implementations use two primitive families, depending on the parity of one parameter. To avoid exposing code symbols, let the coprime parameters be \(a,b\in \mathbb{Z}_{>0}\) with \(a\ne 2b\), and let the common scale be \(t\ge 1\). Step 1: Primitive Parameterization When \(b\) is odd, the primitive pair is $$X_{\mathrm{odd}}=b(b+4a)^3,\qquad Y_{\mathrm{odd}}=4a(a-2b)^3.$$ When \(b\) is even, the primitive pair is $$X_{\mathrm{even}}=2b\left(\frac{b}{2}+2a\right)^3,\qquad Y_{\mathrm{even}}=a(a-2b)^3.$$ These two branches are exactly the two branches used by the C++, Python, and Java implementations....
Detailed mathematical approach
Problem Summary
The implementations compute a modular quantity \(H(N)\) for \(N=10^{15}\), with modulus
$$M=1031^3+2.$$
Instead of scanning all integer pairs \((x,y)\) with \(\max(|x|,|y|)\le N\), the solution starts from the algebraic classification used for the nested-radical identities behind the problem. Once that classification is available, the task becomes:
$$H(N)=\sum_{\substack{(x,y)\ \text{admissible}\\ \max(|x|,|y|)\le N}} (|x|+|y|)\pmod{M}.$$
The key point is that every admissible pair comes from a primitive parameter pair together with a square scaling factor, so the search can be reorganized arithmetically rather than geometrically.
Mathematical Approach
The implementations use two primitive families, depending on the parity of one parameter. To avoid exposing code symbols, let the coprime parameters be \(a,b\in \mathbb{Z}_{>0}\) with \(a\ne 2b\), and let the common scale be \(t\ge 1\).
Step 1: Primitive Parameterization
When \(b\) is odd, the primitive pair is
$$X_{\mathrm{odd}}=b(b+4a)^3,\qquad Y_{\mathrm{odd}}=4a(a-2b)^3.$$
When \(b\) is even, the primitive pair is
$$X_{\mathrm{even}}=2b\left(\frac{b}{2}+2a\right)^3,\qquad Y_{\mathrm{even}}=a(a-2b)^3.$$
These two branches are exactly the two branches used by the C++, Python, and Java implementations. The coprimality condition
$$\gcd(a,b)=1$$
eliminates redundant representations, and the excluded line
$$a=2b$$
is the degenerate case where one cubic factor vanishes.
Step 2: Square Scaling Generates the Full Family
Each primitive pair produces an infinite family by multiplying both coordinates by a common square:
$$ (x,y)=\left(t^2X_\ast,\ t^2Y_\ast\right), $$
where \((X_\ast,Y_\ast)\) means the appropriate primitive pair from the odd or even branch.
Therefore, for a fixed primitive seed, the admissible values under the box constraint are determined by
$$t^2\max(|X_\ast|,|Y_\ast|)\le N,$$
so the largest allowed scale is
$$t_{\max}=\left\lfloor \sqrt{\frac{N}{\max(|X_\ast|,|Y_\ast|)}} \right\rfloor.$$
This is why the outer arithmetic work is done only on primitive seeds; all non-primitive solutions are absorbed into a single closed-form sum over \(t\).
Step 3: Cube-Free Kernels Remove the Collapsing Cases
The implementations precompute the cube-free kernel
$$\operatorname{cf}(n)=\prod_{p} p^{v_p(n)\bmod 3},$$
meaning that every full cube \(p^3\) is stripped away, while exponents \(1\) and \(2\) remain. The primitive pair is kept only if the two relevant kernels differ.
For odd \(b\), the required inequality is
$$\operatorname{cf}(b)\ne \operatorname{cf}(4a).$$
For even \(b\), the required inequality is
$$\operatorname{cf}(2b)\ne \operatorname{cf}(a).$$
This filter removes the cases where the two cube-root components collapse to the same cube-free part and the intended nested-radical structure is lost.
Step 4: Sum All Square Multiples at Once
For one primitive seed, every scaled pair contributes
$$|t^2X_\ast|+|t^2Y_\ast|=t^2\bigl(|X_\ast|+|Y_\ast|\bigr).$$
Hence the total contribution of that seed is
$$\bigl(|X_\ast|+|Y_\ast|\bigr)\sum_{t=1}^{t_{\max}} t^2.$$
The implementations use the standard closed form
$$\sum_{t=1}^{u} t^2=\frac{u(u+1)(2u+1)}{6},$$
so the whole scaled family is compressed into one arithmetic term rather than one loop over each individual pair.
Step 5: Root Bounds Make the Search Finite
The fourth-root bound for the outer parameter comes from the fact that the positive coordinate in each branch already grows like \(b^4\). In particular, the formulas imply
$$b\le \left\lfloor (4N)^{1/4}\right\rfloor.$$
Once \(b\) is fixed, the positive coordinate alone gives an upper bound for \(a\). For odd \(b\),
$$b(b+4a)^3\le N \quad\Longrightarrow\quad a\le \frac{\left\lfloor (N/b)^{1/3}\right\rfloor-b}{4}.$$
For even \(b\),
$$2b\left(\frac{b}{2}+2a\right)^3\le N \quad\Longrightarrow\quad a\le \frac{\left\lfloor (N/(2b))^{1/3}\right\rfloor-\frac{b}{2}}{2}.$$
These cube-root bounds drastically reduce the search space before the remaining filters are applied.
Worked Example: \(H(1000)=2535\)
The implementations include a small checkpoint at \(N=1000\). Here
$$\left\lfloor(4N)^{1/4}\right\rfloor=\left\lfloor 4000^{1/4}\right\rfloor=7.$$
Only two primitive seeds survive all conditions and still satisfy the size bound.
First, with \(a=1\) and odd \(b=1\),
$$X_\ast=1(1+4)^3=125,\qquad Y_\ast=4\cdot 1\cdot(1-2)^3=-4.$$
The kernel condition is \( \operatorname{cf}(1)\ne \operatorname{cf}(4)\), so the seed is valid. Its scale limit is
$$t_{\max}=\left\lfloor\sqrt{\frac{1000}{125}}\right\rfloor=2,$$
so its total contribution is
$$ (125+4)(1^2+2^2)=129\cdot 5=645. $$
Second, with \(a=1\) and even \(b=2\),
$$X_\ast=2\cdot 2\left(1+2\right)^3=108,\qquad Y_\ast=1(1-4)^3=-27.$$
Now \( \operatorname{cf}(4)\ne \operatorname{cf}(1)\), and
$$t_{\max}=\left\lfloor\sqrt{\frac{1000}{108}}\right\rfloor=3.$$
This contributes
$$ (108+27)(1^2+2^2+3^2)=135\cdot 14=1890. $$
Adding the two surviving families gives
$$645+1890=2535,$$
which matches the checkpoint in the implementations exactly.
How the Code Works
The implementation begins by computing exact integer fourth roots, cube roots, and square roots, so the search bounds are numerically safe even at very large \(N\).
It then builds a smallest-prime-factor sieve up to the limit needed for all cube-free-kernel lookups. From that sieve it derives the table of cube-free kernels \( \operatorname{cf}(n)\).
Next, it scans the outer parameter up to \( \lfloor(4N)^{1/4}\rfloor\). Odd values use the odd primitive formula; even values use the even primitive formula. For each branch it computes the cube-root upper bound for the inner parameter, then applies, in order, the coprimality test, the excluded diagonal test, the cube-free-kernel inequality, and the final max-norm check.
Whenever a primitive seed survives, the implementation computes \(t_{\max}\), evaluates the closed form for \(1^2+\cdots+t_{\max}^2\), multiplies by \(|X_\ast|+|Y_\ast|\), and adds the result modulo \(M\).
The C++ version parallelizes the outer-parameter loop across several threads, while the Python and Java implementations perform the same arithmetic sequentially. In all three languages, the mathematical work is the same.
Complexity Analysis
Let \(B=\lfloor(4N)^{1/4}\rfloor\). The outer loop runs over \(O(B)=O(N^{1/4})\) values.
For a fixed outer parameter \(b\), the inner limit is \(O((N/b)^{1/3})\). Summing this over all \(b\le B\) gives the rough bound
$$\sum_{b\le B} O\!\left((N/b)^{1/3}\right)=O\!\left(N^{1/3}\sum_{b\le B} b^{-1/3}\right)=O(N^{1/2}).$$
So the arithmetic enumeration is about \(O(N^{1/2})\) candidate checks before constant-factor pruning from coprimality and kernel filters. The sieve for cube-free kernels reaches only the largest parameter size, which is \(O(N^{1/3})\), so its cost is \(O(N^{1/3}\log\log N)\) time and \(O(N^{1/3})\) memory.
In practice, the wall-clock time is better than a naive pair scan by many orders of magnitude, and the threaded C++ implementation further reduces elapsed time without changing the asymptotic count.
Footnotes and References
- Problem page: https://projecteuler.net/problem=880
- Nested radicals: Wikipedia — Nested radical
- p-adic valuation: Wikipedia — p-adic valuation
- Power-free integers and cube-free structure: Wikipedia — Power-free integer
- Sum of squares: Wikipedia — Square pyramidal number
Problem 880 source code
C++
#include <algorithm>
#include <atomic>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <numeric>
#include <string>
#include <thread>
#include <vector>
using namespace std;
static constexpr uint64_t MOD = 1031ULL * 1031ULL * 1031ULL + 2ULL;
static inline __int128 i128_abs(__int128 x) { return x < 0 ? -x : x; }
static uint64_t isqrt_u64(uint64_t n) {
long double xd = sqrt((long double)n);
uint64_t x = (uint64_t)xd;
while ((__int128)(x + 1) * (x + 1) <= n) ++x;
while ((__int128)x * x > n) --x;
return x;
}
static uint64_t icbrt_u64(uint64_t n) {
long double xd = cbrt((long double)n);
uint64_t x = (uint64_t)xd;
auto cube = [](__int128 a) { return a * a * a; };
while (cube((__int128)(x + 1)) <= n) ++x;
while (cube((__int128)x) > n) --x;
return x;
}
static uint64_t i4rt_u64(uint64_t n) {
long double xd = pow((long double)n, 0.25L);
uint64_t x = (uint64_t)xd;
auto pow4 = [](__int128 a) { return a * a * a * a; };
while (pow4((__int128)(x + 1)) <= n) ++x;
while (pow4((__int128)x) > n) --x;
return x;
}
static vector<uint32_t> precompute_cube_free(uint32_t limit) {
vector<uint32_t> spf(limit + 1);
for (uint32_t i = 0; i <= limit; ++i) spf[i] = i;
for (uint32_t i = 2; (uint64_t)i * i <= limit; ++i) {
if (spf[i] == i) {
for (uint64_t j = (uint64_t)i * i; j <= limit; j += i) {
if (spf[(uint32_t)j] == (uint32_t)j) spf[(uint32_t)j] = i;
}
}
}
vector<uint32_t> cf(limit + 1);
cf[0] = 0;
if (limit >= 1) cf[1] = 1;
for (uint32_t n = 2; n <= limit; ++n) {
uint32_t x = n;
uint64_t res = 1;
while (x > 1) {
uint32_t p = spf[x];
uint32_t e = 0;
while (x % p == 0) {
x /= p;
++e;
}
uint32_t r = e % 3;
if (r == 1) res *= p;
else if (r == 2) res *= (uint64_t)p * p;
}
cf[n] = (uint32_t)res;
}
return cf;
}
static inline uint64_t sumsq_mod(uint64_t n) {
// n(n+1)(2n+1)/6
__int128 t = (__int128)n * (n + 1) * (2 * (__int128)n + 1);
t /= 6;
t %= MOD;
return (uint64_t)t;
}
struct PairLL {
long long x;
long long y;
};
static PairLL example_pair_from_pqm(uint64_t p, uint64_t q, uint64_t m) {
__int128 X0 = 0;
__int128 Y0 = 0;
if (q & 1ULL) {
__int128 t = (__int128)q + 4 * (__int128)p;
X0 = (__int128)q * t * t * t;
__int128 u = (__int128)p - 2 * (__int128)q;
Y0 = 4 * (__int128)p * u * u * u;
} else {
__int128 halfq = (__int128)q / 2;
__int128 t = halfq + 2 * (__int128)p;
X0 = 2 * (__int128)q * t * t * t;
__int128 u = (__int128)p - 2 * (__int128)q;
Y0 = (__int128)p * u * u * u;
}
__int128 mm = (__int128)m * m;
__int128 X = X0 * mm;
__int128 Y = Y0 * mm;
if (i128_abs(X) > i128_abs(Y)) swap(X, Y);
return {(long long)X, (long long)Y};
}
static uint64_t solve_H_mod(uint64_t N, int threads) {
uint64_t q_max = i4rt_u64((uint64_t)(4 * (__int128)N));
uint64_t max_p_even2 = 0;
if (N >= 4) {
uint64_t c = icbrt_u64(N / 4);
if (c > 1) max_p_even2 = (c - 1) / 2;
}
uint64_t max_p_odd1 = 0;
{
uint64_t c = icbrt_u64(N);
if (c > 1) max_p_odd1 = (c - 1) / 4;
}
uint64_t max_p = max(max_p_even2, max_p_odd1);
uint32_t limit = (uint32_t)(max<uint64_t>(4 * max_p, 2 * q_max) + 16);
vector<uint32_t> cubeFree = precompute_cube_free(limit);
threads = max(1, threads);
threads = min(threads, 64);
atomic<uint32_t> nextQ(1);
vector<uint64_t> partial(threads, 0);
auto worker = [&](int tid) {
uint64_t local = 0;
while (true) {
uint32_t q = nextQ.fetch_add(1, memory_order_relaxed);
if (q > q_max) break;
if (q & 1U) {
uint64_t Ndiv = N / q;
uint64_t c = icbrt_u64(Ndiv);
if (c <= q) continue;
uint64_t p_max = (c - q) / 4;
if (p_max == 0) continue;
uint32_t dx = cubeFree[q];
for (uint64_t p = 1; p <= p_max; ++p) {
if (std::gcd(p, (uint64_t)q) != 1) continue;
if (p == 2 * (uint64_t)q) continue;
uint64_t t = (uint64_t)q + 4ULL * p;
__int128 t3 = (__int128)t * t * t;
__int128 X0 = (__int128)q * t3;
__int128 u = (__int128)p - 2 * (__int128)q;
__int128 Y0 = 4 * (__int128)p * u * u * u;
__int128 absY = i128_abs(Y0);
__int128 max0 = (X0 > absY) ? X0 : absY;
if (max0 > (__int128)N) continue;
uint32_t dy = cubeFree[4 * p];
if (dx == dy) continue;
uint64_t max0_u = (uint64_t)max0;
uint64_t mmax = isqrt_u64(N / max0_u);
uint64_t s2 = sumsq_mod(mmax);
uint64_t absSum = ((uint64_t)X0 + (uint64_t)absY) % MOD;
local = (local + (uint64_t)((__int128)absSum * s2 % MOD)) % MOD;
}
} else {
uint64_t Ndiv = N / (2ULL * q);
uint64_t c = icbrt_u64(Ndiv);
uint64_t halfq = q / 2ULL;
if (c <= halfq) continue;
uint64_t p_max = (c - halfq) / 2;
if (p_max == 0) continue;
uint32_t dx = cubeFree[2 * q];
for (uint64_t p = 1; p <= p_max; ++p) {
if (std::gcd(p, (uint64_t)q) != 1) continue;
if (p == 2 * (uint64_t)q) continue;
uint64_t t = halfq + 2ULL * p;
__int128 t3 = (__int128)t * t * t;
__int128 X0 = 2 * (__int128)q * t3;
__int128 u = (__int128)p - 2 * (__int128)q;
__int128 Y0 = (__int128)p * u * u * u;
__int128 absY = i128_abs(Y0);
__int128 max0 = (X0 > absY) ? X0 : absY;
if (max0 > (__int128)N) continue;
uint32_t dy = cubeFree[p];
if (dx == dy) continue;
uint64_t max0_u = (uint64_t)max0;
uint64_t mmax = isqrt_u64(N / max0_u);
uint64_t s2 = sumsq_mod(mmax);
uint64_t absSum = ((uint64_t)X0 + (uint64_t)absY) % MOD;
local = (local + (uint64_t)((__int128)absSum * s2 % MOD)) % MOD;
}
}
}
partial[tid] = local;
};
vector<thread> pool;
pool.reserve(threads);
for (int t = 0; t < threads; ++t) pool.emplace_back(worker, t);
for (auto& th : pool) th.join();
uint64_t ans = 0;
for (uint64_t v : partial) ans = (ans + v) % MOD;
return ans;
}
static void run_selftest() {
{
PairLL p1 = example_pair_from_pqm(1, 1, 1);
assert(p1.x == -4 && p1.y == 125);
}
{
PairLL p2 = example_pair_from_pqm(5, 2, 1);
assert(p2.x == 5 && p2.y == 5324);
}
{
uint64_t got = solve_H_mod(1000, 1);
assert(got == 2535);
}
}
int main(int argc, char** argv) {
uint64_t N = 1000000000000000ULL;
int threads = (int)std::thread::hardware_concurrency();
if (threads <= 0) threads = 4;
bool selftest = false;
for (int i = 1; i < argc; ++i) {
string s = argv[i];
if (s == "--selftest") {
selftest = true;
} else if (s.rfind("--threads=", 0) == 0) {
threads = stoi(s.substr(10));
threads = max(1, threads);
} else {
N = stoull(s);
}
}
if (selftest) run_selftest();
uint64_t ans = solve_H_mod(N, threads);
cout << ans << "\n";
return 0;
}
Python
import math
import sys
MOD = 1031 * 1031 * 1031 + 2
def isqrt_u64(n):
return math.isqrt(n)
def icbrt_u64(n):
if n < 0:
return -icbrt_u64(-n)
x = int(math.floor(n ** (1.0/3.0)))
while (x + 1) ** 3 <= n:
x += 1
while x ** 3 > n:
x -= 1
return x
def i4rt_u64(n):
x = int(math.floor(n ** 0.25))
while (x + 1) ** 4 <= n:
x += 1
while x ** 4 > n:
x -= 1
return x
def precompute_cube_free(limit):
spf = list(range(limit + 1))
for i in range(2, math.isqrt(limit) + 1):
if spf[i] == i:
for j in range(i * i, limit + 1, i):
if spf[j] == j:
spf[j] = i
cf = [0] * (limit + 1)
if limit >= 1:
cf[1] = 1
for n in range(2, limit + 1):
x = n
res = 1
while x > 1:
p = spf[x]
e = 0
while x % p == 0:
x //= p
e += 1
r = e % 3
if r == 1:
res *= p
elif r == 2:
res *= p * p
cf[n] = res
return cf
def sumsq_mod(n):
t = n * (n + 1) * (2 * n + 1) // 6
return t % MOD
def solve_H_mod(N):
q_max = i4rt_u64(4 * N)
max_p_even2 = 0
if N >= 4:
c = icbrt_u64(N // 4)
if c > 1:
max_p_even2 = (c - 1) // 2
max_p_odd1 = 0
c = icbrt_u64(N)
if c > 1:
max_p_odd1 = (c - 1) // 4
max_p = max(max_p_even2, max_p_odd1)
limit = max(4 * max_p, 2 * q_max) + 16
cube_free = precompute_cube_free(limit)
ans = 0
for q in range(1, q_max + 1):
if q % 2 != 0:
Ndiv = N // q
c = icbrt_u64(Ndiv)
if c <= q:
continue
p_max = (c - q) // 4
if p_max == 0:
continue
dx = cube_free[q]
for p in range(1, p_max + 1):
if math.gcd(p, q) != 1:
continue
if p == 2 * q:
continue
t = q + 4 * p
t3 = t * t * t
X0 = q * t3
u = p - 2 * q
Y0 = 4 * p * u * u * u
absY = abs(Y0)
max0 = X0 if X0 > absY else absY
if max0 > N:
continue
dy = cube_free[4 * p]
if dx == dy:
continue
mmax = isqrt_u64(N // max0)
s2 = sumsq_mod(mmax)
absSum = (X0 + absY) % MOD
ans = (ans + absSum * s2) % MOD
else:
Ndiv = N // (2 * q)
c = icbrt_u64(Ndiv)
halfq = q // 2
if c <= halfq:
continue
p_max = (c - halfq) // 2
if p_max == 0:
continue
dx = cube_free[2 * q]
for p in range(1, p_max + 1):
if math.gcd(p, q) != 1:
continue
if p == 2 * q:
continue
t = halfq + 2 * p
t3 = t * t * t
X0 = 2 * q * t3
u = p - 2 * q
Y0 = p * u * u * u
absY = abs(Y0)
max0 = X0 if X0 > absY else absY
if max0 > N:
continue
dy = cube_free[p]
if dx == dy:
continue
mmax = isqrt_u64(N // max0)
s2 = sumsq_mod(mmax)
absSum = (X0 + absY) % MOD
ans = (ans + absSum * s2) % MOD
return ans % MOD
def solve():
return str(solve_H_mod(10**15))
if __name__ == "__main__":
print(solve())
Java
import java.math.BigInteger;
public class Euler880 {
static final long MOD = 1031L * 1031L * 1031L + 2L;
static long isqrtU64(long n) {
long x = (long) Math.sqrt(n);
while (new java.math.BigInteger(String.valueOf(x + 1)).pow(2)
.compareTo(new java.math.BigInteger(String.valueOf(n))) <= 0)
++x;
while (new java.math.BigInteger(String.valueOf(x)).pow(2)
.compareTo(new java.math.BigInteger(String.valueOf(n))) > 0)
--x;
return x;
}
static long icbrtU64(long n) {
long x = (long) Math.cbrt(n);
while (new java.math.BigInteger(String.valueOf(x + 1)).pow(3)
.compareTo(new java.math.BigInteger(String.valueOf(n))) <= 0)
++x;
while (new java.math.BigInteger(String.valueOf(x)).pow(3)
.compareTo(new java.math.BigInteger(String.valueOf(n))) > 0)
--x;
return x;
}
static long i4rtU64(long n) {
long x = (long) Math.pow(n, 0.25);
while (new java.math.BigInteger(String.valueOf(x + 1)).pow(4)
.compareTo(new java.math.BigInteger(String.valueOf(n))) <= 0)
++x;
while (new java.math.BigInteger(String.valueOf(x)).pow(4)
.compareTo(new java.math.BigInteger(String.valueOf(n))) > 0)
--x;
return x;
}
static int[] precomputeCubeFree(int limit) {
int[] spf = new int[limit + 1];
for (int i = 0; i <= limit; ++i)
spf[i] = i;
for (int i = 2; (long) i * i <= limit; ++i) {
if (spf[i] == i) {
for (long j = (long) i * i; j <= limit; j += i) {
if (spf[(int) j] == (int) j)
spf[(int) j] = i;
}
}
}
int[] cf = new int[limit + 1];
if (limit >= 1)
cf[1] = 1;
for (int n = 2; n <= limit; ++n) {
int x = n;
long res = 1;
while (x > 1) {
int p = spf[x];
int e = 0;
while (x % p == 0) {
x /= p;
++e;
}
int r = e % 3;
if (r == 1)
res *= p;
else if (r == 2)
res *= (long) p * p;
}
cf[n] = (int) res;
}
return cf;
}
static long sumsqMod(long n) {
BigInteger bn = BigInteger.valueOf(n);
BigInteger s = bn.multiply(bn.add(BigInteger.ONE))
.multiply(bn.multiply(BigInteger.valueOf(2)).add(BigInteger.ONE));
s = s.divide(BigInteger.valueOf(6));
return s.mod(BigInteger.valueOf(MOD)).longValue();
}
static long gcd(long a, long b) {
if (b == 0)
return a;
return gcd(b, a % b);
}
// I will use BigInteger to translate __int128 ops since intermediate values top
// 128 bit.
static long solveHMod(long N) {
long qMax = i4rtU64(4L * N);
long maxPEven2 = 0;
if (N >= 4) {
long c = icbrtU64(N / 4);
if (c > 1)
maxPEven2 = (c - 1) / 2;
}
long maxPOdd1 = 0;
long c_odd = icbrtU64(N);
if (c_odd > 1)
maxPOdd1 = (c_odd - 1) / 4;
long maxP = Math.max(maxPEven2, maxPOdd1);
int limit = (int) (Math.max(4L * maxP, 2L * qMax) + 16);
int[] cubeFree = precomputeCubeFree(limit);
long ans = 0;
BigInteger bN = BigInteger.valueOf(N);
for (long q = 1; q <= qMax; ++q) {
if (q % 2 != 0) {
long ndiv = N / q;
long c = icbrtU64(ndiv);
if (c <= q)
continue;
long pMax = (c - q) / 4;
if (pMax == 0)
continue;
int dx = cubeFree[(int) q];
for (long p = 1; p <= pMax; ++p) {
if (gcd(p, q) != 1)
continue;
if (p == 2L * q)
continue;
long t = q + 4L * p;
BigInteger bt = BigInteger.valueOf(t);
BigInteger bt3 = bt.pow(3);
BigInteger x0 = BigInteger.valueOf(q).multiply(bt3);
long u = p - 2L * q;
BigInteger bu = BigInteger.valueOf(u);
BigInteger bu3 = bu.pow(3);
BigInteger y0 = BigInteger.valueOf(4L * p).multiply(bu3);
BigInteger absY = y0.abs();
BigInteger max0 = x0.compareTo(absY) > 0 ? x0 : absY;
if (max0.compareTo(bN) > 0)
continue;
int dy = cubeFree[(int) (4 * p)];
if (dx == dy)
continue;
long mmax = isqrtU64(bN.divide(max0).longValue());
long s2 = sumsqMod(mmax);
long absSum = x0.add(absY).mod(BigInteger.valueOf(MOD)).longValue();
ans = (ans + (absSum * s2) % MOD) % MOD;
}
} else {
long ndiv = N / (2L * q);
long c = icbrtU64(ndiv);
long halfq = q / 2L;
if (c <= halfq)
continue;
long pMax = (c - halfq) / 2;
if (pMax == 0)
continue;
int dx = cubeFree[(int) (2 * q)];
for (long p = 1; p <= pMax; ++p) {
if (gcd(p, q) != 1)
continue;
if (p == 2L * q)
continue;
long t = halfq + 2L * p;
BigInteger bt = BigInteger.valueOf(t);
BigInteger bt3 = bt.pow(3);
BigInteger x0 = BigInteger.valueOf(2L * q).multiply(bt3);
long u = p - 2L * q;
BigInteger bu = BigInteger.valueOf(u);
BigInteger bu3 = bu.pow(3);
BigInteger y0 = BigInteger.valueOf(p).multiply(bu3);
BigInteger absY = y0.abs();
BigInteger max0 = x0.compareTo(absY) > 0 ? x0 : absY;
if (max0.compareTo(bN) > 0)
continue;
int dy = cubeFree[(int) p];
if (dx == dy)
continue;
long mmax = isqrtU64(bN.divide(max0).longValue());
long s2 = sumsqMod(mmax);
long absSum = x0.add(absY).mod(BigInteger.valueOf(MOD)).longValue();
ans = (ans + (absSum * s2) % MOD) % MOD;
}
}
}
return ans;
}
public static String solve() {
return Long.toString(solveHMod(1000000000000000L));
}
public static void main(String[] args) {
System.out.println(solve());
}
}