Day 100-Q4 | 杜教篩（Du's Sieve）

問題

オイラーのトーシェント関数 $\varphi(i)$ について、

$$S(N) = \sum_{i=1}^{N} \varphi(i) \pmod{998244353}$$

を求めよ。$N$ が非常に大きいため素朴な篩では間に合わない。杜教篩（Du's Sieve）を用いて $O(N^{2/3})$ 時間で $S(N)$ を求めよ。

入力形式

制約

$1 \le N \le 10^{11}$

入出力例

入力例1

出力例1

$\varphi(1)+\cdots+\varphi(10)=1+1+2+2+4+2+6+4+6+4=32$。入力例2: $N=1$ → 出力 $1$。

概念図

ヒント（段階的開示）

ヒント1: 方向性

$S(N)$を直接計算するのに$O(N)$かかるなら、$N$より小さい$S(\lfloor N/d\rfloor)$の値だけを使って再帰的に表現できないか考える。$\lfloor N/d\rfloor$が取りうる相異なる値は$O(\sqrt N)$種類しかない。

ヒント2: アプローチ

$\sum_{d\mid n}\varphi(d)=n$（ディリクレ畳み込み $\varphi*\mathbf{1}=\mathrm{Id}$）を$n=1,\dots,N$で足し合わせ整理すると $\sum_{d=1}^N S(\lfloor N/d\rfloor)=N(N+1)/2$。これを解いて $S(N)=\frac{N(N+1)}{2}-\sum_{d=2}^N S(\lfloor N/d\rfloor)$。右辺は$\lfloor N/d\rfloor$の値ごとにブロック化して$O(\sqrt N)$項で計算する。小さい$n$は線形篩で前計算し、大きい$n$はメモ化再帰する。

ヒント3: 誘導（コード骨格）

memo = {}
def S(n):
    if n <= N0:
        return prefix[n]
    if n in memo:
        return memo[n]
    res = n * (n + 1) // 2 % MOD
    d = 2
    while d <= n:
        nd = n // d
        d2 = n // nd
        res -= (d2 - d + 1) * S(nd)
        d = d2 + 1
    memo[n] = res % MOD
    return memo[n]

n*(n+1)//2のmod計算には$2$の逆元 pow(2, MOD-2, MOD) を使う。

模範解答 (Python)

import sys

MOD = 998244353
INV2 = pow(2, MOD - 2, MOD)


def solve():
    N = int(sys.stdin.readline())

    N0 = int(N ** (2 / 3)) + 10
    if N0 > N:
        N0 = N

    phi = list(range(N0 + 1))
    primes = []
    is_comp = [False] * (N0 + 1)
    for i in range(2, N0 + 1):
        if not is_comp[i]:
            primes.append(i)
            phi[i] = i - 1
        for p in primes:
            if i * p > N0:
                break
            is_comp[i * p] = True
            if i % p == 0:
                phi[i * p] = phi[i] * p
                break
            else:
                phi[i * p] = phi[i] * (p - 1)

    prefix = [0] * (N0 + 1)
    for i in range(1, N0 + 1):
        prefix[i] = (prefix[i - 1] + phi[i]) % MOD

    memo = {}

    def S(n):
        if n <= N0:
            return prefix[n]
        if n in memo:
            return memo[n]
        res = (n % MOD) * ((n + 1) % MOD) % MOD * INV2 % MOD
        d = 2
        while d <= n:
            nd = n // d
            d2 = n // nd
            res = (res - (d2 - d + 1) * S(nd)) % MOD
            d = d2 + 1
        res %= MOD
        memo[n] = res
        return res

    print(S(N) % MOD)


solve()

計算量: $O(N^{2/3})$（$N_0\approx N^{2/3}$の前計算 + メモ化再帰）。

Step-by-Step 解説

1閾値$N_0$の決定と線形篩
$N_0\approx N^{2/3}$以下は線形篩で$\varphi$を直接計算し前綴和テーブルを作る。

2再帰関数$S(n)$の設計
$n\le N_0$なら前計算済みテーブルを参照。そうでなければ再帰式を使い、メモ化で重複計算を防ぐ。

3ブロックによる和の高速化
$\lfloor n/d\rfloor$が同じ値を取る区間をまとめて処理し、区間数を$O(\sqrt n)$に抑える。

4計算量の見積もり
$N_0$を適切に選ぶことで総計算量が$O(N^{2/3})$に収まる。

5mod計算の注意点
$2$の逆元を使い、割り算をmod上で正しく扱う。

よくあるミス

ミス	原因	正しい書き方
整数除算のまま`n*(n+1)//2 % MOD`とする	mod計算に統一しておらず他言語移植や更なる高速化に弱い	`(n%MOD)((n+1)%MOD)%MODINV2%MOD`を使う
$N_0$を$N^{1/2}$など小さく取りすぎる	前計算は速くなるが再帰の呼び出し回数が増えかえって遅くなる	$N_0\approx N^{2/3}$程度に設定する
`memo`に負の値を保存する	引き算の結果を正規化せず保存	保存前に必ず`% MOD`する
線形篩のループ条件を誤り配列外参照	`if i*p > N0: break`を篩の内側ループ先頭に置き忘れる	内側ループ先頭で必ず範囲チェックする

次のステップ

発展: 同じ枠組みでメビウス関数$\mu$の前綴和$M(N)$を求める（$\mu*\mathbf{1}=\varepsilon$から同様の再帰式）
Min-25篩（Day014既出）との使い分け（素数上の値が多項式的ならMin-25、畳み込みが単純なら杜教篩）を整理する
次回予告: Held-Karpの1-木下界とLagrange緩和（劣勾配法によるTSP下界の逐次改善）

自己評価

完全理解おおむね理解要復習

Day 100-Q4 — 杜教篩（Du's Sieve）

問題