問題
長さ $N$ の数列 $A_1,\dots,A_N$ と $B_1,\dots,B_N$ が与えられる。
$$C_k = \sum_{\substack{1 \le i,j \le N \\ \mathrm{lcm}(i,j) = k}} A_i B_j \pmod{998244353}$$
を、$k=1,2,\dots,N$ について求めよ（この畳み込みを LCM畳み込み と呼ぶ）。
入力形式
N
A_1 A_2 ... A_N
B_1 B_2 ... B_N制約
$1 \le N \le 10^6$
$0 \le A_i, B_i < 998244353$
入出力例
入力例1
4
1 1 0 0
1 0 1 0出力例1
1 1 1 0入力例2
2
1 0
0 1出力例2
0 1$A_1=1$のみ非零、$B_2=1$のみ非零。寄与は$(i,j)=(1,2)$の1通りのみで$\mathrm{lcm}(1,2)=2$、よって$C_1=0,\ C_2=1$。
概念図: 約数方向のゼータ変換とメビウス変換
ヒント（段階的開示）
ヒント1: 方向性
$\gcd$ 畳み込み（$\sum_{\gcd(i,j)=k}A_iB_j$）は「$k$の倍数方向」への総和（ゼータ変換）→ 各点積 → メビウス逆変換、という手順で高速化できることが知られている。LCM畳み込みも同じ枠組みに載せられないか、$\mathrm{lcm}(i,j)=k$ という条件を「$i,j$ がともに $k$ の約数である」という条件に言い換えられないか考えてみる。
ヒント2: アプローチ
もし $i \mid k$ かつ $j \mid k$ であれば、必ず $\mathrm{lcm}(i,j) \mid k$ が成り立つ。つまり「$i,j$ がともに $k$ の約数」という集合を全て足し合わせると、それは「$\mathrm{lcm}(i,j)=d$ を満たす全ての $d\mid k$ についての $C_d$ の和」に一致する。これは通常の約数総和のゼータ変換（$F[k]=\sum_{d\mid k}A_d$）そのものであり、素数ごとに $O(N\log\log N)$ で計算できる。
ヒント3: 誘導（コード骨格）
# F[k] = sum_{d|k} A[d] を計算する高速ゼータ変換（約数方向）
# 各素数 p ごとに、k = p, 2p, 3p, ... の順（昇順）で F[k] += F[k//p]
#
# G も同様に計算し、H[k] = F[k] * G[k] とすると
# H[k] = sum_{i|k, j|k} A_i B_j = sum_{d|k} C[d]
#
# H から C を復元するには「約数総和の逆変換」＝メビウス変換を行う。
# 各素数 p ごとに、k = (Nをpで割った最大の倍数), ..., 2p, p の順（降順）で H[k] -= H[k//p]模範解答 (Python)
import sys
MOD = 998244353
def sieve_primes(n):
is_composite = [False] * (n + 1)
primes = []
for i in range(2, n + 1):
if not is_composite[i]:
primes.append(i)
for j in range(i * i, n + 1, i):
is_composite[j] = True
return primes
def zeta_divisor(A, N, primes):
# A[k] <- sum_{d | k} A[d]（昇順に多倍数へ伝播）
for p in primes:
for k in range(p, N + 1, p):
A[k] = (A[k] + A[k // p]) % MOD
def mobius_divisor(A, N, primes):
# zeta_divisor の逆変換（降順に多倍数から差し引く）
for p in primes:
start = (N // p) * p
for k in range(start, p - 1, -p):
A[k] = (A[k] - A[k // p]) % MOD
def solve():
data = sys.stdin.read().split()
N = int(data[0])
A = [0] + [int(x) % MOD for x in data[1:1 + N]]
B = [0] + [int(x) % MOD for x in data[1 + N:1 + 2 * N]]
primes = sieve_primes(N)
zeta_divisor(A, N, primes)
zeta_divisor(B, N, primes)
H = [0] * (N + 1)
for k in range(1, N + 1):
H[k] = A[k] * B[k] % MOD
mobius_divisor(H, N, primes)
print(*H[1:N + 1])
solve()計算量: エラトステネスの篩 $O(N\log\log N)$、ゼータ・メビウス変換それぞれ $O(N\log\log N)$。入力例1・stress test（200ケースの乱択ブルートフォース比較）で検証済み。
Step-by-Step 解説
1なぜ「約数」方向のゼータ変換なのか
$\gcd$畳み込みでは「$k$の倍数」方向を使うが、LCM畳み込みでは向きが逆。$i\mid k,\ j\mid k \Rightarrow \mathrm{lcm}(i,j)\mid k$ という包含関係を使うため「$k$の約数」方向のゼータ変換を使う。
$\gcd$畳み込みでは「$k$の倍数」方向を使うが、LCM畳み込みでは向きが逆。$i\mid k,\ j\mid k \Rightarrow \mathrm{lcm}(i,j)\mid k$ という包含関係を使うため「$k$の約数」方向のゼータ変換を使う。
2素数ごとのSOS DP的伝播
$F[k]=\sum_{d\mid k}A_d$ を素直に計算すると $O(N^2)$。素因数分解の構造を使い、素数 $p$ を1つずつ処理して押し込んでいくと各素数で $O(N/p)$ の更新のみで済み、全体 $O(N\log\log N)$ になる。
$F[k]=\sum_{d\mid k}A_d$ を素直に計算すると $O(N^2)$。素因数分解の構造を使い、素数 $p$ を1つずつ処理して押し込んでいくと各素数で $O(N/p)$ の更新のみで済み、全体 $O(N\log\log N)$ になる。
3点ごとの積
$F,G$ が求まったら $H[k]=F[k]G[k] \bmod p$。これは $\sum_{d\mid k}C[d]$ を意味する。
$F,G$ が求まったら $H[k]=F[k]G[k] \bmod p$。これは $\sum_{d\mid k}C[d]$ を意味する。
4メビウス変換で復元
素数ごとに降順に $H[k]\mathrel{-}=H[k/p]$ を適用する（昇順だと二重に引きすぎる）。
素数ごとに降順に $H[k]\mathrel{-}=H[k/p]$ を適用する（昇順だと二重に引きすぎる）。
5計算量の確認
$N\le10^6$ でも $O(N\log\log N)$ なら余裕を持って高速に動作する。
$N\le10^6$ でも $O(N\log\log N)$ なら余裕を持って高速に動作する。
よくあるミス
| ミス | 原因 | 正しい書き方 |
|---|---|---|
| メビウス変換のループを昇順で回してしまう | ゼータ変換と同じ順序で良いと錯覚する | メビウス変換は必ず降順（大きい$k$から）で処理する |
| $\gcd$畳み込みと同じ「倍数方向」ゼータ変換を使ってしまう | LCMとGCDの向きの違いを混同する | LCM畳み込みは「約数方向」のゼータ・メビウス変換を使う |
メビウス変換の開始点を単純にrange(N,p-1,-p)としてしまう | $N$自体が$p$の倍数とは限らないことを見落とす | start = (N // p) * pとして$N$以下で最大の$p$の倍数から始める |
| $A_i,B_i$の添字を0始まりのまま扱ってしまう | 数論的な畳み込みは1始まりが自然であることを忘れる | 配列の先頭にダミー要素を足して1始まりにする |
次のステップ
- 発展: GCD畳み込み（倍数方向のゼータ・メビウス変換）と組み合わせて、$\gcd$と$\mathrm{lcm}$両方の条件が絡む数え上げ問題を解いてみる
- 次回予告: 二部グラフの最小頂点被覆（König's theorem）