學了好久，終於基本弄明白了
推薦兩個博客：
戳我
戳我
再推薦幾本書：
《ACM/ICPC算法基礎訓練教程》
《組合數學》（清華大學出版社）
《高中數學選修》
預備知識
復數方面
找數學老師去
\[i^{2}=-1，i為虛數的單位 \]
坐標系上縱軸就是虛數軸，復數就是這上面的點
三種表示法:
$$一般：a + bi，a為實部，b為虛部$$
$$指數：e^{i\theta}坐標系上的模長$$
$$三角：模長(cos\theta + i sin \theta)$$
運算：
加減法：實部虛部分別相加
乘法：$$(a + bi) * (c + di) = ac + adi + bci + bdi^{2}
= ac-bd+(ad+bc)i$$
歐拉公式
\[e^{ix} = cosx + isinx(就是指數表示和三角表示) \]
\[特別的e^{i\pi} = -1 \]
多項式
\[系數表示法：A(x) = \Sigma _{k=0}^{n - 1} a_kx^k \]
\[點值表示法：對於所有的x_k，求出它們對應的A(x)，設為y_k \]
\[則可以用\{(x_0, y_0), (x_1, y_1), ......, (x_n-1, y_n-1)\} 表示這個多項式 並且是唯一確定的\]
單位復數根
\[n次單位復數根\omega^{n} = 1，n次單位復數根剛好有n個對應e^{\frac{2k\pi i}{n}}，其中k=0到n - 1 \]
三個性質：
消去引理：
$$n, d, k為正整數，則\omega^{dk}{dn}=\omega^{k}{n}$$
$$證明：套e^{\frac{2k\pi i}{n}} 即可$$
折半引理：
$$n為大於零的偶數，則(\omega^{k+\frac{n}{2}}{n})^{2}=\omega^{2k+n}{n}=\omega^{2k}{n}\omega^{n}{n}=(\omega^{k}{n})^{2}$$
求和引理：
大於1的整數n，和不被n整除的非負整數k，有
$$\Sigma^{n-1}{j=0}(\omega^{k}_{n})^{j}=0$$
證明可以用等比數列求和公式得到（很簡單的，手推一遍就好）
Rader排序
其實就是二進制數位翻轉
正題
DFT
對於k=0~n-1，定義：
\[y_k=A(\omega^{k}_{n}) = \Sigma^{n-1}_{j=0} a_j(\omega^{k}_{n})^j \]
\[得到的y稱為a的離散傅里葉變換，記作y=DFT_n(a) (這里的y，a指的是所有的y_k, a_k，即向量y,a) \]
逆DFT
\[就是DFT的逆變換，求出向量a，記為DFT^{-1} \]
假設得到了向量y
\[對於y_k = \Sigma^{n-1}_{i=0}a_i(\omega^{k})^i \]
\[有a_k = \frac{1}{n}\Sigma^{n-1}_{i=0}y_i(\omega^{-k})^i \]
\[證明：a_k=\frac{1}{n}\Sigma^{n-1}_{i=0}y_i(\omega^{-k})^i=\frac{1}{n}\Sigma^{n-1}_{i=0}( \Sigma^{n-1}_{j=0}a_j(\omega^{k})^j)(\omega^{-k})^i=\frac{1}{n} \Sigma^{n-1}_{i=0}a_i(\Sigma^{n-1}_{j=0}(\omega^{j-k})^i) \]
\[可以用等比數列求和出上面的就是a_k（當j\ne k是括號里的為0，當j=k時為1） \]
FFT
上面已經把DFT和逆DFT搞定了，兩個幾乎是一樣的
所以求多項式的積（卷積）可以用DFT轉換成點值表示，就可以O(n)，一一相乘，得到積的多項式的點值表示，最后用逆DFT得到系數表示
復雜度瓶頸在於怎樣快速求解DFT（逆DFT和DFT方法一樣）
FFT就是一個O(nlogn)求解DFT的方法
首先把A(x)分成奇數項和偶數項記作
\[A^{[0]}(x) = a_0 + a_2x + a_4x^2 + ... + a_{n-2}x^{\frac{n}{2} - 1} \]
\[A^{[1]}(x) = a_1 + a_3x + a_5x^2 + ... + a_{n-1}x^{\frac{n}{2} - 1} \]
\[顯然A(x) = A^{[0]}(x^2) + xA^{[1]}(x^2) \]
那么
\[A(\omega^k_n)=A^{[0]}((\omega^k_n)^2) + \omega^k_n A^{[1]}((\omega^k_n)^2)=A^{[0]}(\omega^k_{\frac{n}{2}}) + \omega^k_n A^{[1]}(\omega^k_{\frac{n}{2}}) \]
\[因為\omega^{\frac{n}{2}}_{n}=\omega_{2}=e^{k\pi i}= cos k\pi + i sin k\pi = -1 \]
\[所以A(\omega^{k+\frac{n}{2}}_n)=A^{[0]}(\omega^k_{\frac{n}{2}}) - \omega^k_n A^{[1]}(\omega^k_{\frac{n}{2}}) \]
這稱為蝴蝶操作
於是對每個y值的求解可以通過分組求出，若遞歸變成處理子任務，這樣復雜度就成了O(nlogn)
這樣不停地分組，最后就相當於Rader排序了一番，所以也可以變成非遞歸的
注意每次都要把多項式補成2的冪，便於FFT
遞歸寫可能好理解一些，但不好寫
還有一些東西什么的，其實記一記就好了其實自己說不清
系統的復數complex代碼
# include <bits/stdc++.h>
# define RG register
# define IL inline
# define Fill(a, b) memset(a, b, sizeof(a))
using namespace std;
typedef long long ll;
const int _(3e6 + 10);
const double Pi = acos(-1);
IL ll Read(){
char c = '%'; ll x = 0, z = 1;
for(; c > '9' || c < '0'; c = getchar()) if(c == '-') z = -1;
for(; c >= '0' && c <= '9'; c = getchar()) x = x * 10 + c - '0';
return x * z;
}
int n, m, r[_], l;
complex <double> a[_], b[_];
IL void FFT(complex <double> *P, int opt){
for(RG int i = 0; i < n; ++i) if(i < r[i]) swap(P[i], P[r[i]]); //Rader排序
for(RG int i = 1; i < n; i <<= 1){
complex <double> W(cos(Pi / i), opt * sin(Pi / i)); //旋轉因子
for(RG int p = i << 1, j = 0; j < n; j += p){
complex <double> w(1, 0);
for(RG int k = 0; k < i; ++k, w *= W){
complex <double> X = P[j + k], Y = w * P[j + k + i];
P[j + k] = X + Y; P[j + k + i] = X - Y; //蝴蝶操作
}
}
}
}
int main(RG int argc, RG char *argv[]){
n = Read(); m = Read();
for(RG int i = 0; i <= n; ++i) a[i] = Read();
for(RG int i = 0; i <= m; ++i) b[i] = Read();
m += n;
for(n = 1; n <= m; n <<= 1) ++l;//補成2的冪
for(RG int i = 0; i < n; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));//Rader排序預處理
FFT(a, 1); FFT(b, 1); //DFT
for(RG int i = 0; i < n; ++i) a[i] = a[i] * b[i]; //點值直接相乘
FFT(a, -1); //逆DFT
for(RG int i = 0; i <= m; ++i) printf("%d ", (int)(a[i].real() / n + 0.5));
return 0;
}
或者可以自己定義complex，用復數運算
struct Complex{
double real, image;
IL Complex(){ real = image = 0; }
IL Complex(RG double a, RG double b){ real = a; image = b; }
IL Complex operator +(RG Complex B){ return Complex(real + B.real, image + B.image); }
IL Complex operator -(RG Complex B){ return Complex(real - B.real, image - B.image); }
IL Complex operator *(RG Complex B){ return Complex(real * B.real - image * B.image, real * B.image + image * B.real); }
}
NTT(快速數論變換)
前置技能原根
設\(g\)為\(p\)(質數)的原根
則\(e^{\frac{2\pi i}{n}}\equiv\omega_n\equiv g^{\frac{p-1}{n}}(mod \ p)\)
帶進去就好了
Reverse的那個不會證明
\(UOJ\)的模板
# include <bits/stdc++.h>
# define RG register
# define IL inline
# define Fill(a, b) memset(a, b, sizeof(a))
using namespace std;
typedef long long ll;
const int Zsy(998244353);
const int Phi(998244352);
const int G(3);
const int _(4e5 + 5);
IL ll Input(){
RG ll x = 0, z = 1; RG char c = getchar();
for(; c < '0' || c > '9'; c = getchar()) z = c == '-' ? -1 : 1;
for(; c >= '0' && c <= '9'; c = getchar()) x = (x << 1) + (x << 3) + (c ^ 48);
return x * z;
}
int n, m, N, l, r[_], A[_], B[_];
IL int Pow(RG ll x, RG ll y){
RG ll ret = 1;
for(; y; y >>= 1, x = x * x % Zsy)
if(y & 1) ret = ret * x % Zsy;
return ret;
}
IL void NTT(RG int *P, RG int opt){
for(RG int i = 0; i < N; ++i) if(r[i] < i) swap(P[r[i]], P[i]);
for(RG int i = 1; i < N; i <<= 1){
RG int W = Pow(G, Phi / (i << 1));
if(opt == -1) W = Pow(W, Zsy - 2);
for(RG int j = 0, p = i << 1; j < N; j += p){
RG int w = 1;
for(RG int k = 0; k < i; ++k, w = 1LL * w * W % Zsy){
RG int X = P[k + j], Y = 1LL * w * P[k + j + i] % Zsy;
P[k + j] = (X + Y) % Zsy, P[k + j + i] = (X - Y + Zsy) % Zsy;
}
}
}
}
int main(RG int argc, RG char* argv[]){
n = Input(), m = Input();
for(RG int i = 0; i <= n; ++i) A[i] = Input();
for(RG int i = 0; i <= m; ++i) B[i] = Input();
for(n += m, N = 1; N <= n; N <<= 1) ++l;
for(RG int i = 0; i < N; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
NTT(A, 1); NTT(B, 1);
for(RG int i = 0; i < N; ++i) A[i] = 1LL * A[i] * B[i] % Zsy;
NTT(A, -1);
RG int inv = Pow(N, Zsy - 2);
for(RG int i = 0; i <= n; ++i) printf("%lld ", 1LL * A[i] * inv % Zsy);
return 0;
}