連載しょぼんコーダー 2 包除

第2回は, 約数系包除を使うある問題について, 典型的な方法と, 包除原理を用いたもう1通りの解法を解説します.

次の問題を2通りで解いてみます.

例題2（想定Diff 1500）
自然数 $N$ が与えられます. 集合 $\{1, 2, ..., N\}$ の要素の数が $2$ 以上である部分集合 $\{x_1, x_2, ..., x_k\}$ について, $x_1, x_2, ..., x_k$ の最大公約数が $1$ であるものは何通りありますか？ただし, 答えは非常に大きくなる可能性があるので $1000000007$ で割った余りを求めてください.
時間制限： $2$ 秒, メモリ制限： $1$ GB
制約： $2 \leq N \leq 262143$
入力：4　→　出力：10
$\{1, 2, 3, 4\}, \{1, 2, 3\}, \{1, 2, 4\},$
$\{2, 3, 4\}, \{1, 3, 4\}, \{1, 2\},$
$\{2, 3\}, \{3, 4\}, \{1, 3\}, \{1, 4\}$ の $10$ 通りあります.

まずは, DP(動的計画法)を使う方法を考えてみます. 典型テクニックとして,
$DP(i)$ ： $\{1, 2, ..., N\}$ の部分集合であって, 要素の数が $2$ 以上で, そのすべての要素の最大公約数がちょうど $i$ であるものの個数
とします. 1つ目の解法は「約数系包除」と呼ばれているものです.

まず, $\{1, 2, ..., N\}$ の部分集合 $X$ であって, 要素の数が $2$ 以上で, 最大公約数が $i$ の倍数のものを数えることを考えます.
最大公約数が $i$ の倍数になるための必要十分条件は, $X$ の要素がすべて $i$ の倍数になることです. これを示します.
まず, 最大公約数が $i$ の倍数ならば, 最大公約数 $g$ は整数 $m$ を用いて $g = mi$ と表すことができます. 任意の $X$ の要素 $a$ について, $a = ki$ となる整数 $k$ があります. これは, $a$ が $i$ の倍数であることを示します.
次に, $X$ の要素がすべて $i$ の倍数ならば, $i$ は $X$ のすべての要素を割り切ります. よって, 最大公約数は $i$ の倍数になります.

今, $\{1, 2, ..., N\}$ の要素であって, $i$ の倍数であるものは $\lfloor N/i \rfloor$ 個存在します.
そのうち, 要素の数が $2$ 以上の部分集合は, $m = \lfloor N/i \rfloor$ とすると, 要素数 $m$ の部分集合は $2^m$ 個で, まず空集合 $1$ 個を除き, 次に要素数 $1$ の集合 $m$ 個を除きます.
残ったのは $2^m - m - 1$ 個で, これが $\{1, 2, ..., N\}$ の部分集合 $X$ であって, 要素の数が $2$ 以上で, 最大公約数が $i$ の倍数であるものの個数です.

さて, $DP(i)$ を求めます. 今, $DP(i) + DP(2i) + DP(3i) + \cdots = 2^m - m - 1$ ですから, 上から順番に更新することで $DP(i) = 2^m - m - 1 - DP(2i) - DP(3i) - \cdots$ で更新できます.（もちろん, 再帰関数を使ってもよいです.）
このとき, $m=0, 1$ でも偶然かうまくいってます. なので場合分けせずにそのまま計算します.
$j > N \Longrightarrow DP(j) = 0$ なので, $N$ 未満を考慮すればよいので有限回でプログラムを終わらせることができます.

ここで, 各 $i$ について, $2^m$ を計算するのにそれぞれ $O(\log m)$ かかります. $m = \lfloor N/i \rfloor$ で, $1, 2, \cdots, N$ についてこれを計算するので, $O(\sum^{N}_{i=1} \log \lfloor N/i \rfloor) = O(N)$ です. $2i, 3i, \cdots, mi$ を計算すればよいのですが, $O(\sum^{N}_{i=1} \lfloor N/i \rfloor) = O(N \log N)$ です. これが主にボトルネックになります.

よって, 結局 $DP(1)$ を求めるのに $O(N \log N)$ かかります.

mod = 1000000007
n = int(input())
#dp[i] = gcd(X)がiであるようなものの数.
t = n//2
dp = [0 for i in range(t+1)]
for i in range(t,0,-1):
    dp[i] = (pow(2,n//i,mod)-n//i-1)%mod
    for j in range(2*i,t+1,i):
        dp[i] = (dp[i] - dp[j]) % mod
print(dp[1])

次に, 包除原理を使って問題を解きます. 包除原理は, いくつかの部分集合の和集合の要素数について, それらの積集合の要素数を利用して求めるテクニックで, ダブルカウントを防ぎたいときに有利です.
高校数学の数Aで集合の要素をやった方は, 次の式に見覚えがあるかもしれません.

$|A \cup B| = |A| + |B| - |A \cap B|$
$|A \cup B \cup C| = |A| + |B| + |C|$ $- |A \cap B| - |A \cap C| - |B \cap C| - |A \cap B \cap C|$

これで2個の場合は分かるのではないかと思います. しかし一般の場合に拡張しようとすると大変です. 個々のエリアの要素について, その要素を含む集合をすべて着目します.
たとえば, 2個の集合の場合, $A - B$ の要素については, $A$ しか寄与しません . それに対し $A \cap B$ の要素については, $A, B$ の2つが寄与します. 普通に足すと, $A \cap B$ の各要素が重複してカウントされ, 求めた結果がでません.
そこで, $A \cap B$ をマイナスの寄与と考えることで, 寄与が $1+1-1 = 1$ 個分になり, ダブルカウントを防ぐことができます.

いま, $x$ がちょうど $n$ 個の集合 $A_1, A_2, \cdots, A_n$ に含まれる要素であるとして, 寄与をちょうど1個にする方法を考えます.
1個の積集合は $A_1, A_2, \cdots, A_n$ の $n = \binom{n}{1}$ 個です. ただし, $\binom{n}{k} = {}_nC_{k}$ は二項係数です.
2個の積集合は $A_1 \cap A_2$ などの $\frac{n(n-1)}{2} = \binom{n}{2}$ 個です. $n$ 個の集合から $2$ 個選ぶことを考えるとわかりやすいと思います.
同様にして, $k$ 個の積集合は全部で $\binom{n}{k}$ 個です. ここで, $k$ が偶数のときプラスに寄与, $k$ が奇数のときマイナスに寄与すると考えます. このとき, $x$ がカウントされる回数は

$\binom{n}{1} - \binom{n}{2} + \binom{n}{3} + \cdots + (-1)^{k-1} \binom{n}{k} + \cdots + (-1)^{n-1} \binom{n}{n}$

となります. 二項定理より $\sum^{n}_{k=0} (-1)^k \binom{n}{k} = (1-1)^n = 0$ となります.
よって, $\sum^{n}_{k=1} (-1)^{k-1} \binom{n}{k} = - \sum^{n}_{k=0} (-1)^k \binom{n}{k} + \binom{n}{0} = 1$ であり, なんとすべての寄与を合計すると $1$ となります. これが包除原理の根本にある仕組みです.

これで, $|A \cup B \cup C|$ の式で $|A \cap B \cap C|$ の寄与がプラスである理由が分かったと思います. 同様に4個の場合の式も作れて,
$|A \cup B \cup C \cup D| = |A| + |B| + |C| + |D|$
$- |A \cap B| - |A \cap C| - |A \cap D|$ $- |B \cap C| - |B \cap D| - |C \cap D| + |A \cap B \cap C|$
$+ |A \cap B \cap D| + |A \cap C \cap D| + |B \cap C \cap D| - |A \cap B \cap C \cap D|$
となります.

そろそろ本題として, もとの問題の答えを求めてみましょう. すべての要素数 $2$ 以上の部分集合の個数から, 最大公約数が $2$ 以上であり, 要素数 $2$ 以上の部分集合の個数を除けばよいです.
後者を数え上げることを考えてみましょう. 1は素数を約数に持ちませんが, 2以上の整数は必ず素数を約数に持ちます. 最大公約数が素数 $p$ の倍数（すなわち, 要素すべてが $p$ の倍数）であり, 要素数 $2$ 以上の部分集合すべての集合を $A_p$ とします.
求めたいのは $|A_2 \cup A_3 \cup A_5 \cup A_7 \cup A_{11} \cup \cdots|$ です. $p > N$ なら, $A_{p} = \emptyset$ ですので, 実際は有限個の和集合となります.

実際に包除原理を愚直に適用すると, $N$ 以下の最大の素数を $p_m$ , 素数の個数を $\pi(N)$ として, $\{A_2 \cup A_3 \cup A_5 \cup A_7 \cup A_{11} \cup \cdots \cup A_{p_m}\}$ の空でない部分集合をすべて考える必要があるので, $O(2^\pi(N))$ かかり間に合いません. TLEコードは以下です. 前半ではエラトステネスのふるいをしていて, 後半ではDFSのように探索しています.

n = int(input())

sieve = []
flag = [0] * (n+1)

for i in range(2, n+1):
    if flag[i] == 0:
        sieve.append(i)
        for j in range(i, n+1, i):
            flag[j] = 1

mod = 10**9 + 7
ans = (pow(2, n, mod) - n - 1) % mod

size = len(sieve)
mada = [(1, 0, 0, 1)] #(現在の値, sieveをここまで見た, 探索するか, パリティ)
while mada:
    i, ind, tosearch, parity = mada.pop()
    if tosearch:
        ans = (ans + parity * (pow(2, n//i, mod) - n//i - 1)) % mod
    if ind < size:
        mada.append((i * sieve[ind], ind + 1, 1, - parity))
        mada.append((i, ind + 1, 0, parity))

print(ans)

そこで, $A_{p_1} \cap A_{p_2} \cap \cdots \cap A_{p_k}$ は要素の数が $2$ 以上で, 最大公約数が $p_1p_2 \cdots p_k$ の倍数であるものの個数です. よって, $p_1p_2 \cdots p_k > N$ の場合は $|A_{p_1} \cap A_{p_2} \cap \cdots \cap A_{p_k}| = 0$ となり, 考慮する必要がありません.
よって, $p_i$ の積集合をだんだん探索していき, その積が $N$ を超えた瞬間にストップすることを考えます. これを枝刈りといいます.
また, $p_i \lt p_j$ で, 正整数 $A$ に対し, $N \lt Ap_i$ ならば $N \lt Ap_j$ が成立します. よって積が $N$ を超えた瞬間,「 $p_i$ を選択しないという選択肢」も排除可能です. これで時間計算量をさらに短縮することができます.

これを踏まえると, 次のようになります. 時間計算量はわかりませんが, エラトステネスのふるいで $O(N \log \log N)$ , DFSでおそらく $O(N \log N)$ 程度です.

n = int(input())

sieve = []
flag = [0] * (n+1)

for i in range(2, n+1):
    if flag[i] == 0:
        sieve.append(i)
        for j in range(i, n+1, i):
            flag[j] = 1

mod = 10**9 + 7
ans = (pow(2, n, mod) - n - 1) % mod

size = len(sieve)
mada = [(1, 0, 0, 1)] #(現在の値, sieveをここまで見た, 探索するか, パリティ)
while mada:
    i, ind, tosearch, parity = mada.pop()
    if tosearch:
        ans = (ans + parity * (pow(2, n//i, mod) - n//i - 1)) % mod
    if ind < size:
        if i * sieve[ind] <= n:
            mada.append((i * sieve[ind], ind + 1, 1, - parity))
            mada.append((i, ind + 1, 0, parity))

print(ans)