Problem 565: Divisibility of Sum of Divisors

Project Euler Problem 565 Solution

EulerSolve provides an optimized solution for Project Euler Problem 565, Divisibility of Sum of Divisors, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For $N=10^{11}$ and the prime $p=2017$, we must compute $$S(N,p)=\sum_{\substack{1\le n\le N\\ p\mid \sigma(n)}} n,$$ where $\sigma(n)$ is the sum of positive divisors of $n$. A brute-force scan would require evaluating $\sigma(n)$ for every $n\le N$, which is far too expensive. The solution instead identifies the exact prime-power exponents that force divisibility by $2017$, rewrites the task as a weighted union of exact-valuation events, and then sums that union with inclusion-exclusion. Mathematical Approach Write the prime factorization of $n$ as $$n=\prod_q q^{a_q}.$$ Because the divisor-sum function is multiplicative, $$\sigma(n)=\prod_q \sigma\!\left(q^{a_q}\right).$$ Therefore $2017\mid \sigma(n)$ if and only if at least one prime $q$ satisfies $2017\mid \sigma(q^{a_q})$. Step 1: Reduce the problem to trigger prime powers Define the set of trigger pairs $$B=\left\{(q,e): e\ge 1,\ q^e\le N,\ 2017\mid \sigma(q^e)\right\},$$ where $q$ is prime. Then the target sum is the weighted union $$S(N,2017)=\sum_{\substack{1\le n\le N\\ \exists (q,e)\in B:\ v_q(n)=e}} n,$$ where $v_q(n)$ denotes the exponent of $q$ in $n$. Exact valuations matter: the condition depends on the precise exponent $a_q$, not merely on divisibility by $q^e$....

Detailed mathematical approach

Problem Summary

For $N=10^{11}$ and the prime $p=2017$, we must compute

$$S(N,p)=\sum_{\substack{1\le n\le N\\ p\mid \sigma(n)}} n,$$

where $\sigma(n)$ is the sum of positive divisors of $n$. A brute-force scan would require evaluating $\sigma(n)$ for every $n\le N$, which is far too expensive. The solution instead identifies the exact prime-power exponents that force divisibility by $2017$, rewrites the task as a weighted union of exact-valuation events, and then sums that union with inclusion-exclusion.

Mathematical Approach

Write the prime factorization of $n$ as

$$n=\prod_q q^{a_q}.$$

Because the divisor-sum function is multiplicative,

$$\sigma(n)=\prod_q \sigma\!\left(q^{a_q}\right).$$

Therefore $2017\mid \sigma(n)$ if and only if at least one prime $q$ satisfies $2017\mid \sigma(q^{a_q})$.

Step 1: Reduce the problem to trigger prime powers

Define the set of trigger pairs

$$B=\left\{(q,e): e\ge 1,\ q^e\le N,\ 2017\mid \sigma(q^e)\right\},$$

where $q$ is prime.

Then the target sum is the weighted union

$$S(N,2017)=\sum_{\substack{1\le n\le N\\ \exists (q,e)\in B:\ v_q(n)=e}} n,$$

where $v_q(n)$ denotes the exponent of $q$ in $n$. Exact valuations matter: the condition depends on the precise exponent $a_q$, not merely on divisibility by $q^e$.

The prime $q=2017$ never contributes, because

$$\sigma(2017^e)=1+2017+\cdots+2017^e\equiv 1 \pmod{2017}.$$

Step 2: Characterize the bad exponents with multiplicative order

For a prime $q\neq 2017$,

$$\sigma(q^e)=1+q+\cdots+q^e=\frac{q^{e+1}-1}{q-1}.$$

If $q\not\equiv 1 \pmod{2017}$, the denominator is invertible modulo $2017$, so

$$2017\mid \sigma(q^e)\iff q^{e+1}\equiv 1 \pmod{2017}.$$

Let $\operatorname{ord}_{2017}(q)$ be the multiplicative order of $q$ modulo $2017$. Then

$$2017\mid \sigma(q^e)\iff \operatorname{ord}_{2017}(q)\mid (e+1),$$

so the bad exponents are exactly

$$e=m\operatorname{ord}_{2017}(q)-1,\qquad m\ge 1.$$

Since $2017-1=2016=2^5\cdot 3^2\cdot 7$, every such order divides $2016$.

The residue class $q\equiv 1 \pmod{2017}$ behaves differently, because then

$$\sigma(q^e)\equiv e+1 \pmod{2017}.$$

That would require $e+1$ to be a multiple of $2017$, impossible under $q^e\le 10^{11}$. So primes congruent to $1$ modulo $2017$ contribute nothing.

Step 3: Separate the order-2 branch from the rest

If $q\equiv -1 \pmod{2017}$, then $\operatorname{ord}_{2017}(q)=2$, hence all odd exponents are formally bad.

For the target scale, however, only $e=1$ can occur. The smallest prime with $q\equiv -1 \pmod{2017}$ is $12101$, and

$$12101^3>10^{11}.$$

So the order-2 branch contributes only prime powers of the form $q^1$.

All remaining bad cases come from primes $q\not\equiv \pm 1 \pmod{2017}$ with order at least $3$, so their smallest bad exponent is at least $2$. Consequently such primes matter only when $q\le \sqrt{N}$, which is why this branch is searched only up to $\lfloor \sqrt{N}\rfloor$.

Step 4: Turn each trigger pair into a closed-form weighted sum

For each $(q,e)\in B$, define the event

$$A_{q,e}=\{n\le N:\ v_q(n)=e\}.$$

Every integer in this event has the form $n=q^e m$ with $q\nmid m$ and

$$m\le M=\left\lfloor \frac{N}{q^e}\right\rfloor.$$

If $T(x)=x(x+1)/2$ denotes the triangular-number sum, then

$$\sum_{n\in A_{q,e}} n =q^e\sum_{\substack{m\le M\\ q\nmid m}} m =q^e\left(T(M)-q\,T\!\left(\left\lfloor \frac{M}{q}\right\rfloor\right)\right).$$

This gives every single-event contribution in constant time once $q^e$ is known.

Step 5: Use pairwise inclusion-exclusion, and stop there exactly

Events with the same prime $q$ but different exponents are disjoint, so overlaps only arise from distinct primes. If $(q,e)$ and $(r,f)$ are trigger pairs with $q\neq r$, then

$$n=q^e r^f m,\qquad q\nmid m,\qquad r\nmid m,$$

with

$$m\le X=\left\lfloor \frac{N}{q^e r^f}\right\rfloor.$$

Therefore

$$\sum_{n\in A_{q,e}\cap A_{r,f}} n =q^e r^f\left(T(X)-q\,T\!\left(\left\lfloor \frac{X}{q}\right\rfloor\right)-r\,T\!\left(\left\lfloor \frac{X}{r}\right\rfloor\right)+qr\,T\!\left(\left\lfloor \frac{X}{qr}\right\rfloor\right)\right).$$

Hence

$$S(N,2017)=\sum_{(q,e)\in B}\sum_{n\in A_{q,e}} n-\sum_{\substack{(q,e),(r,f)\in B\\ q<r}}\sum_{n\in A_{q,e}\cap A_{r,f}} n,$$

because triple intersections cannot occur for $N=10^{11}$. The three smallest admissible trigger prime powers are

$$2311^2=5{,}340{,}721,\qquad 229^3=12{,}008{,}989,\qquad 3739^2=13{,}980{,}121,$$

and their product already exceeds $10^{11}$ by many orders of magnitude. So second-order inclusion-exclusion is exact, not approximate.

Worked Example

Two concrete trigger prime powers show the mechanism.

First, $12101\equiv -1 \pmod{2017}$, so order $2$ gives the bad exponent $e=1$, and indeed

$$\sigma(12101)=1+12101=12102=6\cdot 2017.$$

Second, $2311$ has multiplicative order $3$ modulo $2017$, so $e=2$ is bad, and

$$\sigma(2311^2)=1+2311+2311^2=5{,}343{,}033=2649\cdot 2017.$$

For $N=10^{11}$, their overlap uses

$$P=12101\cdot 2311^2=64{,}628{,}064{,}821\le 10^{11},$$

so exactly one multiplier $m=1$ is possible. The pair contribution is therefore precisely $P$, and it must be subtracted once from the two single-event contributions.

How the Code Works

The implementation first generates all primes up to $\lfloor \sqrt{N}\rfloor$. That prime list serves two purposes: it determines multiplicative orders modulo $2017$ for the non-$\pm 1$ branch, and it also drives a sieve on the arithmetic progression $2017k-1$ to collect every prime congruent to $-1$ modulo $2017$ up to $N$.

Next, the C++, Python, and Java implementations enumerate all trigger pairs $(q,e)$. For the order-2 branch they record only exponent $1$. For the other branch they compute the multiplicative order as a divisor of $2016$, list every exponent $e=m\operatorname{ord}_{2017}(q)-1$ with $q^e\le N$, and discard primes that cannot contribute.

Once the trigger list is known, the implementation accumulates all single-event formulas above. It then subtracts all admissible pair overlaps in three groups: order-2 with order-2, order-2 with other triggers, and other triggers with each other. The Python implementation delegates execution to the compiled solver, while the Java implementation mirrors the same arithmetic directly with arbitrary-precision integers.

The solver also verifies the published checkpoints

$$S(10^6,2017)=150850429,\qquad S(10^9,2017)=249652238344557,$$

before evaluating the final $N=10^{11}$ case.

Complexity Analysis

Generating all primes up to $\sqrt{N}$ costs $O(\sqrt{N}\log\log N)$ time and $O(\sqrt{N})$ memory. Sieving the progression $2017k-1$ up to $N$ uses an array of length $\lfloor (N+1)/2017\rfloor$, so that stage needs $O(N/2017)$ memory and about $O((N/2017)\log\log N)$ marking work. The remaining order checks and inclusion-exclusion sums run only over the finite trigger list and its admissible pairs. For the target parameters, those later steps are smaller than the two sieve phases.

Footnotes and References

Problem page: https://projecteuler.net/problem=565
Divisor function: Wikipedia — Divisor function
Multiplicative order: Wikipedia — Multiplicative order
$p$-adic valuation: Wikipedia — $p$-adic valuation
Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle

Problem 565 source code

C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>
#include <tuple>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;
using i64 = std::int64_t;

static std::string to_string_u128(u128 value) {
    if (value == 0) return "0";
    std::string s;
    while (value > 0) {
        const int digit = static_cast<int>(value % 10);
        s.push_back(static_cast<char>('0' + digit));
        value /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

static u128 tri(u64 n) {
    // n*(n+1)/2 in 128-bit.
    return static_cast<u128>(n) * static_cast<u128>(n + 1) / 2;
}

static u64 mod_pow(u64 a, u64 e, u64 mod) {
    u64 r = 1 % mod;
    a %= mod;
    while (e) {
        if (e & 1) r = (u128)r * a % mod;
        a = (u128)a * a % mod;
        e >>= 1;
    }
    return r;
}

static i64 egcd(i64 a, i64 b, i64& x, i64& y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    i64 x1 = 0, y1 = 0;
    const i64 g = egcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

static u64 mod_inv(u64 a, u64 mod) {
    i64 x = 0, y = 0;
    const i64 g = egcd(static_cast<i64>(a), static_cast<i64>(mod), x, y);
    assert(g == 1);
    i64 res = x % static_cast<i64>(mod);
    if (res < 0) res += static_cast<i64>(mod);
    return static_cast<u64>(res);
}

static std::vector<int> sieve_primes(int n) {
    std::vector<bool> is_prime(n + 1, true);
    is_prime[0] = is_prime[1] = false;
    for (int i = 2; 1LL * i * i <= n; ++i) {
        if (!is_prime[i]) continue;
        for (int j = i * i; j <= n; j += i) is_prime[j] = false;
    }
    std::vector<int> primes;
    for (int i = 2; i <= n; ++i) {
        if (is_prime[i]) primes.push_back(i);
    }
    return primes;
}

static u64 order_mod_prime(u64 a, u64 p) {
    // p is prime. Return multiplicative order of a mod p (a != 0 mod p).
    assert(p == 2017);
    const u64 mod = p;
    u64 ord = p - 1;  // 2016 = 2^5 * 3^2 * 7
    const u64 factors[] = {2, 3, 7};
    for (u64 f : factors) {
        while (ord % f == 0 && mod_pow(a, ord / f, mod) == 1) ord /= f;
    }
    return ord;
}

// Sieve primes q <= N with q ≡ -1 (mod p), using a sieve over k where q = p*k - 1.
static std::vector<u64> primes_mod_minus_one(u64 N, u64 p, const std::vector<int>& primes_up_to_sqrtN) {
    assert(p == 2017);
    const u64 k_max = (N + 1) / p;  // p*k - 1 <= N  <=>  k <= (N+1)/p
    std::vector<std::uint8_t> is_comp(k_max + 1, 0);

    for (int r : primes_up_to_sqrtN) {
        if (r == static_cast<int>(p)) continue;
        const u64 rr = static_cast<u64>(r) * static_cast<u64>(r);
        if (rr > N) break;

        const u64 inv = mod_inv(p % static_cast<u64>(r), static_cast<u64>(r));  // k ≡ inv (mod r)
        const u64 k_min = (rr + 1 + p - 1) / p;                                 // ceil((r^2+1)/p)

        u64 k = inv;
        if (k < k_min) {
            const u64 step = static_cast<u64>(r);
            const u64 t = (k_min - k + step - 1) / step;
            k += t * step;
        }
        for (; k <= k_max; k += static_cast<u64>(r)) {
            is_comp[k] = 1;
        }
    }

    std::vector<u64> out;
    out.reserve(2'200'000);
    for (u64 k = 1; k <= k_max; ++k) {
        if (is_comp[k]) continue;
        const u64 q = p * k - 1;
        if (q >= 2 && q <= N) out.push_back(q);
    }
    return out;
}

struct BadPrime {
    u64 q = 0;
    std::vector<int> bad_exps;  // exact exponents e with p|sigma(q^e) within q^e<=N
};

static std::vector<BadPrime> bad_primes_other(u64 N, u64 p, int sqrtN, const std::vector<int>& primes_up_to_sqrtN) {
    std::vector<BadPrime> bad;
    for (int q : primes_up_to_sqrtN) {
        if (q > sqrtN) break;
        if (q == static_cast<int>(p)) continue;
        const u64 a = static_cast<u64>(q) % p;
        if (a == 1 || a == p - 1) continue;  // order 1 is irrelevant here; order 2 handled separately

        const u64 ord = order_mod_prime(a, p);

        // max exponent with q^e <= N
        int emax = 0;
        u64 pw = 1;
        while (pw <= N / static_cast<u64>(q)) {
            pw *= static_cast<u64>(q);
            ++emax;
        }

        std::vector<int> exps;
        for (int m = 1;; ++m) {
            const u64 e = ord * static_cast<u64>(m) - 1;
            if (e > static_cast<u64>(emax)) break;
            exps.push_back(static_cast<int>(e));
        }
        if (exps.empty()) continue;

        bad.push_back(BadPrime{static_cast<u64>(q), std::move(exps)});
    }
    return bad;
}

static u128 sum_exact_bad_exp(u64 N, u64 q, int e) {
    // Sum_{n<=N, v_q(n)=e} n = q^e * sum_{m<=N/q^e, q∤m} m.
    u128 qpow = 1;
    for (int i = 0; i < e; ++i) qpow *= q;
    const u64 M = static_cast<u64>(static_cast<u128>(N) / qpow);
    const u64 Mq = M / q;
    const u128 sum_not_div_q = tri(M) - static_cast<u128>(q) * tri(Mq);
    return qpow * sum_not_div_q;
}

static u128 sum_pair_exact(u64 N, u64 q, int eq, u64 r, int er) {
    // Sum_{n<=N, v_q(n)=eq and v_r(n)=er} n.
    u128 qpow = 1;
    for (int i = 0; i < eq; ++i) qpow *= q;
    u128 rpow = 1;
    for (int i = 0; i < er; ++i) rpow *= r;
    const u128 P = qpow * rpow;
    if (P > N) return 0;
    const u64 M = static_cast<u64>(static_cast<u128>(N) / P);

    const u64 Mq = M / q;
    const u64 Mr = M / r;
    const u64 Mqr = M / (q * r);

    const u128 sum_coprime = tri(M) - static_cast<u128>(q) * tri(Mq) - static_cast<u128>(r) * tri(Mr) +
                             static_cast<u128>(q) * static_cast<u128>(r) * tri(Mqr);
    return P * sum_coprime;
}

static u128 S(u64 N, u64 p) {
    const int sqrtN = static_cast<int>(std::sqrt(static_cast<long double>(N)));
    const auto primes = sieve_primes(sqrtN);

    // Order-2 primes: q ≡ -1 (mod p) with exact exponent 1 (since (min such q)^3 > 1e11).
    const auto p2 = primes_mod_minus_one(N, p, primes);

    // Other bad primes within sqrt(N): exact exponents e = ord(q)*m - 1 within range.
    auto other = bad_primes_other(N, p, sqrtN, primes);

    // Singles
    u128 singles = 0;
    for (u64 q : p2) {
        singles += sum_exact_bad_exp(N, q, 1);
    }
    for (const auto& bp : other) {
        for (int e : bp.bad_exps) singles += sum_exact_bad_exp(N, bp.q, e);
    }

    // Pairs: only need q<r once, but given N=1e11 there are no triple intersections.
    u128 pairs = 0;

    // order2-order2 pairs: q*r <= N, so q must be <= sqrt(N).
    const u64 sqrtN_u = static_cast<u64>(sqrtN);
    const auto it_small_end = std::upper_bound(p2.begin(), p2.end(), sqrtN_u);
    for (auto itq = p2.begin(); itq != it_small_end; ++itq) {
        const u64 q = *itq;
        const u64 lim = N / q;
        auto itr_end = std::upper_bound(itq + 1, p2.end(), lim);
        for (auto itr = itq + 1; itr != itr_end; ++itr) {
            const u64 r = *itr;
            pairs += sum_pair_exact(N, q, 1, r, 1);
        }
    }

    // order2-other pairs: very few due to size constraints (but do it generically).
    for (const auto& bp : other) {
        for (int e : bp.bad_exps) {
            u128 rpow = 1;
            for (int i = 0; i < e; ++i) rpow *= bp.q;
            if (rpow > N) continue;
            const u64 lim = static_cast<u64>(static_cast<u128>(N) / rpow);
            for (u64 q : p2) {
                if (q > lim) break;
                pairs += sum_pair_exact(N, q, 1, bp.q, e);
            }
        }
    }

    // other-other pairs: (empirically none for N=1e11, but keep safe).
    for (std::size_t i = 0; i < other.size(); ++i) {
        for (std::size_t j = i + 1; j < other.size(); ++j) {
            for (int ei : other[i].bad_exps) {
                for (int ej : other[j].bad_exps) {
                    pairs += sum_pair_exact(N, other[i].q, ei, other[j].q, ej);
                }
            }
        }
    }

    return singles - pairs;
}

}  // namespace

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    const u64 p = 2017;

    // Validation points from the statement.
    assert(to_string_u128(S(1'000'000ULL, p)) == "150850429");
    assert(to_string_u128(S(1'000'000'000ULL, p)) == "249652238344557");

    std::cout << to_string_u128(S(100'000'000'000ULL, p)) << '\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.math.BigInteger;
import java.util.ArrayList;
import java.util.List;

public class Euler565 {
    static BigInteger tri(long n) {
        return BigInteger.valueOf(n).multiply(BigInteger.valueOf(n + 1)).divide(BigInteger.valueOf(2));
    }

    static long modPow(long a, long e, long mod) {
        long r = 1 % mod;
        a %= mod;
        while (e > 0) {
            if ((e & 1) != 0)
                r = (long) ((BigInteger.valueOf(r).multiply(BigInteger.valueOf(a))).remainder(BigInteger.valueOf(mod))
                        .longValue());
            a = (long) ((BigInteger.valueOf(a).multiply(BigInteger.valueOf(a))).remainder(BigInteger.valueOf(mod))
                    .longValue());
            e >>= 1;
        }
        return r;
    }

    static long[] egcd(long a, long b) {
        if (b == 0)
            return new long[] { 1, 0, a };
        long[] res = egcd(b, a % b);
        long x1 = res[0];
        long y1 = res[1];
        long g = res[2];
        long x = y1;
        long y = x1 - (a / b) * y1;
        return new long[] { x, y, g };
    }

    static long modInv(long a, long mod) {
        long[] res = egcd(a, mod);
        long x = res[0];
        long r = x % mod;
        if (r < 0)
            r += mod;
        return r;
    }

    static List<Integer> sievePrimes(int n) {
        boolean[] isPrime = new boolean[n + 1];
        for (int i = 2; i <= n; i++)
            isPrime[i] = true;
        for (int i = 2; i * i <= n; i++) {
            if (isPrime[i]) {
                for (int j = i * i; j <= n; j += i) {
                    isPrime[j] = false;
                }
            }
        }
        List<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= n; i++) {
            if (isPrime[i])
                primes.add(i);
        }
        return primes;
    }

    static long orderModPrime(long a, long p) {
        long ord = p - 1;
        long[] factors = { 2, 3, 7 };
        for (long f : factors) {
            while (ord % f == 0 && modPow(a, ord / f, p) == 1) {
                ord /= f;
            }
        }
        return ord;
    }

    static List<Long> primesModMinusOne(long N, long p, List<Integer> primesUpToSqrtN) {
        long kMax = (N + 1) / p;
        byte[] isComp = new byte[(int) (kMax + 1)];

        for (int r : primesUpToSqrtN) {
            if (r == p)
                continue;
            long rr = (long) r * r;
            if (rr > N)
                break;

            long inv = modInv(p % r, r);
            long kMin = (rr + 1 + p - 1) / p;

            long k = inv;
            if (k < kMin) {
                long t = (kMin - k + r - 1) / r;
                k += t * r;
            }
            for (; k <= kMax; k += r) {
                isComp[(int) k] = 1;
            }
        }

        List<Long> out = new ArrayList<>(2200000);
        for (long k = 1; k <= kMax; k++) {
            if (isComp[(int) k] == 0) {
                long q = p * k - 1;
                if (q >= 2 && q <= N)
                    out.add(q);
            }
        }
        return out;
    }

    static class BadPrime {
        long q;
        List<Integer> badExps;

        BadPrime(long q, List<Integer> badExps) {
            this.q = q;
            this.badExps = badExps;
        }
    }

    static List<BadPrime> badPrimesOther(long N, long p, int sqrtN, List<Integer> primes) {
        List<BadPrime> bad = new ArrayList<>();
        for (int qInt : primes) {
            if (qInt > sqrtN)
                break;
            long q = qInt;
            if (q == p)
                continue;
            long a = q % p;
            if (a == 1 || a == p - 1)
                continue;

            long ord = orderModPrime(a, p);
            int emax = 0;
            long pw = 1;
            while (pw <= N / q) {
                pw *= q;
                emax++;
            }

            List<Integer> exps = new ArrayList<>();
            for (int m = 1;; m++) {
                long e = ord * m - 1;
                if (e > emax)
                    break;
                exps.add((int) e);
            }
            if (!exps.isEmpty()) {
                bad.add(new BadPrime(q, exps));
            }
        }
        return bad;
    }

    static BigInteger sumExactBadExp(long N, long q, int e) {
        BigInteger qpow = BigInteger.ONE;
        for (int i = 0; i < e; i++)
            qpow = qpow.multiply(BigInteger.valueOf(q));
        long M = BigInteger.valueOf(N).divide(qpow).longValue();
        long Mq = M / q;
        BigInteger sumNotDivQ = tri(M).subtract(BigInteger.valueOf(q).multiply(tri(Mq)));
        return qpow.multiply(sumNotDivQ);
    }

    static BigInteger sumPairExact(long N, long q, int eq, long r, int er) {
        BigInteger qpow = BigInteger.ONE;
        for (int i = 0; i < eq; i++)
            qpow = qpow.multiply(BigInteger.valueOf(q));
        BigInteger rpow = BigInteger.ONE;
        for (int i = 0; i < er; i++)
            rpow = rpow.multiply(BigInteger.valueOf(r));
        BigInteger P = qpow.multiply(rpow);

        if (P.compareTo(BigInteger.valueOf(N)) > 0)
            return BigInteger.ZERO;

        long M = BigInteger.valueOf(N).divide(P).longValue();
        long Mq = M / q;
        long Mr = M / r;
        long Mqr = M / (q * r);

        BigInteger sumCoprime = tri(M)
                .subtract(BigInteger.valueOf(q).multiply(tri(Mq)))
                .subtract(BigInteger.valueOf(r).multiply(tri(Mr)))
                .add(BigInteger.valueOf(q).multiply(BigInteger.valueOf(r)).multiply(tri(Mqr)));

        return P.multiply(sumCoprime);
    }

    static int upperBound(List<Long> list, int start, long val) {
        int left = start;
        int right = list.size();
        while (left < right) {
            int mid = left + (right - left) / 2;
            if (list.get(mid) > val)
                right = mid;
            else
                left = mid + 1;
        }
        return left;
    }

    static BigInteger computeS(long N, long p) {
        int sqrtN = (int) Math.sqrt(N);
        List<Integer> primes = sievePrimes(sqrtN);

        List<Long> p2 = primesModMinusOne(N, p, primes);
        List<BadPrime> other = badPrimesOther(N, p, sqrtN, primes);

        BigInteger singles = BigInteger.ZERO;
        for (long q : p2) {
            singles = singles.add(sumExactBadExp(N, q, 1));
        }
        for (BadPrime bp : other) {
            for (int e : bp.badExps) {
                singles = singles.add(sumExactBadExp(N, bp.q, e));
            }
        }

        BigInteger pairs = BigInteger.ZERO;

        int itSmallEnd = upperBound(p2, 0, sqrtN);
        for (int i = 0; i < itSmallEnd; i++) {
            long q = p2.get(i);
            long lim = N / q;
            int itrEnd = upperBound(p2, i + 1, lim);
            for (int j = i + 1; j < itrEnd; j++) {
                long r = p2.get(j);
                pairs = pairs.add(sumPairExact(N, q, 1, r, 1));
            }
        }

        for (BadPrime bp : other) {
            for (int e : bp.badExps) {
                BigInteger rpow = BigInteger.ONE;
                for (int i = 0; i < e; i++)
                    rpow = rpow.multiply(BigInteger.valueOf(bp.q));
                if (rpow.compareTo(BigInteger.valueOf(N)) > 0)
                    continue;
                long lim = BigInteger.valueOf(N).divide(rpow).longValue();
                for (long q : p2) {
                    if (q > lim)
                        break;
                    pairs = pairs.add(sumPairExact(N, q, 1, bp.q, e));
                }
            }
        }

        for (int i = 0; i < other.size(); i++) {
            for (int j = i + 1; j < other.size(); j++) {
                for (int ei : other.get(i).badExps) {
                    for (int ej : other.get(j).badExps) {
                        pairs = pairs.add(sumPairExact(N, other.get(i).q, ei, other.get(j).q, ej));
                    }
                }
            }
        }

        return singles.subtract(pairs);
    }

    public static String solve() {
        return computeS(100000000000L, 2017).toString();
    }

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