ABC 353 / Convolution (予習)

May 12, 2024

ABC 353

ABC 353 に参加しました。何とかレーティングをキープできました。

A 問題

y+0 のように +0 を入れると read が整数型に推論される技を知りました:

main=interact$show.f.map read.words;f(_:h:x)=head$[i|(i,y)<-zip[2..]x,y+0>h]++[-1]

B 問題

% 演算子を定義し、関数呼び出し時の空白の数を減らしましょう:

main=interact$show.(0%).tail.map read.words
a%[_]=min 1a
a%(k:y:r)|a+y>k=1+0%(k:y:r)|0<1=(a+y)%(k:r)

C 問題

AA をソートします。 mod\bmod と言わず 10810^8 を引くことにすれば、一括処理できます。

ii を固定したとき、和が 10810^8 を超え始める点を見つけるためには 2 分探索または尺取り法を使います。 Haskell で実装する尺取り法の必勝法……を考え始めることの検討を加速します。

D 問題

ゴリ押しするか簡単な解法を考え抜きます。

私はゴリ押ししました。日記もゴリ押しに一塩です。

E 問題

これは輪を通る弦の組み合わせの数をカウントする問題だと思いました。先頭の文字でグループ分けしつつ、 1 文字ずつ剥いて行けば解けます。

2024-05-12-rings.webp
Figure 1: 貪欲に数えて行けば良い

なお文字列ソートの計算量は O(NL1(logN+logL1))O(N L_1 (\log N + \log L_1)) 程度 らしいです 。今回は文字列長の総和 L3105L \le 3 \cdot 10^5 のため N=L1=LN = L_1 = \sqrt L とすれば O(LlogL)O(L \log L) ? 速いですね〜 (読もう……)

G 問題

EDPC - Z と同じ問題に見えました。絶対値を外すため、 CHT に 2 本ずつ直線を追加していけば解けるはずです。 CHT は未習得の高度典型なので解けませんでした。

ところが絶対値を外して現れた式は、単なる 2 つの Max であってセグメント木 2 本に載ります。解けるべき問題でした。

なお CHT に似た何かがある? みたいなので、 CHT に挑む時までに gksato さんの提出を読み込んでみようと思います。

Convolution (予習)

高度典型: FFT/NTT/Convolution の勉強ノートです (未実装) 。

内容は誤りだらけです 。真剣に読まないでください。

三角基底による関数の展開

まずはフーリエ変換をエアプします。 フーリエ解析―基礎と応用 を参考にしました。

ボトムアップかつ天下り的に行きます。

冪級数展開 (Taylor 展開)

三角関数のように無限回微分できる関数は、冪級数の形に展開できることが知られています:

f(x)=n=0f(n)(x0)n!(xx0)nx0=0=n=0f(n)(0)n!xn\begin{aligned} f(x) &= \sum \limits_{n=0}^{\infty} \frac {f^{(n)}(x_0)} {n!} (x - x_0)^{n} \bigg|_{x_0 = 0} \\ &= \sum \limits_{n=0}^{\infty} \frac {f^{(n)}(0)} {n!} x^{n} \end{aligned}

三角関数と指数関数を x0=0x_0 = 0 の回りで冪級数展開すると、オイラーの公式 eix=cosx+isinxe^{ix} = \cos x + i \sin x が得られます。

三角基底 (複素三角基底?) {eix}\{e^{ix}\} は完全直行基底であり、ある種の f(x)f(x) が三角基底の線型結合で表せることが知られています。

フーリエ級数 (FS)

フーリエ変換の基本は周期関数です。周期 TT の周期関数 f(x)f(x) は周期 T,2T,T, 2T, \dots の三角基底の線型結合で表せるものとします。係数 cnc_n を複素フーリエ級数と呼びます:

f(x)n=cnei2nπTx:=n=cneiωnx\begin{aligned} f(x) &\sim \sum_{n=-\infty}^{\infty} c_n e^{i \frac {2 n \pi} {T} x} \\ :&= \sum_{n=-\infty}^{\infty} c_n e^{i \omega_n x} \end{aligned}

三角基底は直行基底であり、基底に対する成分 (フーリエ級数 cnc_n) は内積に比例します:

(f(x),eiωx)=(cneiωx,eiωx)=cn(eiωx,eiωx)=cneiωx2=cnT2T2eiωxeiωxdx=cnT\begin{aligned} (f(x), e^{i\omega x}) &= (c_n e^{i \omega x}, e^{i \omega x}) \\ &= c_n (e^{i \omega x}, e^{i \omega x}) \\ &= c_n \| e^{i \omega x} \|^2 \\ &= c_n \int_{\frac{-T}{-2}}^{\frac{T}{2}} e^{i\omega x} e^{-i\omega x} \mathrm{d}x \\ &= c_n T \end{aligned}

cnc_nf(x)f(x) に代入すると以下の形になります:

f(x)=n=1T(f(x),eiωnx)eiωnxf(x) = \sum_{n=-\infty}^{\infty} \frac {1} {T} (f(x), e^{i\omega_n x}) e^{i \omega_n x}

フーリエ変換 (FT)

T=2πΔωT = \frac {2\pi} {\Delta \omega} \rightarrow \infty とすれば、非周期関数 f(x)f(x) を連続な基底ベクトル {eiωx}ω\{e^{i\omega x}\}_{\omega} で展開できます (逆フーリエ変換 (IFT)) 。

f(x)=n=Δω2π(f(x),eiωnx)eiωnx=dω2π(f(x),eiωx)eiωx:=dω2πF[f](ω)eiωx\begin{aligned} f(x) &= \sum_{n=-\infty}^{\infty} \frac {\Delta \omega} {2 \pi} (f(x), e^{i\omega_n x}) e^{i\omega_n x} \\ &= \int_{-\infty}^{\infty} \frac {\mathrm{d}\omega} {2 \pi} (f(x), e^{i\omega x}) e^{i\omega x} \\ :&= \int_{-\infty}^{\infty} \frac {\mathrm{d}\omega} {2 \pi} \mathcal{F}[f](\omega) e^{i\omega x} \end{aligned}

フーリエ変換 F[f](ω):=(f(x),eiωx)\mathcal{F}[f](\omega) := (f(x), e^{i\omega x})f(x)f(x)eiωxe^{i\omega x} 成分に相当します。指数関数の微積分は簡単なので、 {eiωnx}\{e^{i\omega_n x}\} を基底としたのは幸先良さそうです。

三角基底による離散関数の展開

デジタル信号処理へ寄ります。ここからは ビギナーズ デジタルフーリエ変換 および小野測器の 計測コラムem138 添付資料 を参考にエアプします。

時間離散フーリエ変換

関数 f(x)f(x)τ\tau の間隔で離散化し、離散関数 fτ(x)f_{\tau}(x) を得ます。 f(x)f(x) に周期的デルタ関数 δτ(x)\delta_{\tau}(x) をかけることで、 f(x)f(x) の離散化を表現します:

δτ(x):=nZδ(xnτ)fτ(x):=δτ(x)f(x)\begin{aligned} \delta_{\tau} (x) &:= \sum_{n \in \mathbb{Z}} \delta(x - n \tau) \\ f_{\tau}(x) &:= \delta_{\tau} (x) f(x) \end{aligned}

fτ(x)f_{\tau}(x) のフーリエ変換も離散化されています:

F[fτ](ω)=(fτ(x),eiωx)=fτ(x)eiωxdx=nNf(nτ)eiωnτ\begin{aligned} \mathcal{F}[f_{\tau}](\omega) &= (f_{\tau}(x), e^{i\omega x}) \\ &= \int_{-\infty}^{\infty} f_{\tau}(x) e^{-i\omega x} \mathrm{d}x \\ &= \sum_{n \in \mathbb{N}} f(n \tau) e^{-i \omega n \tau} \end{aligned}

さらに時間が離散化されたことから位相に周期性が生じており、 F[fτ(x)](ω)\mathcal{F}[f_{\tau}(x)](\omega) は周期 2πτ\frac {2 \pi} {\tau} の周期関数になりました (F[fτ](ω)=F[fτ](ω+2πτ)\mathcal{F}[f_{\tau}](\omega) = \mathcal{F}[f_{\tau}](\omega + \frac {2 \pi} {\tau}))。

改めて fτ(nτ)f_{\tau}(n \tau) を三角基底で展開すると、次の式を得ます:

fτ(nτ)=πτπτdω2πF[fτ](ω)eiωnτ\begin{aligned} f_{\tau}(n \tau) &= \int_{-\frac {\pi} {\tau}}^{\frac {\pi} {\tau}} \frac {\mathrm{d}\omega} {2 \pi} \mathcal{F}[f_{\tau}](\omega) e^{i\omega n \tau} \end{aligned}

離散フーリエ級数 (DFS)

時間信号 x(t)x(t) の一部を時間幅 TT で切り取り、 τ\tau の間隔で離散化したとします。この信号列を周期 TT の周期関数に拡張すると、やはり複素フーリエ級数の和の形に展開できます。

x(kτ)=nZcneiωnkτ=nZcnei2nπTkτ\begin{aligned} x(k \tau) &= \sum_{n \in \mathbb{Z}} c_n e^{i \omega_n k \tau} \\ &= \sum_{n \in \mathbb{Z}} c_n e^{i \frac {2n\pi} {T} k \tau} \end{aligned}

nZn \in \mathbb{Z} とありますが、時間が離散化されたことによって位相が周期的になり、 x(kτ)x(k \tau){eiωnx}n[0,N1]Z\{e^{i\omega_n x}\}_{n \in [0, N - 1] \cap \mathbb{Z}} のみにより展開されるはずです。実際 nn+rNn \rightarrow n + rN に分解して確かめられます (天才だ……):

x(kτ)=n[0,N1]ZrZcn+rNeiωn+rNkτ=n[0,N1]ZeiωnkτrZcn+rNcn=0(n>N)=n[0,N1]Zeiωnkτcn\begin{aligned} x(k \tau) &= \sum_{n \in [0, N - 1] \cap \mathbb{Z}} \sum_{r \in \mathbb{Z}} c_{n + rN} e^{i \omega_{n + rN} k \tau} \\ &= \sum_{n \in [0, N - 1] \cap \mathbb{Z}} e^{i \omega_{n} k \tau} \sum_{r \in \mathbb{Z}} c_{n + rN} \bigg|_{c_n = 0 (|n| > N)} \\ &= \sum_{n \in [0, N - 1] \cap \mathbb{Z}} e^{i \omega_{n} k \tau} c_n \end{aligned}

上記の 2 ~ 3 行で帯域制限を設け、 x(t)x(t) が N 次の高調波成分までしか持たないとした場合、 rZcn+rN=cn\sum_{r \in \mathbb{Z}} c_{n + rN} = c_n から非常に簡素な式に整理できました。 N 個 (N=TτN = \frac T {\tau}) の信号 {x(nτ)}n\{x(n\tau)\}_n を整理すると、以下の行列で書けます:

[x(0)x(τ)x(TτT)]=[1111eiω1eiω1(N1)1eiωN1eiωN1(N1)][c0c1cN1]\begin{aligned} [x(0)x(τ)x(TτT)]\begin{bmatrix} x(0) \\ x(\tau) \\ \vdots \\ x(\frac {T - \tau} T) \end{bmatrix} &= \begin{bmatrix} 1 & 1 & \dots & 1 \\ 1 & e^{i\omega_1} & \dots & e^{i{\omega_1(N-1)}} \\ \vdots & \vdots & & \vdots \\ 1 & e^{i\omega_{N-1}} & \dots & e^{i{\omega_{N-1}(N-1)}} \\ \end{bmatrix} [c0c1cN1]\begin{bmatrix} c_0 \\ c_1 \\ \vdots \\ c_{N-1} \end{bmatrix} \end{aligned}

よって時間信号 x(t)x(t) に対して 1. その一部を時間幅 TT で切り抜き 2. τ\tau の間隔で離散化し 3. 周期 TT の周期関数に拡張し 4. 帯域制限を行うと、 N 回のサンプリング結果から複素フーリエ級数が分かり、 AC の新作が発表されます。

離散フーリエ変換 (DFT)

上記の x(nτ):=xnx(n\tau) := x_n に対し、なぜか改めて離散フーリエ変換 XkX_k を定義します (なぜ……?):

Xk:=n[0,N1]Zxnei2πNk:=n[0,N1]ZxnWnk\begin{aligned} X_k &:= \sum_{n \in [0, N - 1] \cap \mathbb{Z}} x_n e^{-i \frac {2\pi} {N}k} \\ &:= \sum_{n \in [0, N - 1] \cap \mathbb{Z}} x_n W^{nk} \end{aligned}

この気持ちは勉強不足のため理解できていません。フーリエ変換の方が重要なので、単純な式にしたかった?

高速フーリエ変換 (FFT)

(2-radix) DFT の再帰的な定義

DFT の高速計算 (FFT) を小野測器の 計測コラム emm140 号用 から学びます。 P2 の図から汲み取れる通り、 8 点 DFT の出力 {X8,k}k[0,7]\{X_{8, k}\}_{k \in [0, 7]} は、それぞれ 4 点 DFT の和に分解できます:

X8,k({xk}k[0,7]Z)=n[0,7]xnWnk=n[0,2,4,6]xnWnk+n[0,2,4,6]xn+1W(n+1)k=n[0,2,4,6]xnWnk+Wkn[1,3,5,7]xnWn=n[0,2,4,6]xnWnk+Wkn[1,3,5,7]xnWn=X4,k(x0,x2,x4,x6)+WkX4,k(x1,x3,x5,x7)\begin{aligned} X_{8,k} (\{x_k\}_{k \in [0, 7] \cap \mathbb{Z}}) &= \sum_{n \in [0, 7]} x_{n} W^{nk} \\ &= \sum_{n \in [0, 2, 4, 6]} x_{n} W^{nk} + \sum_{n \in [0, 2, 4, 6]} x_{n+1} W^{(n+1)k} \\ &= \sum_{n \in [0, 2, 4, 6]} x_{n} W^{nk} + W_k \sum_{n \in [1, 3, 5, 7]} x_{n} W^{n} \\ &= \sum_{n \in [0, 2, 4, 6]} x_{n} W^{nk} + W_k \sum_{n \in [1, 3, 5, 7]} x_{n} W^{n} \\ &= X_{4,k}({x_0, x_2, x_4, x_6}) + W_k X_{4,k}({x_1, x_3, x_5, x_7}) \end{aligned}

8 点 DFT の出力を以下の信号流れ図にまとめます。 {xn}n\{x_n\}_n が上下に 2 分割されており、再帰的に O(NlogN)O(N \log N) で計算できることが予想できます。

2024-05-12-butterfly-8.webp

バタフライ演算の部分を丁寧に図示すると以下です。黒点を接続とし、加算器と乗算器を明示しています:

2024-05-12-butterfly-8'.webp
Figure 2: 信号流れ図が読めなかったので

これを簡略化し、また W84=eiπ=1W_8^4 = e^{i \pi} = -1 を代入すると、前の図になります。読めないよ〜〜

バタフライ演算

8 点 DFT を再帰的に展開すると、 x0,..,x7x_0, .., x_7x0,x4,x2,x6,x1,x5,x3,x7x_0, x_4, x_2, x_6, x_1, x_5, x_3, x_7 の並びになります:

2024-05-12-butterfly-8-4-2-1.webp

したがって以下の手順で高速に DFT を計算できます。

  1. 添字の置換を (一括して) 行う
  2. バタフライ演算を繰り返し適用する

添字のソート方法

人が数列をソートするときは、最も大きな位から順にソートすることが多いです。 (2-radix) FFT においては (2 進数表記で) 小さな位から順番にソートされていくことになります。したがって添字をビット反転 (例: 0b1100 -> 0b0011) した値を基準に {xn}n\{x_n\}_n を昇順ソートしたことになります。

2024-05-12-dft-indices.webp
Figure 3: 2-radix FFT における引数の置換の追跡 (検索すれば同じ図が出てきます)

NTT (数論変換)

バタフライ演算の実装にあたり、 WnkW_n^k をどう計算するか。ここで FFT は精度が悪いらしい ので、競プロでは eiωne^{i\omega_n} でなく mod998244353\bmod 998244353 の世界で直行基底を定義してフーリエ変換を行うようです。

感想

くう〜疲れましたw まだ道半ばです。たぶん実装は遅延セグメント木より簡単そうかな……