Day 178-Q4 — 2026-10-14 | f(p^e)=p xor e の総和を強力数法で（g=φ に 2 だけ補正を掛けて g(p)=f(p) を合わせ、h を強力数の上だけで回す。min_25 篩と同じ答え・同じ速さ）

📋 コピーimport sys

MOD = 10**9 + 7


def main():
    N = int(sys.stdin.readline())
    # ---- 1) φ の前計算（篩）: T ≈ 0.5·N^(2/3)（下の測定で決めた）
    T = min(N, max(100, int(0.5 * N ** (2 / 3))))
    phi = list(range(T + 1))
    comp = bytearray(T + 1)
    for p in range(2, T + 1):
        if not comp[p]:
            comp[p * p::p] = b"\x01" * len(range(p * p, T + 1, p)) if p * p <= T else b""
            for j in range(p, T + 1, p):
                phi[j] -= phi[j] // p
    pre = [0] * (T + 1)
    acc = 0
    for i in range(1, T + 1):
        acc += phi[i]
        pre[i] = acc % MOD
    # ---- 2) 杜教篩: m > T の S_φ(m) を m = N//k の k で memo
    big = {}
    inv2 = (MOD + 1) // 2

    def S_phi(m):
        if m <= T:
            return pre[m]
        k = N // m
        if k in big:
            return big[k]
        # 小さい m から順に計算すれば再帰は浅い（下で k の大きい順に呼ぶ）
        res = m % MOD * ((m + 1) % MOD) % MOD * inv2 % MOD
        d = 2
        while d <= m:
            q = m // d
            nd = m // q
            res -= (nd - d + 1) * S_phi(q)
            d = nd + 1
        res %= MOD
        big[k] = res
        return res

    for k in range(N // (T + 1), 0, -1):       # m = N//k が大きくなる順 → 再帰は 1 段で済む
        if N // k > T:
            S_phi(N // k)

    def S_g(m):                                 # g = φ * u、u(2) = 2 だけ非零 → S_g(m) = S_φ(m) + 2 S_φ(m//2)
        return (S_phi(m) + 2 * S_phi(m // 2)) % MOD

    # ---- 3) 素数 p <= √N ごとに h(p^e) を作る
    r = int(N ** 0.5)
    while r * r > N:
        r -= 1
    while (r + 1) * (r + 1) <= N:
        r += 1
    sieve = bytearray([1]) * (r + 1)
    primes = []
    for i in range(2, r + 1):
        if sieve[i]:
            primes.append(i)
            sieve[i * i::i] = bytearray(len(range(i * i, r + 1, i)))
    hcoef = []
    for p in primes:
        pe = [1]
        while pe[-1] <= N // p:
            pe.append(pe[-1] * p)
        E = len(pe) - 1
        if p == 2:
            g = [1] + [3 if e == 1 else pe[e] for e in range(1, E + 1)]
        else:
            g = [1] + [pe[e] - pe[e - 1] for e in range(1, E + 1)]
        h = [1, 0]
        for e in range(2, E + 1):
            v = p ^ e
            for i in range(1, e + 1):
                v -= g[i] * h[e - i]
            h.append(v % MOD)
        hcoef.append(h)
    # ---- 4) 強力数を DFS で列挙し Σ h(n) S_g(N//n)
    ans = S_g(N)
    stack = [(1, 1, 0)]                         # (n, h(n), 次に使える素数の添字)
    while stack:
        n, hn, i = stack.pop()
        while i < len(primes):
            p = primes[i]
            if n * p * p > N:
                break
            h = hcoef[i]
            m = n * p * p
            e = 2
            while m <= N:
                hv = hn * h[e] % MOD
                if hv:
                    ans += hv * S_g(N // m)
                    stack.append((m, hv, i + 1))
                m *= p
                e += 1
            i += 1
    print(ans % MOD)


main()