AGC038-C の高速ゼータ変換と高速メビウス変換解法
以下にわかりやすい記事があった。うまく理解できず苦戦した。今もわかっていないが今後のために行間として考えたことをメモ
GCD(Ai, Aj) = xであるような全てのAi, AjについてAi * Aj の和を求めることができれば解が求まることが上記記事よりわかる。
まず、f(x)を(x=Aiであるようなiの数) * Ai と定義する。実装としては1からMax(Ai)までのxについてf(x)をそれぞれ求めるだけで準備できる。 そうすると、f(x) * f(y) は Ai=x, Aj=y となるような全ての(Ai, Aj)に関するAi*Ajの和になる。
次にXの倍数に関する高速ゼータ変換によってF(X)を求める。 F(X)はf(x)の和の形になる。例えばF(1)はf(1)+f(2)+...+f(Max(Ai))であり、F(3)はf(3)+f(6)+f(9)... となる。 これは普通にX<=Max(Ai)について計算すると計算量はMax(Ai) * log(Max(Ai))になる。 高速ゼータ変換を利用すると、Max(Ai) * loglog(Max(Ai)) になる？
次にH(X) = F(X) * F(X) を計算する。単にそれぞれのXについてF(X) * F(X)を計算するだけなので計算量はMax(Ai)
H(X) を高速メビウス変換して、h(x)を得る。 高速メビウス変換は高速ゼータ変換の逆変換であり、h(x) は xの約数Xについてu(x, X) * H(X)の和を取ったものとなる。u(x, X)はメビウス関数という簡単に求められるらしいがよくわかっていない。 h(x)はgcd(x1, x2) = xとなる全てのx1, x2についてf(x1) * f(x2)の和を取ったものになる。 なぜそうなるかは冒頭記事参照。 x<=Max(gcd(x1, x2))についてh(x)を計算すればいい。計算時間はMax(Ai) * loglog(Max(Ai)) になる？
これでh(x)を見ればGCD(Ai, Aj) = xであるような全てのAi, AjについてAi * Aj の和がわかり解が求まる。
わかっていない事として、高速化の方法、なぜ逆変換できるのか(できる条件とナイーブな計算方法)、GCDやLCMじゃないもの一般の集合などについて応用例はどんなものかなどがある。 特に前半2つについてわからないと、まだこの問題は解けない。（とりあえずライブラリコピペしてしまう事で一応解けるけど。）
FFTやNTTとも関連が深そうだし研究も進んでそう。いい教科書が欲しい。
