Problem 715: Sextuplet Norms

Project Euler Problem 715 Solution

EulerSolve provides an optimized solution for Project Euler Problem 715, Sextuplet Norms, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The target quantity is a summatory arithmetic function $G(n)$, evaluated modulo $10^9+7$, with the Project Euler input at $n=10^{12}$. The native implementations show that the summand is multiplicative and that its local behavior depends on whether an odd prime is congruent to $1$ or $3$ modulo $4$. The relevant Dirichlet character is $$\chi_4(m)=\begin{cases} 0,& 2\mid m,\\ 1,& m\equiv 1\pmod 4,\\ -1,& m\equiv 3\pmod 4. \end{cases}$$ Reading the prime and prime-power contributions from the implementations gives the multiplicative term $$a(1)=1,\qquad a(p^e)=p^{3e}\left(1-\frac{\chi_4(p)}{p^3}\right)=p^{3e}-\chi_4(p)p^{3e-3}.$$ Therefore the problem is to compute $$G(n)=\sum_{m\le n} a(m)\pmod{10^9+7}.$$ A direct loop to $10^{12}$ is not viable, so the implementation uses a compressed prime-prefix sieve together with a DFS over admissible prime-power products. Mathematical Approach The entire method can be reconstructed from two facts visible in the implementations: a prime first appears with weight $p^3-\chi_4(p)$, and every extra copy of the same prime multiplies the contribution by another factor $p^3$....

Detailed mathematical approach

Problem Summary

The target quantity is a summatory arithmetic function $G(n)$, evaluated modulo $10^9+7$, with the Project Euler input at $n=10^{12}$. The native implementations show that the summand is multiplicative and that its local behavior depends on whether an odd prime is congruent to $1$ or $3$ modulo $4$.

The relevant Dirichlet character is

$$\chi_4(m)=\begin{cases} 0,& 2\mid m,\\ 1,& m\equiv 1\pmod 4,\\ -1,& m\equiv 3\pmod 4. \end{cases}$$

Reading the prime and prime-power contributions from the implementations gives the multiplicative term

$$a(1)=1,\qquad a(p^e)=p^{3e}\left(1-\frac{\chi_4(p)}{p^3}\right)=p^{3e}-\chi_4(p)p^{3e-3}.$$

Therefore the problem is to compute

$$G(n)=\sum_{m\le n} a(m)\pmod{10^9+7}.$$

A direct loop to $10^{12}$ is not viable, so the implementation uses a compressed prime-prefix sieve together with a DFS over admissible prime-power products.

Mathematical Approach

The entire method can be reconstructed from two facts visible in the implementations: a prime first appears with weight $p^3-\chi_4(p)$, and every extra copy of the same prime multiplies the contribution by another factor $p^3$.

Step 1: Identify the Local Prime Weight

For a prime $p$, the first occurrence contributes

$$w(p)=p^3-\chi_4(p).$$

So the three residue classes behave as follows:

$$w(2)=8,\qquad w(p)=p^3-1\ \text{if }p\equiv 1\pmod 4,\qquad w(p)=p^3+1\ \text{if }p\equiv 3\pmod 4.$$

This is the arithmetic signature of the problem: primes $1\bmod 4$ and $3\bmod 4$ are treated differently, exactly through the character $\chi_4$.

Step 2: Extend from Primes to Prime Powers

Once a branch already contains $p$, increasing the exponent by $1$ multiplies the current local contribution by $p^3$. Therefore

$$a(p^e)=w(p)\,p^{3(e-1)}=(p^3-\chi_4(p))p^{3(e-1)}=p^{3e}-\chi_4(p)p^{3e-3}.$$

In particular, the factor for $p=2$ is simply $2^{3e}$, because $\chi_4(2)=0$. For odd primes the correction term is exactly one unit of $p^{3e-3}$, with sign determined by $p\bmod 4$.

Step 3: Rebuild the Global Multiplicative Function

Independent prime branches multiply, so for

$$n=\prod_{p^e\parallel n} p^e$$

we get

$$a(n)=\prod_{p^e\parallel n} a(p^e)=n^3\prod_{p\mid n}\left(1-\frac{\chi_4(p)}{p^3}\right).$$

This is the closed form encoded by the DFS. An equivalent divisor-sum identity is

$$a(n)=\sum_{d\mid n}\mu(d)\chi_4(d)\left(\frac{n}{d}\right)^3,$$

because only squarefree divisors survive the Möbius factor. That identity is not used directly in the code, but it confirms the prime-power formula and makes the multiplicativity transparent.

Step 4: Separate the Prime Contribution

The summatory function can be decomposed into the contribution of $1$, the contribution of primes, and the remaining composite branches. The prime part is

$$\sum_{p\le x} a(p)=\sum_{p\le x}\bigl(p^3-\chi_4(p)\bigr).$$

For that reason the implementation builds two prime-prefix tables:

$$S_3(x)=\sum_{p\le x} p^3,\qquad C_4(x)=\sum_{p\le x}\chi_4(p).$$

Then the prime contribution is their difference

$$S(x)=S_3(x)-C_4(x)=\sum_{p\le x}\bigl(p^3-\chi_4(p)\bigr).$$

Both tables start from easy full-integer prefixes and then sieve away composite contributions prime by prime. For cubes, the starting formula is

$$\sum_{m=2}^{x} m^3=\left(\frac{x(x+1)}{2}\right)^2-1.$$

For the mod-$4$ character, the full-integer prefix is periodic:

$$\sum_{m=2}^{x}\chi_4(m)=\begin{cases} 0,& x\equiv 1,2\pmod 4,\\ -1,& x\equiv 0,3\pmod 4. \end{cases}$$

Step 5: Use Floor-Division Compression and DFS

The values $\left\lfloor n/i\right\rfloor$ repeat heavily, so the implementation stores prime-prefix information only on the compressed domain

$$\mathcal{D}(n)=\left\{1,2,\dots,\left\lfloor\frac{n}{\lfloor\sqrt n\rfloor}\right\rfloor\right\}\cup\left\{\left\lfloor\frac{n}{i}\right\rfloor:1\le i\le \lfloor\sqrt n\rfloor\right\}.$$

This reduces the table size from $O(n)$ to $O(\sqrt n)$. After the prime-only prefixes are available on that domain, a DFS enumerates multiplicative branches in increasing prime order. Ordering the primes this way prevents double-counting. At each branch:

$$\text{first copy of }p \Rightarrow p^3-\chi_4(p),\qquad \text{each extra copy of }p \Rightarrow \times p^3.$$

The prime-prefix tables instantly add the cases where a branch is extended by one larger prime, while recursion handles products containing several larger primes.

Worked Example: $G(10)$

The local formula gives

$$\begin{aligned} a(1)&=1,\\ a(2)&=2^3=8,\\ a(3)&=3^3\left(1+\frac{1}{3^3}\right)=28,\\ a(4)&=4^3=64,\\ a(5)&=5^3\left(1-\frac{1}{5^3}\right)=124,\\ a(6)&=a(2)a(3)=224,\\ a(7)&=7^3\left(1+\frac{1}{7^3}\right)=344,\\ a(8)&=8^3=512,\\ a(9)&=9^3\left(1+\frac{1}{3^3}\right)=756,\\ a(10)&=a(2)a(5)=992. \end{aligned}$$

Summing these values yields

$$G(10)=1+8+28+64+124+224+344+512+756+992=3053,$$

which matches the checkpoint used by the native implementations.

How the Code Works

The C++ and Java implementations first generate all primes up to $\lfloor\sqrt n\rfloor$. They then build two compressed prime-prefix tables: one for $\sum_{p\le x} p^3$, and one for $\sum_{p\le x}\chi_4(p)$. Each table starts from the corresponding prefix over all integers and applies a prime-sieving transform so that only prime contributions remain.

After subtracting those tables, the implementation has fast access to $\sum_{p\le x}(p^3-\chi_4(p))$ for every needed compressed argument $x$. A DFS then enumerates composite multiplicative branches in increasing prime order. When a new prime enters a branch, the factor is $p^3-\chi_4(p)$; when the same prime is repeated, the factor gains another multiplier $p^3$. The prefix tables let the implementation add all one-more-prime extensions in bulk, while recursion handles deeper products.

The final answer is

$$G(n)=1+\sum_{p\le n}\bigl(p^3-\chi_4(p)\bigr)+\text{composite DFS contribution}\pmod{10^9+7}.$$

The Python implementation is intentionally thin: it compiles and launches the native solver when needed, then returns the parsed numeric result.

Complexity Analysis

The compressed floor-division domain has size $O(\sqrt n)$, so the memory usage is $O(\sqrt n)$. The two prime-prefix transforms run on that compressed domain rather than on all integers up to $n$; in standard Min_25-style analysis this is roughly $O(n^{3/4}/\log n)$ work. The DFS over admissible prime-power products is much smaller than a direct scan to $n$ and is practical for the required input size. Overall the method is decisively sublinear in $n$ and is designed specifically for $n=10^{12}$.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=715
Dirichlet character: Wikipedia — Dirichlet character
Multiplicative function: Wikipedia — Multiplicative function
Möbius function: Wikipedia — Möbius function
Min_25 sieve overview: OI Wiki — Min_25 sieve

Problem 715 source code

C++

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

namespace {

using i64 = long long;
using i128 = __int128_t;
using u64 = std::uint64_t;

constexpr i64 kMod = 1'000'000'007LL;
constexpr i64 kDefaultN = 1'000'000'000'000LL;
constexpr i64 kCheckpointN1 = 10LL;
constexpr i64 kCheckpointExpected1 = 3'053LL;
constexpr i64 kCheckpointN2 = 100'000LL;
constexpr i64 kCheckpointExpected2 = 157'612'967LL;

struct Options {
    i64 n = kDefaultN;
    bool run_checkpoints = true;
};

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, i64& value) {
    if (arg.rfind(prefix, 0U) != 0U) return false;
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) return false;

    unsigned long long parsed = 0ULL;
    for (const char ch : tail) {
        if (ch < '0' || ch > '9') return false;
        const unsigned long long digit = static_cast<unsigned long long>(ch - '0');
        if (parsed > (std::numeric_limits<unsigned long long>::max() - digit) / 10ULL) return false;
        parsed = parsed * 10ULL + digit;
    }
    if (parsed > static_cast<unsigned long long>(std::numeric_limits<i64>::max())) return false;
    value = static_cast<i64>(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-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--n=", options.n)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n >= 1;
}

inline i64 mod_norm(i64 x) {
    x %= kMod;
    if (x < 0) x += kMod;
    return x;
}

inline i64 mod_add(i64 a, i64 b) {
    i64 s = a + b;
    if (s >= kMod) s -= kMod;
    return s;
}

inline i64 mod_mul(i64 a, i64 b) {
    return static_cast<i64>((static_cast<u64>(a) * static_cast<u64>(b)) % static_cast<u64>(kMod));
}

std::vector<int> sieve_primes(const int n) {
    std::vector<char> is_prime(static_cast<std::size_t>(n + 1), 1);
    if (n >= 0) is_prime[0] = 0;
    if (n >= 1) is_prime[1] = 0;
    for (int i = 2; 1LL * i * i <= n; ++i) {
        if (!is_prime[static_cast<std::size_t>(i)]) continue;
        for (int j = i * i; j <= n; j += i) is_prime[static_cast<std::size_t>(j)] = 0;
    }

    std::vector<int> primes;
    for (int i = 2; i <= n; ++i) {
        if (is_prime[static_cast<std::size_t>(i)]) primes.push_back(i);
    }
    return primes;
}

std::pair<std::vector<i64>, std::vector<i64>> sum_prime_cubes(
    const i64 n,
    const int L,
    const int XL,
    const std::vector<int>& primes
) {
    std::vector<i64> V(static_cast<std::size_t>(XL), 0);
    i64 s = -1;
    for (int i = 1; i < XL; ++i) {
        const i64 ii = i;
        s = mod_norm(s + mod_mul(mod_mul(ii % kMod, ii % kMod), ii % kMod));
        V[static_cast<std::size_t>(i)] = s;
    }

    std::vector<i64> bigV(static_cast<std::size_t>(L + 1), 0);
    const i64 mod2 = kMod * 2LL;
    for (int i = 1; i <= L; ++i) {
        const i64 x = (n / i) % mod2;
        const i64 y = static_cast<i64>((static_cast<i128>(x) * (x + 1LL) % mod2) / 2);
        bigV[static_cast<std::size_t>(i)] = mod_norm(mod_mul(y % kMod, y % kMod) - 1);
    }

    for (const int p : primes) {
        const i64 sp = V[static_cast<std::size_t>(p - 1)];
        const i64 p_sq = 1LL * p * p;
        const i64 p3 = mod_mul(p % kMod, mod_mul(p % kMod, p % kMod));
        const i64 y = n / p;
        const int iL = static_cast<int>(std::min<i64>(y / p, L));

        for (int i = 1; i <= iL; ++i) {
            const i64 z = y / i;
            const i64 v = (z < XL) ? V[static_cast<std::size_t>(z)] : bigV[static_cast<std::size_t>(i * p)];
            bigV[static_cast<std::size_t>(i)] = mod_norm(
                bigV[static_cast<std::size_t>(i)] - mod_mul(p3, mod_norm(v - sp)));
        }

        for (int x = XL - 1; x > 0; --x) {
            if (x < p_sq) break;
            const int z = x / p;
            V[static_cast<std::size_t>(x)] = mod_norm(
                V[static_cast<std::size_t>(x)] - mod_mul(p3, mod_norm(V[static_cast<std::size_t>(z)] - sp)));
        }
    }

    return {std::move(V), std::move(bigV)};
}

std::pair<std::vector<i64>, std::vector<i64>> prime_balance_mod4(
    const i64 n,
    const int L,
    const int XL,
    const std::vector<int>& primes
) {
    std::vector<i64> V(static_cast<std::size_t>(XL), 0);
    for (int i = 3; i < XL; i += 4) V[static_cast<std::size_t>(i)] = -1;
    for (int i = 4; i < XL; i += 4) V[static_cast<std::size_t>(i)] = -1;

    std::vector<i64> bigV(static_cast<std::size_t>(L + 1), 0);
    for (int i = 1; i <= L; ++i) {
        if (((n / i - 1) % 4) >= 2) bigV[static_cast<std::size_t>(i)] = -1;
    }

    for (std::size_t idx = 1; idx < primes.size(); ++idx) {
        const int p = primes[idx];
        const i64 sp = V[static_cast<std::size_t>(p - 1)];
        const i64 y = n / p;
        const int iL = static_cast<int>(std::min<i64>(y / p, L));
        const i64 f = 2 - (p % 4);

        for (int i = 1; i <= iL; ++i) {
            const i64 z = y / i;
            if (z < XL) {
                bigV[static_cast<std::size_t>(i)] -= f * (V[static_cast<std::size_t>(z)] - sp);
            } else {
                bigV[static_cast<std::size_t>(i)] -= f * (bigV[static_cast<std::size_t>(i * p)] - sp);
            }
        }

        const i64 p_sq = 1LL * p * p;
        for (int x = XL - 1; x > 0; --x) {
            if (x < p_sq) break;
            const int z = x / p;
            V[static_cast<std::size_t>(x)] -= f * (V[static_cast<std::size_t>(z)] - sp);
        }
    }

    return {std::move(V), std::move(bigV)};
}

struct DfsContext {
    const std::vector<int>& primes;
    const std::vector<i64>& p2;
    const std::vector<i64>& p3;
    int L = 0;
    const std::vector<i64>& V;
    const std::vector<i64>& bigV;

    i64 dfs_branch(const int i, const i64 n0, const i64 x0) const {
        const i64 p_sq = p2[static_cast<std::size_t>(i)];
        if (n0 < p_sq) return 0;

        const i64 p = primes[static_cast<std::size_t>(i)];
        const i64 p_cubed = p3[static_cast<std::size_t>(i)];
        i64 n = n0 / p;
        i64 x = x0 * p;
        i64 f = mod_norm(p_cubed + (p % 4) - 2);
        i64 res = 0;

        const i64 pref = (x <= L) ? bigV[static_cast<std::size_t>(x)] : V[static_cast<std::size_t>(n)];
        res = mod_add(res, mod_mul(f, mod_norm(pref - V[static_cast<std::size_t>(p)])));

        while (true) {
            if (n > p_sq) {
                res = mod_add(res, mod_mul(f, dfs(i + 1, n, x)));
            }
            n /= p;
            if (n == 0) break;

            x *= p;
            f = mod_mul(f, p_cubed);
            res = mod_add(res, f);

            if (n > p) {
                const i64 pref2 = (x <= L) ? bigV[static_cast<std::size_t>(x)] : V[static_cast<std::size_t>(n)];
                res = mod_add(res, mod_mul(f, mod_norm(pref2 - V[static_cast<std::size_t>(p)])));
            }
        }
        return res;
    }

    i64 dfs(const int i0, const i64 n0, const i64 x0) const {
        i64 res = 0;
        for (int i = i0; i < static_cast<int>(primes.size()); ++i) {
            if (n0 < p2[static_cast<std::size_t>(i)]) break;
            res = mod_add(res, dfs_branch(i, n0, x0));
        }
        return res;
    }
};

i64 solve(const i64 n) {
    const int L = static_cast<int>(std::sqrt(static_cast<long double>(n) + 0.5L));
    const int XL = static_cast<int>(n / L) + 1;
    const std::vector<int> primes = sieve_primes(L);

    auto fut_prime_balance = std::async(std::launch::async, [&]() {
        return prime_balance_mod4(n, L, XL, primes);
    });
    auto [V1, bigV1] = sum_prime_cubes(n, L, XL, primes);
    auto [V2, bigV2] = fut_prime_balance.get();

    std::vector<i64> V(static_cast<std::size_t>(XL), 0);
    for (int i = 0; i < XL; ++i) {
        V[static_cast<std::size_t>(i)] = mod_norm(
            V1[static_cast<std::size_t>(i)] - mod_norm(V2[static_cast<std::size_t>(i)]));
    }

    std::vector<i64> bigV(static_cast<std::size_t>(L + 1), 0);
    for (int i = 0; i <= L; ++i) {
        bigV[static_cast<std::size_t>(i)] = mod_norm(
            bigV1[static_cast<std::size_t>(i)] - mod_norm(bigV2[static_cast<std::size_t>(i)]));
    }

    std::vector<i64> p2(primes.size(), 0);
    std::vector<i64> p3(primes.size(), 0);
    for (std::size_t i = 0; i < primes.size(); ++i) {
        const i64 p = primes[i];
        p2[i] = p * p;
        p3[i] = mod_mul(p % kMod, mod_mul(p % kMod, p % kMod));
    }

    const DfsContext ctx{primes, p2, p3, L, V, bigV};

    int root_end = 0;
    while (root_end < static_cast<int>(primes.size()) && p2[static_cast<std::size_t>(root_end)] <= n) {
        ++root_end;
    }

    if (root_end == 0) {
        return mod_norm(1 + bigV[1]);
    }

    unsigned int threads = std::thread::hardware_concurrency();
    if (threads == 0) threads = 1;
    threads = std::min<unsigned int>(threads, static_cast<unsigned int>(root_end));
    if (threads == 1) {
        return mod_norm(1 + bigV[1] + ctx.dfs(0, n, 1));
    }

    std::atomic<int> next_idx(0);
    std::vector<i64> partial(static_cast<std::size_t>(threads), 0);
    std::vector<std::thread> workers;
    workers.reserve(static_cast<std::size_t>(threads));

    for (unsigned int tid = 0; tid < threads; ++tid) {
        workers.emplace_back([&, tid]() {
            i64 local = 0;
            while (true) {
                const int i = next_idx.fetch_add(1, std::memory_order_relaxed);
                if (i >= root_end) break;
                local = mod_add(local, ctx.dfs_branch(i, n, 1));
            }
            partial[static_cast<std::size_t>(tid)] = local;
        });
    }
    for (auto& th : workers) th.join();

    i64 dfs_total = 0;
    for (const i64 v : partial) dfs_total = mod_add(dfs_total, v);
    return mod_norm(1 + bigV[1] + dfs_total);
}

bool run_checkpoints() {
    if (solve(kCheckpointN1) != kCheckpointExpected1) {
        std::cerr << "Validation failed: G(10)\n";
        return false;
    }
    if (solve(kCheckpointN2) != kCheckpointExpected2) {
        std::cerr << "Validation failed: G(10^5)\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_checkpoints && !run_checkpoints()) {
        return 1;
    }
    std::cout << solve(options.n) << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.util.ArrayList;
import java.util.List;
import java.util.Arrays;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import java.util.concurrent.Future;
import java.util.concurrent.Callable;

public class Euler715 {
    static final long kMod = 1000000007L;

    static long modNorm(long x) {
        x %= kMod;
        if (x < 0)
            x += kMod;
        return x;
    }

    static long modAdd(long a, long b) {
        long s = a + b;
        if (s >= kMod)
            s -= kMod;
        return s;
    }

    static long modMul(long a, long b) {
        return (a * b) % kMod;
    }

    static int[] sievePrimes(int n) {
        byte[] isPrime = new byte[n + 1];
        Arrays.fill(isPrime, (byte) 1);
        if (n >= 0)
            isPrime[0] = 0;
        if (n >= 1)
            isPrime[1] = 0;

        for (int i = 2; (long) i * i <= n; ++i) {
            if (isPrime[i] != 0) {
                for (int j = i * i; j <= n; j += i) {
                    isPrime[j] = 0;
                }
            }
        }

        int count = 0;
        for (int i = 2; i <= n; ++i) {
            if (isPrime[i] != 0)
                count++;
        }

        int[] primes = new int[count];
        int idx = 0;
        for (int i = 2; i <= n; ++i) {
            if (isPrime[i] != 0)
                primes[idx++] = i;
        }
        return primes;
    }

    static class Pair<T, U> {
        T first;
        U second;

        Pair(T first, U second) {
            this.first = first;
            this.second = second;
        }
    }

    static Pair<long[], long[]> sumPrimeCubes(long n, int L, int XL, int[] primes) {
        long[] V = new long[XL];
        long s = -1;
        for (int i = 1; i < XL; ++i) {
            long ii = i % kMod;
            s = modNorm(s + modMul(modMul(ii, ii), ii));
            V[i] = s;
        }

        long[] bigV = new long[L + 1];
        long mod2 = kMod * 2L;
        for (int i = 1; i <= L; ++i) {
            long x = (n / i) % mod2;
            long y = (x * (x + 1L) % mod2) / 2L;
            bigV[i] = modNorm(modMul(y % kMod, y % kMod) - 1);
        }

        for (int p : primes) {
            long sp = V[p - 1];
            long pSq = (long) p * p;
            long pMod = p % kMod;
            long p3 = modMul(pMod, modMul(pMod, pMod));
            long y = n / p;
            int iL = (int) Math.min(y / p, L);

            for (int i = 1; i <= iL; ++i) {
                long z = y / i;
                long v = (z < XL) ? V[(int) z] : bigV[i * p];
                bigV[i] = modNorm(bigV[i] - modMul(p3, modNorm(v - sp)));
            }

            for (int x = XL - 1; x > 0; --x) {
                if (x < pSq)
                    break;
                int z = x / p;
                V[x] = modNorm(V[x] - modMul(p3, modNorm(V[z] - sp)));
            }
        }

        return new Pair<>(V, bigV);
    }

    static Pair<long[], long[]> primeBalanceMod4(long n, int L, int XL, int[] primes) {
        long[] V = new long[XL];
        for (int i = 3; i < XL; i += 4)
            V[i] = -1;
        for (int i = 4; i < XL; i += 4)
            V[i] = -1;

        long[] bigV = new long[L + 1];
        for (int i = 1; i <= L; ++i) {
            if (((n / i - 1) % 4) >= 2)
                bigV[i] = -1;
        }

        for (int p : primes) {
            long sp = V[p - 1];
            long y = n / p;
            int iL = (int) Math.min(y / p, L);
            long f = 2 - (p % 4);

            for (int i = 1; i <= iL; ++i) {
                long z = y / i;
                if (z < XL) {
                    bigV[i] -= f * (V[(int) z] - sp);
                } else {
                    bigV[i] -= f * (bigV[i * p] - sp);
                }
            }

            long pSq = (long) p * p;
            for (int x = XL - 1; x > 0; --x) {
                if (x < pSq)
                    break;
                int z = x / p;
                V[x] -= f * (V[z] - sp);
            }
        }

        return new Pair<>(V, bigV);
    }

    static class DfsContext {
        int[] primes;
        long[] p2;
        long[] p3;
        int L;
        long[] V;
        long[] bigV;

        long dfsBranch(int i, long n0, long x0) {
            long pSq = p2[i];
            if (n0 < pSq)
                return 0;

            long p = primes[i];
            long pCubed = p3[i];
            long n = n0 / p;
            long x = x0 * p;
            long f = modNorm(pCubed + (p % 4) - 2);
            long res = 0;

            long pref = (x <= L) ? bigV[(int) x] : V[(int) n];
            res = modAdd(res, modMul(f, modNorm(pref - V[(int) p])));

            while (true) {
                if (n > pSq) {
                    res = modAdd(res, modMul(f, dfs(i + 1, n, x)));
                }
                n /= p;
                if (n == 0)
                    break;

                x *= p;
                f = modMul(f, pCubed);
                res = modAdd(res, f);

                if (n > p) {
                    long pref2 = (x <= L) ? bigV[(int) x] : V[(int) n];
                    res = modAdd(res, modMul(f, modNorm(pref2 - V[(int) p])));
                }
            }
            return res;
        }

        long dfs(int i0, long n0, long x0) {
            long res = 0;
            for (int i = i0; i < primes.length; ++i) {
                if (n0 < p2[i])
                    break;
                res = modAdd(res, dfsBranch(i, n0, x0));
            }
            return res;
        }
    }

    public static String solve() {
        long n = 1000000000000L;
        int L = (int) Math.sqrt(n + 0.5);
        int XL = (int) (n / L) + 1;
        int[] primes = sievePrimes(L);

        Pair<long[], long[]> res1 = sumPrimeCubes(n, L, XL, primes);
        long[] V1 = res1.first;
        long[] bigV1 = res1.second;

        Pair<long[], long[]> res2 = primeBalanceMod4(n, L, XL, primes);
        long[] V2 = res2.first;
        long[] bigV2 = res2.second;

        long[] V = new long[XL];
        for (int i = 0; i < XL; ++i) {
            V[i] = modNorm(V1[i] - modNorm(V2[i]));
        }

        long[] bigV = new long[L + 1];
        for (int i = 0; i <= L; ++i) {
            bigV[i] = modNorm(bigV1[i] - modNorm(bigV2[i]));
        }

        long[] p2 = new long[primes.length];
        long[] p3 = new long[primes.length];
        for (int i = 0; i < primes.length; ++i) {
            long p = primes[i];
            p2[i] = p * p;
            long pMod = p % kMod;
            p3[i] = modMul(pMod, modMul(pMod, pMod));
        }

        DfsContext ctx = new DfsContext();
        ctx.primes = primes;
        ctx.p2 = p2;
        ctx.p3 = p3;
        ctx.L = L;
        ctx.V = V;
        ctx.bigV = bigV;

        int rootEnd = 0;
        while (rootEnd < primes.length && p2[rootEnd] <= n) {
            rootEnd++;
        }

        if (rootEnd == 0) {
            return Long.toString(modNorm(1 + bigV[1]));
        }

        int threads = Runtime.getRuntime().availableProcessors();
        threads = Math.min(threads, rootEnd);

        ExecutorService executor = Executors.newFixedThreadPool(threads);
        List<Future<Long>> futures = new ArrayList<>();

        for (int i = 0; i < rootEnd; i++) {
            final int branchI = i;
            futures.add(executor.submit(() -> ctx.dfsBranch(branchI, n, 1L)));
        }

        long dfsTotal = 0;
        try {
            for (Future<Long> f : futures) {
                dfsTotal = modAdd(dfsTotal, f.get());
            }
        } catch (Exception e) {
            e.printStackTrace();
        }
        executor.shutdown();

        return Long.toString(modNorm(1 + bigV[1] + dfsTotal));
    }

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