Senの競技プログラミング備忘録

こけた問題を自分用の解説で載せる。けんちょんさんのブログを目指したい。質的にも量的にも。こけた問題だけに限定するけど

高速フーリエ変換の超丁寧な解説! Atcoder社の解説を詳解!

はじめに

 学生最強コンのAが高速フーリエ変換そのものであると聞き、高速フーリエ変換を勉強しました。ネット上のフーリエ級数、変換の資料をあさりつつ、tsutajさん、Atcoder社の公式解説をいろいろ眺めて理解に努めようとしたが、頭が弱いため完全理解には6時間かかりました
 けんちょんさんがネットの集合知をまとめてくれることを敬虔に祈念しつつ、自分の勉強も兼ねて、これから高速フーリエ変換を勉強する方々の助けになれるようにこの記事を書きます。

この記事はAtcoder社の公式の高速フーリエ変換のスライドを見ながら、ページごとの詳解をするスタンスを取っています。スライドを同時に見ることを *強く* 推奨します!
また、スライドの後ろの発展的な話題については触れません。初心者向けの記事である上に僕はまだ理解できてません。

まだ始まってもないのですが、ソースコードを大変参考になったtsutajさん、この解説のベースとなったAtcoderの公式解説さん、理解が進まない頭の悪い僕に一生懸命説明したhotman78さんに強い感謝の気持ちをここで示します。

最低限の数学的知識

  • 数3レベルの複素数
  • 数ⅠAレベルの場合の数(もいらないかも)
  • 数学ⅠAレベルの多項式

これさえわかれば、問題ありません。
また、この記事ではアルゴリズムと実装の理屈が分かって使えればよいので、数学的な厳密性は全く保証しません。

解説

これから、スライドの(Xページの)Y行目という表現を多用します。Y行目は、数式、タイトルも一行として数えます。
例えば、スライドの3ページでは、計7行あり、1行目にタイトル、4行目に数式があります。

1ページ

がんばるぞい!
ところでTypical Contestの難易度差激しすぎませんか?

2ぺージ

この問題の類題は2019年の学生最強コン決勝Aにも出ました。高速フーリエ変換を貼るだけで終わるみたいですが、もっと効率のいい解法が存在します。それについてもいずれかきます。(To do)

3ページ

主菜(A)の価格 iが分かれば、副菜(B)の価格 K - iもわかります。 i K - iが同じ時、Aは A_i種類だけ、Bは B_{K-i}種類だけ、それぞれ自由に選んでよいので、スライドのとおりの式になります。
後は、スライドの式のとおり、 iを1から Kまでループを回せば集計できます。

4, 5, 6ページ

唐突に多項式の係数にしてますが、スライドのとおり、二つの多項式の掛け算にすると、掛けてできた多項式の係数は、なんと先ほどの3ページの kごとの計算や集計の結果となるのだ!
ちなみに、この際 A_0 = 0, B_0 = 0とする。こうすると、主菜や副菜のうち一品だけで Kになるのが答えに入らない。なぜなら、係数の掛け算で定数項が0なので、何を掛けても定数項は増えないし、ほかの項に影響を与えない。
つまり、定石みたいなものである!

 A_0 = 0, A_1 = 1, A_2 = 4, B_0 = 0, B _1 = 5, B_2 = 3
ならば、
 (4*x^2+1*x+0)(3*x^2+5*x+0) = 12*x^4+23*x^3+5*x^2+0*x+0
となる。スライド3の方法と一致することを確認できる。

愚直に書けると、多項式の次数を M Nとすると、 O(MN)の計算量になる。
これを O((M+N){log(M+N)})にするテクニックが、高速フーリエ変換である!

7ページ

スライドにある通りです。

8ページ

スライド7ページのとおりに、うまく n+1個以上の点を選び、そこからうまく補間すれば、掛けた式の係数が分かりそう!

9, 10ぺージ

急に複素数が出てきます!なぜ1のn乗根を持つのかというと、
10ページに書いてある通り

  •  i乗( 0 <= i < n)間では異なる iからは異なる \zeta_n^iが得られる。(3行目)
  • 四行目のあるような式を計算すると、キレイサッパリします!なぜその式を計算するのかは、後々わかります。

つまり、簡単に言うと、 x = \zeta_n^iをn個用意して、多項式に代入する形を取れば、計算上都合がいいのです!

11ページ

9ページに続いて、急によくわからないものを定義し始めました。スライドに書いてある通り、今まで考えてた f(x) x = \zeta_n^iの計 n点を代入すると、 n個の f(\zeta_n^i)になる。
 i番目の f(x)に代入した結果が x^iの係数になるようなまた別の式を作っています。
急にそんなことをやりだす理由としては、こうやるとうまくいく!と今考えたほうがいいです。

ちなみに、この離散フーリエ変換(DFT)の形は僕の見てた参考資料と少し違います。競プロでFFTを扱う分には、フーリエ変換のほかの兄弟も知る必要はおそらくないので、気になる人は参考資料のサイトをご覧下さい。このアルゴリズムを理解する分には、こうやっておくとなんか計算がすごく楽!という理解で構いません。

12ページ

先ほど(11ページで)作ったばっかりの \hat{f}(x)の式は、 f(\zeta_n^i)x^iの足し合わせで表してました。
ここで、大もとで求めたかった f(x)はn次の多項式であることを思い出す。
多項式ならば、11ページでの \hat{f}(x)のように、 x^i * c_iと、係数 c_iを掛けたものを足し合わせた形で書けます(それはそう)

11ページのように、 f(x)もその足し算で書いてみました。この際、それぞれの x^iにかかってる係数はわからない(というかこれが知りたい)ので、ひとまず c_iとスライドの中ではおいてます。

そして、12ページでの式変形の結果、どうやらもともと4行目で奥の()の中の足し算に入っていた c_iを、外にくくりだせたようです。
次のソースコードの中の場所のビフォーアフターと同じです。

for(int i = 0; i < A; i++){
    //今はココ
    for(int j = 0; j < B; j++){
        //もともとココ

    }
}

数学では、一つ外のループに出すことで式変形がうまくいく場合があります。今回もそのようです。
次のページで、なぜ1のn乗根をわざわざ代入したのか、なぜくくりだしたのか、解決します!

13, 14ページ

非実数である x^n = 1を満たす xは、
 x^n - 1 = (x-1)(1+x+x^2+...+x^n-1)を満たすので、 0 = 1+x+x^2+...+x^n-1が成り立ちます。この時の xは実数でない1の n乗根のいずれでも成り立つ。
3行目で j \neq kの時、[\zeta_n^{j - k}]は必ず1のn乗根のうちのどれかになります。
そして、1の非実数のn乗根は先ほど述べた累乗を足した式を必ず満たします。3行目では、奥に入ってるΣでは (\zeta_n^{j-k})^i i = 0つまり (\zeta_n^{j-k})^0 = 1から (\zeta_n^{j-k})^{n- 1})まで足し合わせるので、先ほどのn乗根の満たすべき式の右辺そのものになって、0になります!

 j=kの時は、0乗なので1になり、それが n回足し算されるので、最終的にはスライドのとおり、 n*c_iと、係数の n倍になります。つまり、 nで割れば f(x)の係数がわかる!この変換を離散フーリエ逆変換といいます(inverse DFT)。

いままでの変な式や代入する xを選んだのは、うまく選ぶとこうなるからです。

15ページ

 \hat{f}(x)の定義で、 f(\zeta_n^i)とするところを、 f(x)=g(x) * h(x)で書き直して、 g(\zeta_n^i) * h(\zeta_n^i)としてます。
スライドのとおり、どうやらフーリエ変換 g(x), h(x)を先に変換しても、掛ければ f(x)の変換になります!
この性質は、次頁アルゴリズムで重要な役割を果たします。

16. 17, 18ページ

流れのまとめです。
ここから離散フーリエ変換と逆変換が高速でできれば幸せだなぁーとわかります。そこが時間計算量ボトルネックですので。

フーリエ変換は、最終的にフーリエ変換ソースコード \zeta_n^iを-1倍して、最後に nで割るだけなので、本質的にはフーリエ変換を頭よくやります。

19, 20ページ

 1*x^3+3*x^2+5*x^1+7を、 x^i iの偶奇によって、スライド通りに分けると
 f_0(x)=3*x^2+7 と  f_1(x) = 1*x^3 + 5*xとなります。
 f_1(x) xでくくりだすと、二つの式は x^2=tで置換をしてみると
 F_0(t) = 3*t+7 と  F_1(t) = 1*t+5
とかけます。
スライドの6行目の式のように、 f(x)の値を知るには、どうやら F_0(x^2) F_1(x^2)の値が分かればいいとなります。
この時点ではまだ意味不明のようだが、次の2,3枚スライドで意味が分かります。

そして、 x = \zeta_n^iをすべて n通り代入すると、スライド20ページにある値が必要となります。

21ページ

なんということでしょう!2行目のように式をいじくると、 n乗根のおかげで、 \zeta_{2n}^2 = \zeta_nになります!実をいうと、三角関数ので表すとまったく同じものであるともわかります!

これのおかげで、20ページで求めたかった2つの数列は7~8行目のように変形できます。この変形のうまあじは次のスライドで伝わります。

22ページ

 n乗根の周期性で、周期が n n/2になったので、なんと20ページで求めたかった数列の前半と後半は実は同じだったのです!
例を挙げると
 f(\zeta_4^0) f(\zeta_4^2) f(\zeta_4^4) f(\zeta_4^6)は、
 f(\zeta_2^0) f(\zeta_2^1) f(\zeta_2^2) f(\zeta_2^3) と同じで、
それは
 f(\zeta_2^0) f(\zeta_2^1) f(\zeta_2^0) f(\zeta_2^1) と同じである。そしてこれらは、2つで周期になってる。

23, 24ページ

22ページの性質を利用すると、
 2^n個の値を代入した値を計算するためには、2つの 2^{n-1}の値の数列がわかれば O(n)で計算できます。」
例えば、8個の点を代入した値の列がほしいのなら、4個の点を代入した列を2つ、求めればよい。そうすると、次図のように完全二分木みたいな感じの再帰関数になります。(次図は 1*x^2+2*x-3の離散フーリエ変換の例)
詳しい図の解説は次回の記事でします。

f:id:Sen_comp:20191002094632p:plain

25ページ

具体的な実装となります。僕もこれを実質写経しました。

26~30ページ

見なかったことにしてください。おそらく初心者向きでは全くありません。

最後に

ここまで読んでいただきお疲れさまでした。ここまでで5555字だそうです。
上記の解説だけで、ピンとくる人は少ないと思うので、続編として実装例と、具体的な値を入れた時の挙動の解説をします。

次編-> 
sen-comp.hatenablog.com

何か誤字脱字、論理上のミスなどあったらご指摘お願いします。