Day 045-Q2 — 格子上ゼータ・メビウス変換(GCD集計)

2026-05-29 赤色 Master / Phase 8+ ★★★★★★★★★ Möbius Inversion + 線形篩 + GCD集計

問題

$N$ 個の正整数 $a_1, \dots, a_N$($1 \le a_i \le M$)が与えられる。 全ての組 $(i, j)$($i < j$)に対して $\gcd(a_i, a_j) = g$ となる組の個数 $f(g)$ を、 全ての $g = 1, 2, \dots, M$ について同時に求めよ。

制約

$1 \le N \le 2 \times 10^5$
$1 \le M \le 10^6$
$1 \le a_i \le M$
時間制限: 2秒
$O(M \log M)$ で解くこと

入出力例

入力例 1

5 10
2 4 6 4 2

出力例 1

g=1: 0
g=2: 4
g=3: 0
g=4: 1
g=5: 0
g=6: 1
g=7: 0
g=8: 0
g=9: 0
g=10: 0

概念図: ゼータ変換とメビウス反転のフロー

freq[g] a_i=g の個数 ゼータ変換 O(M log M) cnt[g] g の倍数の個数 二項係数 C(cnt,2) h[g] g の倍数GCDの組数 メビウス反転 O(M log M) f[g] GCD=g の組数 メビウス反転公式: h[g] = Σ_{k≥1} f[g·k] (g の倍数を GCD として持つ組の合計) 逆に: f[g] = Σ_{k≥1} μ(k) · h[g·k] (メビウス関数 μ を使った反転) μ(1)=1, μ(p1·p2···pk)=(-1)^k(相異なる素数積), μ(n)=0(二乗因子あり) μ は線形篩で O(M) 計算 → 全体 O(M log M)

ヒント(段階的開示)

ヒント1: 方向性
$f(g)$ を直接数えるのではなく、まず「$g$ の倍数である $a_i$ の個数 $cnt[g]$」を数え、 $g$ の倍数全体の組数から、より大きい公約数を持つ組を引くメビウス反転を使います。
ヒント2: アプローチ
  1. $cnt[g]$ = $a_i$ が $g$ の倍数であるものの個数(ゼータ変換)
  2. $h(g) = \binom{cnt[g]}{2}$ = 公約数が $g$ の倍数である組の数
  3. $f(g) = \sum_{k \ge 1} \mu(k) \cdot h(gk)$(メビウス反転公式)
ヒント3: 実装骨格
# 線形篩でμを計算 O(M)
# ゼータ変換: cnt[g] = Σ_{g|k} freq[k]  O(M log M)
for g in range(1, M+1):
    for k in range(g, M+1, g):
        cnt[g] += freq[k]
# h[g] = C(cnt[g], 2)
# メビウス反転: f[g] = Σ_k μ(k) * h[g*k]  O(M log M)
for g in range(1, M+1):
    for k in range(1, M//g+1):
        if mu[k] != 0:
            f[g] += mu[k] * h[g*k]

模範解答 (Python)

import sys

def main():
    data = sys.stdin.read().split()
    idx = 0
    N, M = int(data[idx]), int(data[idx+1]); idx += 2
    a = [int(data[idx+i]) for i in range(N)]; idx += N

    # メビウス関数の線形篩
    mu = [0] * (M + 1)
    mu[1] = 1
    primes = []
    is_composite = [False] * (M + 1)
    for i in range(2, M + 1):
        if not is_composite[i]:
            primes.append(i)
            mu[i] = -1
        for p in primes:
            if i * p > M:
                break
            is_composite[i * p] = True
            if i % p == 0:
                mu[i * p] = 0
                break
            else:
                mu[i * p] = -mu[i]

    freq = [0] * (M + 1)
    for x in a:
        freq[x] += 1

    # ゼータ変換: cnt[g] = Σ freq[k] for k divisible by g
    cnt = [0] * (M + 1)
    for g in range(1, M + 1):
        for k in range(g, M + 1, g):
            cnt[g] += freq[k]

    h = [cnt[g] * (cnt[g] - 1) // 2 for g in range(M + 1)]

    # メビウス反転: f[g] = Σ μ(k) * h[g*k]
    f = [0] * (M + 1)
    for g in range(1, M + 1):
        for k in range(1, M // g + 1):
            if mu[k] != 0:
                f[g] += mu[k] * h[g * k]

    out = []
    for g in range(1, M + 1):
        out.append(f"g={g}: {f[g]}")
    print('\n'.join(out))

main()

Step-by-Step 解説

1頻度配列とゼータ変換
$freq[x]$ = $a_i = x$ の個数。$cnt[g] = \sum_{g|k} freq[k]$ は「$g$ の倍数を値として持つ要素の個数」。外ループ $g$、内ループ $g$ の倍数 $k$ で $O(M \log M)$ で計算できる(調和級数)。
2二項係数 $h[g]$
$h[g] = \binom{cnt[g]}{2}$ は「GCDが $g$ の倍数になる組の数」(GCDが厳密に $g$ かはまだ不明)。
3メビウス関数の線形篩
線形篩で $\mu(n)$ を $O(M)$ で計算。$\mu(n) = 0$ の項は計算をスキップして高速化。
4メビウス反転
$f[g] = \sum_{k \ge 1} \mu(k) \cdot h[gk]$ により、「GCDがちょうど $g$ になる組の数」が $O(M \log M)$ で全 $g$ について一斉に求まる。
5全体計算量
線形篩: $O(M)$、ゼータ変換: $O(M \log M)$、メビウス反転: $O(M \log M)$。合計 $O(M \log M)$。

計算量

線形篩: $O(M)$
ゼータ変換: $O(M \log M)$(調和級数)
メビウス反転: $O(M \log M)$
空間: $O(M)$

よくあるミス

ミス原因正しい書き方
μ(n)=0 の項を計算不要な乗算if mu[k] != 0: でスキップ
cnt 配列を 0-indexedg=1 の添字がずれるcnt = [0] * (M + 1) で1-indexed
ゼータ変換の方向間違い因数方向と倍数方向を逆に外ループ g, 内ループ gの倍数 k
h[0] を計算してしまう範囲外参照range(1, M+1) から開始

次のステップ

  • 発展問題: $\sum_{i=1}^{N} \sum_{j=1}^{N} \gcd(a_i, a_j)$ を $O(M \log M)$ で計算する
  • 類題: LCM convolution、AND/OR/XOR convolution との関係
  • 応用: 数論関数の Dirichlet 積と高速計算

自己評価