Problem 468: Smooth Divisors of Binomial Coefficients

Project Euler Problem 468 Solution

EulerSolve provides an optimized solution for Project Euler Problem 468, Smooth Divisors of Binomial Coefficients, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For an integer $x$ and a smoothness bound $b$, define the largest $b$-smooth divisor by $$S_b(x)=\prod_{p^\alpha \parallel x,\ p \le b} p^\alpha.$$ Only prime powers whose prime base is at most $b$ are kept. For example, $2100=2^2\cdot 3\cdot 5^2\cdot 7$, so $S_4(2100)=2^2\cdot 3=12$. Problem 468 asks for $$F(n)=\sum_{r=0}^{n}\sum_{b=1}^{n} S_b\!\left(\binom{n}{r}\right) \pmod{M},\qquad M=1{,}000{,}000{,}993.$$ A naive approach would refactor every binomial coefficient for every value of $b$. The successful method instead derives a closed form for one fixed coefficient and then updates that value incrementally as $r$ moves across the row. Mathematical Approach The solution has two layers. First, for one fixed integer $x$, it converts $\sum_b S_b(x)$ into a weighted sum of prime-prefix products. Second, it exploits the recurrence between consecutive binomial coefficients so that only a few prime exponents change from one step to the next. Step 1: Rewrite the total as row contributions For each $r$, write $$x_r=\binom{n}{r},\qquad T_r=\sum_{b=1}^{n} S_b(x_r).$$ Then the required answer is simply $$F(n)=\sum_{r=0}^{n} T_r \pmod{M}.$$ Every prime divisor of $\binom{n}{r}$ is at most $n$, so it is enough to consider the primes up to $n$....

Detailed mathematical approach

Problem Summary

For an integer $x$ and a smoothness bound $b$, define the largest $b$-smooth divisor by

$$S_b(x)=\prod_{p^\alpha \parallel x,\ p \le b} p^\alpha.$$

Only prime powers whose prime base is at most $b$ are kept. For example, $2100=2^2\cdot 3\cdot 5^2\cdot 7$, so $S_4(2100)=2^2\cdot 3=12$.

Problem 468 asks for

$$F(n)=\sum_{r=0}^{n}\sum_{b=1}^{n} S_b\!\left(\binom{n}{r}\right) \pmod{M},\qquad M=1{,}000{,}000{,}993.$$

A naive approach would refactor every binomial coefficient for every value of $b$. The successful method instead derives a closed form for one fixed coefficient and then updates that value incrementally as $r$ moves across the row.

Mathematical Approach

The solution has two layers. First, for one fixed integer $x$, it converts $\sum_b S_b(x)$ into a weighted sum of prime-prefix products. Second, it exploits the recurrence between consecutive binomial coefficients so that only a few prime exponents change from one step to the next.

Step 1: Rewrite the total as row contributions

For each $r$, write

$$x_r=\binom{n}{r},\qquad T_r=\sum_{b=1}^{n} S_b(x_r).$$

Then the required answer is simply

$$F(n)=\sum_{r=0}^{n} T_r \pmod{M}.$$

Every prime divisor of $\binom{n}{r}$ is at most $n$, so it is enough to consider the primes up to $n$.

Step 2: Turn the sum over $b$ into weighted prime gaps

Let the primes up to $n$ be

$$2=p_1<p_2<\cdots<p_t\le n,$$

and introduce the sentinel

$$p_{t+1}=n+1.$$

Write the factorization of the current binomial coefficient as

$$x_r=\prod_{i=1}^{t} p_i^{\alpha_i(r)},$$

where some exponents may be zero. If $p_i \le b < p_{i+1}$, then the allowed primes are exactly $p_1,\dots,p_i$, so

$$S_b(x_r)=\prod_{j=1}^{i} p_j^{\alpha_j(r)}.$$

Therefore $S_b(x_r)$ is constant on every interval between consecutive primes, and

$$T_r=1+\sum_{i=1}^{t}(p_{i+1}-p_i)\prod_{j=1}^{i} p_j^{\alpha_j(r)}.$$

The coefficient $p_{i+1}-p_i$ is the number of integers $b$ for which the smooth divisor has exactly that prime-prefix product. This is why prime gaps become the natural weights in the implementation.

Step 3: Update the prime exponents from one binomial coefficient to the next

Consecutive coefficients satisfy the standard ratio

$$\binom{n}{r+1}=\binom{n}{r}\cdot \frac{n-r}{r+1}.$$

So the exponent vector $(\alpha_i(r))$ does not need to be recomputed from scratch. It changes only by the prime powers appearing in $n-r$ and $r+1$.

For each integer $m\le n$, the implementation precomputes a decomposition

$$m=p^{e}u,\qquad p \nmid u,$$

where $p$ is the smallest prime dividing $m$. One lookup then reveals four useful pieces of information: which prime is involved, the prime-power block $p^e$, its modular inverse, and the remaining cofactor $u$. Repeating this stripping process factors any number into distinct prime-power blocks in a small number of table lookups.

Because $M$ is prime and all relevant primes satisfy $p\le n<M$, each block $p^e$ is invertible modulo $M$, so division by the denominator is handled as multiplication by modular inverses.

Step 4: Maintain the weighted prefix products with a segment tree

For each prime $p_i$, define the current local factor

$$A_i(r)=p_i^{\alpha_i(r)},\qquad w_i=p_{i+1}-p_i.$$

A leaf stores the pair

$$A_i(r),\qquad w_iA_i(r).$$

For an interval $[L,R]$, store

$$P_{L,R}=\prod_{k=L}^{R} A_k(r),$$

$$Q_{L,R}=\sum_{i=L}^{R} w_i\prod_{k=L}^{i} A_k(r).$$

If the interval is split into left and right halves, the merge rule is

$$P=P_{\text{left}}P_{\text{right}},\qquad Q=Q_{\text{left}}+P_{\text{left}}Q_{\text{right}}.$$

At the root this gives

$$Q_{1,t}=\sum_{i=1}^{t} w_i\prod_{j=1}^{i} A_j(r),$$

hence

$$T_r=1+Q_{1,t}\pmod{M}.$$

Only the leaves belonging to primes dividing $n-r$ or $r+1$ change from one row to the next, so each transition touches only a small part of the tree.

Step 5: Use symmetry of the binomial row

Since

$$\binom{n}{r}=\binom{n}{n-r},$$

we also have $T_r=T_{n-r}$. Therefore it is enough to sum from $r=0$ to $\lfloor n/2\rfloor$, doubling every non-central term. When $n$ is even, the middle coefficient is added only once.

Worked Example: $n=11$ and $r=4$

Here

$$x_4=\binom{11}{4}=330=2\cdot 3\cdot 5\cdot 11.$$

The primes up to $11$ are $2,3,5,7,11$, and the sentinel is $12$. Because the exponent of $7$ is zero, the prime-prefix products are

$$2,\quad 2\cdot 3=6,\quad 2\cdot 3\cdot 5=30,\quad 30,\quad 30\cdot 11=330.$$

Thus

$$T_4=1+(3-2)\cdot 2+(5-3)\cdot 6+(7-5)\cdot 30+(11-7)\cdot 30+(12-11)\cdot 330=525.$$

Directly listing the smooth divisors for $b=1,\dots,11$ gives

$$1,\,2,\,6,\,6,\,30,\,30,\,30,\,30,\,30,\,30,\,330,$$

whose sum is again $525$.

The next coefficient is

$$\binom{11}{5}=\binom{11}{4}\cdot \frac{7}{5}=462,$$

so only the prime-power contributions attached to $5$ and $7$ change. This illustrates why the row-to-row update is cheap. Summing the whole row symmetrically yields the checkpoint

$$F(11)=3132.$$

How the Code Works

The C++, Python, and Java implementations follow the same numerical strategy. They first generate all primes up to $n$, record the gap from each prime to the next one, and precompute for every integer up to $n$ how to peel off one smallest-prime-power block together with the modular factors needed to apply or remove that block.

The main loop starts from $\binom{n}{0}=1$, where every local factor is $1$. For the current $r$, the implementation reads the segment-tree root, adds $1$, and obtains $T_r$. It then adds that value to the running total, using binomial symmetry to double non-central terms.

To advance from $r$ to $r+1$, the implementation factors $n-r$ and $r+1$ through the precomputed stripping tables, combines all multiplicative changes prime by prime, and updates only the affected leaves. The C++ and Java versions contain the full sieve-and-tree algorithm directly, while the Python implementation delegates to the same compiled computation path.

Complexity Analysis

Let $\pi(n)$ be the number of primes up to $n$. The preprocessing tables for the sieve and the stripped factorizations require $O(n)$ memory, and the segment tree needs $O(\pi(n))$ more.

Across all integers up to $n$, the total amount of prime-stripping work is $O(n\log\log n)$. Each distinct prime touched during a row transition causes one segment-tree update, and each such update costs $O(\log \pi(n))$. So the overall running time is about

$$O\!\bigl(n\log\log n\log \pi(n)\bigr),$$

which in practice behaves like a near-$O(n\log n)$ method. The memory usage is

$$O(n+\pi(n)).$$

Footnotes and References

Problem page: https://projecteuler.net/problem=468
Smooth number: Wikipedia — Smooth number
Binomial coefficient: Wikipedia — Binomial coefficient
Prime factorization: Wikipedia — Prime factorization
Segment tree: Wikipedia — Segment tree

Problem 468 source code

C++

#include <algorithm>
#include <chrono>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <vector>

namespace {

using u32 = std::uint32_t;
using u64 = std::uint64_t;

constexpr u32 kMod = 1'000'000'993U;
constexpr u32 kDefaultN = 11'111'111U;

struct Options {
    u32 n = kDefaultN;
    bool run_checkpoints = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

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

inline u32 add_mod(u32 a, u32 b) {
    const u32 sum = a + b;
    return (sum >= kMod) ? (sum - kMod) : sum;
}

u32 mod_pow(u32 base, u32 exponent) {
    u32 result = 1U;
    u32 b = base;
    u32 e = exponent;
    while (e > 0U) {
        if ((e & 1U) != 0U) {
            result = mul_mod(result, b);
        }
        b = mul_mod(b, b);
        e >>= 1U;
    }
    return result;
}

bool parse_u32_after_prefix(const std::string& arg, const char* prefix, u32& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const u32 digit = static_cast<u32>(c - '0');
        parsed = parsed * 10ULL + static_cast<u64>(digit);
        if (parsed > static_cast<u64>(std::numeric_limits<u32>::max())) {
            return false;
        }
    }

    value = static_cast<u32>(parsed);
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u32 parsed = 0;
    if (!parse_u32_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(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 (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }

        u32 parsed_u32 = 0;
        if (parse_u32_after_prefix(arg, "--n=", parsed_u32)) {
            options.n = parsed_u32;
            continue;
        }

        unsigned parsed_unsigned = 0U;
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
            options.requested_threads = parsed_unsigned;
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    if (options.n == 0U) {
        std::cerr << "--n must be >= 1.\n";
        return false;
    }

    return true;
}

u64 largest_b_smooth_divisor(u32 b, u64 x) {
    u64 result = 1ULL;
    u64 value = x;

    for (u64 p = 2ULL; p * p <= value; ++p) {
        if (value % p != 0ULL) {
            continue;
        }

        u64 prime_power = 1ULL;
        while (value % p == 0ULL) {
            value /= p;
            prime_power *= p;
        }

        if (p <= static_cast<u64>(b)) {
            result *= prime_power;
        }
    }

    if (value > 1ULL && value <= static_cast<u64>(b)) {
        result *= value;
    }

    return result;
}

struct SieveData {
    u32 n = 0U;
    std::vector<u32> weights;
    std::vector<u32> strip_next;
    std::vector<u32> strip_pos;
    std::vector<u32> strip_mul;
    std::vector<u32> strip_inv;
};

SieveData build_sieve_data(u32 n) {
    SieveData data;
    data.n = n;
    std::vector<u32> spf(static_cast<std::size_t>(n) + 1ULL, 0U);
    std::vector<int> prime_pos(static_cast<std::size_t>(n) + 1ULL, -1);
    std::vector<u32> primes;

    primes.reserve(static_cast<std::size_t>(n / 10U));

    for (u32 i = 2U; i <= n; ++i) {
        if (spf[i] == 0U) {
            spf[i] = i;
            primes.push_back(i);
        }

        for (const u32 p : primes) {
            const u64 composite = static_cast<u64>(i) * static_cast<u64>(p);
            if (composite > static_cast<u64>(n)) {
                break;
            }
            spf[static_cast<std::size_t>(composite)] = p;
            if (p == spf[i]) {
                break;
            }
        }
    }

    std::vector<u32> inv_prime_by_pos(primes.size());
    data.weights.resize(primes.size());
    for (std::size_t i = 0; i < primes.size(); ++i) {
        const u32 p = primes[i];
        prime_pos[p] = static_cast<int>(i);
        inv_prime_by_pos[i] = mod_pow(p, kMod - 2U);
        const u32 next_prime = (i + 1U < primes.size()) ? primes[i + 1U] : (n + 1U);
        data.weights[i] = next_prime - p;
    }

    data.strip_next.assign(static_cast<std::size_t>(n) + 1ULL, 1U);
    data.strip_pos.assign(static_cast<std::size_t>(n) + 1ULL, 0U);
    data.strip_mul.assign(static_cast<std::size_t>(n) + 1ULL, 1U);
    data.strip_inv.assign(static_cast<std::size_t>(n) + 1ULL, 1U);

    for (u32 x = 2U; x <= n; ++x) {
        const u32 p = spf[x];
        const u32 pos = static_cast<u32>(prime_pos[p]);
        const u32 inv_p = inv_prime_by_pos[static_cast<std::size_t>(pos)];

        u32 value = x;
        u32 mul = 1U;
        u32 inv_mul = 1U;
        do {
            value /= p;
            mul = mul_mod(mul, p);
            inv_mul = mul_mod(inv_mul, inv_p);
        } while (value % p == 0U);

        data.strip_next[x] = value;
        data.strip_pos[x] = pos;
        data.strip_mul[x] = mul;
        data.strip_inv[x] = inv_mul;
    }

    return data;
}

class SegmentTree {
  public:
    explicit SegmentTree(const std::vector<u32>& weights) {
        leaf_count_ = static_cast<int>(weights.size());
        size_ = 1;
        while (size_ < leaf_count_) {
            size_ <<= 1;
        }

        nodes_.assign(static_cast<std::size_t>(2 * size_), Node{1U, 0U});

        for (int i = 0; i < leaf_count_; ++i) {
            nodes_[static_cast<std::size_t>(size_ + i)].sum =
                weights[static_cast<std::size_t>(i)] % kMod;
        }

        for (int node = size_ - 1; node >= 1; --node) {
            pull(node);
        }
    }

    void multiply_leaf(int leaf_index, u32 factor) {
        int node = size_ + leaf_index;
        Node& leaf = nodes_[static_cast<std::size_t>(node)];
        leaf.prod = mul_mod(leaf.prod, factor);
        leaf.sum = mul_mod(leaf.sum, factor);

        for (node >>= 1; node > 0; node >>= 1) {
            pull(node);
        }
    }

    u32 root_sum() const {
        return nodes_[1].sum;
    }

  private:
    struct Node {
        u32 prod;
        u32 sum;
    };

    void pull(int node) {
        const int left = node << 1;
        const int right = left | 1;
        const Node& lhs = nodes_[static_cast<std::size_t>(left)];
        const Node& rhs = nodes_[static_cast<std::size_t>(right)];
        Node& out = nodes_[static_cast<std::size_t>(node)];
        out.prod = mul_mod(lhs.prod, rhs.prod);
        out.sum = add_mod(lhs.sum, mul_mod(lhs.prod, rhs.sum));
    }

    int leaf_count_ = 0;
    int size_ = 0;
    std::vector<Node> nodes_;
};

class SolverState {
  public:
    explicit SolverState(const SieveData& data) : data_(data), tree_(data.weights) {}

    u32 current_row_contribution() const {
        const u32 v = tree_.root_sum();
        return (v + 1U == kMod) ? 0U : (v + 1U);
    }

    void apply_ratio(u32 numerator, u32 denominator) {
        u32 positions[16];
        u32 factors[16];
        int used = 0;

        auto append = [&](u32 value, bool divide) {
            while (value > 1U) {
                const u32 pos = data_.strip_pos[value];
                const u32 factor = divide ? data_.strip_inv[value] : data_.strip_mul[value];

                int slot = -1;
                for (int i = 0; i < used; ++i) {
                    if (positions[i] == pos) {
                        slot = i;
                        break;
                    }
                }

                if (slot < 0) {
                    positions[used] = pos;
                    factors[used] = factor;
                    ++used;
                } else {
                    factors[slot] = mul_mod(factors[slot], factor);
                }
                value = data_.strip_next[value];
            }
        };

        append(numerator, false);
        append(denominator, true);

        for (int i = 0; i < used; ++i) {
            if (factors[i] != 1U) {
                tree_.multiply_leaf(static_cast<int>(positions[i]), factors[i]);
            }
        }
    }

  private:
    const SieveData& data_;
    SegmentTree tree_;
};

u32 sum_symmetric(const SieveData& data, u32 n) {
    SolverState state(data);
    u32 total = 0U;
    const u32 half = n / 2U;

    for (u32 r = 0U;; ++r) {
        const u32 contribution = state.current_row_contribution();
        if (r == n - r) {
            total = add_mod(total, contribution);
        } else {
            total = add_mod(total, add_mod(contribution, contribution));
        }

        if (r == half) {
            break;
        }
        state.apply_ratio(n - r, r + 1U);
    }

    return total;
}

u32 solve_mod(const SieveData& data,
              u32 n,
              bool allow_multithreading,
              unsigned requested_threads) {
    (void)allow_multithreading;
    (void)requested_threads;
    return sum_symmetric(data, n);
}

u32 solve_for_n(u32 n, bool allow_multithreading, unsigned requested_threads) {
    SieveData data = build_sieve_data(n);
    return solve_mod(data, n, allow_multithreading, requested_threads);
}

bool expect_equal_u64(const char* label, u64 actual, u64 expected) {
    if (actual == expected) {
        return true;
    }

    std::cerr << "Checkpoint failed for " << label << ": got " << actual << ", expected "
              << expected << '\n';
    return false;
}

bool expect_equal_u32(const char* label, u32 actual, u32 expected) {
    return expect_equal_u64(label, static_cast<u64>(actual), static_cast<u64>(expected));
}

bool run_checkpoints() {
    bool ok = true;

    ok = ok && expect_equal_u64("S_1(10)", largest_b_smooth_divisor(1U, 10ULL), 1ULL);
    ok = ok && expect_equal_u64("S_4(2100)", largest_b_smooth_divisor(4U, 2100ULL), 12ULL);
    ok = ok && expect_equal_u64("S_17(2496144)", largest_b_smooth_divisor(17U, 2'496'144ULL),
                                5'712ULL);

    ok = ok && expect_equal_u32("F(11)", solve_for_n(11U, false, 1U), 3'132U);
    ok = ok &&
         expect_equal_u32("F(1111) mod M", solve_for_n(1'111U, false, 1U), 706'036'312U);
    ok = ok &&
         expect_equal_u32("F(111111) mod M", solve_for_n(111'111U, false, 1U), 22'156'169U);

    if (!ok) {
        std::cerr << "At least one checkpoint failed.\n";
    }
    return ok;
}

}  // namespace

int main(int argc, char** argv) {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints()) {
        return 1;
    }

    const auto start = std::chrono::steady_clock::now();
    const SieveData data = build_sieve_data(options.n);
    const u32 answer =
        solve_mod(data, options.n, options.allow_multithreading, options.requested_threads);
    const auto end = std::chrono::steady_clock::now();

    const std::chrono::duration<long double> elapsed = end - start;
    std::cout << "F(" << options.n << ") mod " << kMod << " = " << answer << '\n';
    std::cout << "Elapsed: " << elapsed.count() << " seconds\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;

public class Euler468 {

    static final long MOD = 1000000993L;

    static long mulMod(long a, long b) {
        return (a * b) % MOD;
    }

    static long modPow(long base, long exponent) {
        long result = 1;
        long b = base % MOD;
        long e = exponent;
        while (e > 0) {
            if ((e & 1) != 0) {
                result = mulMod(result, b);
            }
            b = mulMod(b, b);
            e >>= 1;
        }
        return result;
    }

    static class SegmentTree {
        int leafCount;
        int size;
        long[] prod;
        long[] sum;

        SegmentTree(int[] weights) {
            leafCount = weights.length;
            size = 1;
            while (size < leafCount) {
                size <<= 1;
            }
            prod = new long[2 * size];
            sum = new long[2 * size];
            for (int i = 0; i < 2 * size; i++)
                prod[i] = 1;
            for (int i = 0; i < 2 * size; i++)
                sum[i] = 0;

            for (int i = 0; i < leafCount; i++) {
                sum[size + i] = weights[i] % MOD;
            }

            for (int node = size - 1; node >= 1; node--) {
                pull(node);
            }
        }

        void pull(int node) {
            int left = node << 1;
            int right = left | 1;
            prod[node] = mulMod(prod[left], prod[right]);
            sum[node] = (sum[left] + mulMod(prod[left], sum[right])) % MOD;
        }

        void multiplyLeaf(int leafIndex, long factor) {
            int node = size + leafIndex;
            prod[node] = mulMod(prod[node], factor);
            sum[node] = mulMod(sum[node], factor);

            for (node >>= 1; node > 0; node >>= 1) {
                pull(node);
            }
        }

        long rootSum() {
            return sum[1];
        }
    }

    static class SieveData {
        int n;
        int[] weights;
        int[] stripNext;
        int[] stripPos;
        long[] stripMul;
        long[] stripInv;
    }

    static SieveData buildSieveData(int n) {
        SieveData data = new SieveData();
        data.n = n;
        int[] spf = new int[n + 1];
        int[] primePos = new int[n + 1];
        for (int i = 0; i <= n; i++)
            primePos[i] = -1;

        List<Integer> primes = new ArrayList<>(n / 10);
        for (int i = 2; i <= n; i++) {
            if (spf[i] == 0) {
                spf[i] = i;
                primes.add(i);
            }
            for (int p : primes) {
                long composite = (long) i * p;
                if (composite > n)
                    break;
                spf[(int) composite] = p;
                if (p == spf[i])
                    break;
            }
        }

        int pCount = primes.size();
        long[] invPrimeByPos = new long[pCount];
        data.weights = new int[pCount];
        for (int i = 0; i < pCount; i++) {
            int p = primes.get(i);
            primePos[p] = i;
            invPrimeByPos[i] = modPow(p, MOD - 2);
            int nextPrime = (i + 1 < pCount) ? primes.get(i + 1) : (n + 1);
            data.weights[i] = nextPrime - p;
        }

        data.stripNext = new int[n + 1];
        data.stripPos = new int[n + 1];
        data.stripMul = new long[n + 1];
        data.stripInv = new long[n + 1];

        for (int x = 2; x <= n; x++) {
            int p = spf[x];
            int pos = primePos[p];
            long invP = invPrimeByPos[pos];

            int value = x;
            long mul = 1;
            long invMul = 1;
            do {
                value /= p;
                mul = mulMod(mul, p);
                invMul = mulMod(invMul, invP);
            } while (value % p == 0);

            data.stripNext[x] = value;
            data.stripPos[x] = pos;
            data.stripMul[x] = mul;
            data.stripInv[x] = invMul;
        }

        return data;
    }

    public static String solve() {
        int n = 11111111;
        SieveData data = buildSieveData(n);
        SegmentTree tree = new SegmentTree(data.weights);

        long total = 0;
        int half = n / 2;
        int r = 0;

        int[] positions = new int[16];
        long[] factors = new long[16];

        while (true) {
            long v = tree.rootSum();
            long contribution = (v + 1 == MOD) ? 0 : (v + 1);

            if (r == n - r) {
                total = (total + contribution) % MOD;
            } else {
                total = (total + contribution * 2) % MOD;
            }

            if (r == half)
                break;

            int used = 0;
            int val = n - r;
            while (val > 1) {
                int pos = data.stripPos[val];
                long factor = data.stripMul[val];
                int slot = -1;
                for (int i = 0; i < used; i++) {
                    if (positions[i] == pos) {
                        slot = i;
                        break;
                    }
                }
                if (slot < 0) {
                    positions[used] = pos;
                    factors[used] = factor;
                    used++;
                } else {
                    factors[slot] = mulMod(factors[slot], factor);
                }
                val = data.stripNext[val];
            }

            val = r + 1;
            while (val > 1) {
                int pos = data.stripPos[val];
                long factor = data.stripInv[val];
                int slot = -1;
                for (int i = 0; i < used; i++) {
                    if (positions[i] == pos) {
                        slot = i;
                        break;
                    }
                }
                if (slot < 0) {
                    positions[used] = pos;
                    factors[used] = factor;
                    used++;
                } else {
                    factors[slot] = mulMod(factors[slot], factor);
                }
                val = data.stripNext[val];
            }

            for (int i = 0; i < used; i++) {
                if (factors[i] != 1) {
                    tree.multiplyLeaf(positions[i], factors[i]);
                }
            }
            r++;
        }

        return Long.toString(total);
    }

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