【luogu P3803】【模板】多項式乘法（FFT）
【模板】多項式乘法（FFT）
題目鏈接：luogu P3803
題目大意
給你兩個多項式，要你求這兩個多項式乘起來得到的多項式。（卷積）
思路
係數表示法
就是我們一般來表示一個多項式的方法：
\(A(x)=a_1x^k+a_2x^{k-1}+…+a_k\)
或者可以這樣表示：
\(A(x)=\sum\limits_{i=1}^{k}a_i\times x_i\)
那你很容易看到，用來做這道題用係數表示法來做是 \(O(n^2)\) 的。
點值表示法
假設我們已經知道了這個多項式，那我們把 \(n\) 個不同的 \(x\) 帶入，會得到 \(n\) 個取值 \(y\)。
由於 \(n\) 個點可以確定一個 \(n\) 次多項式，那我們就可以根據這 \(n\) 個點確定這個多項式。
那你如果選了 \(n\) 個點 \((x_1,y_1),…,(x_n,y_n)\)，那就有：
\(y_i=\sum\limits_{j=0}^{n-1}a_j\times x_i^j\)
但用它計算還是 \(O(n^2)\)，你選點 \(O(n)\)，對於每個點計算也是 \(O(n)\)。
考慮優化
至於係數表示法，它的係數固定，你一改其它都要改，似乎很難弄。
但是第二種點值法似乎沒有特別大的限制讓它難以優化。
但是怎麼優化呢？這就要看點前置知識才可以繼續了。
前置知識
向量
普通的量只有大小，但是向量就是有方向又有大小的量。
在幾何裡面它其實就是一個箭頭。
圓的弧度制
等於半徑長的圓弧所對的圓心角叫做1弧度的角，叫做弧度，用 rad 表示。
用它做單位來量角的制度就是弧度制。
\(1^{\circ}=\frac{\pi}{180}rad\)
\(180^{\circ}=\pi rad\)
複數
設 \(a,b\) 為實數，\(i^2=-1\)，那形如 \(a+bi\) 的數就叫做複數。
\(i\) 就是其中的複數單位。
複數其實可以在一個二維平面表示出來，\(x\) 軸表示實數 \(a\)，\(y\) 軸（當然在原點的時候不是虛數）表示虛數 \(bi\)。
然後從 \((0,0)\) 到 \((a,b)\) 的向量就表示了複數 \(a+bi\)。
模長：就是那個點到原點的距離，勾股就可以得到。
幅角：從 x 正半軸以逆時針為正方向到一直向量的轉角。
有關複數的運演算法則
那由於複數在平面上是向量，那我們其實可以用向量的加減法則來弄。（就是用那個平行四邊形定則）
然後乘法就稍微複雜一點。
幾何的意義就是，複數相乘，模長也相乘，然後幅角就會相加。
其實我們這裡要的是代數意義，我們化一下式子：
\((a+bi)*(c+di)\)
\(=ac+a*di+b*ci+bi*di\)
\(=ac+i(ad+bc)-bd\)（\(i^2=-1\)）
\(=(ac-bd)+(ad+bc)i\)
單位根
在複平面上，以原點為圓心，\(1\) 為半徑作圓，所得的圓叫單位圓。
以圓點為起點，圓的 \(n\) 等分點為終點，做 \(n\) 個向量，設幅角為正且最小的向量對應的複數為 \(\omega_n\)，稱為 \(n\)次單位根。
那剩下的那 \(n-1\) 個就是 \(\omega_n^2,\omega_n^3,…,\omega_n^n\)。
然後要注意的就是 \(\omega_n^0\) 和 \(\omega_n^n\) 都是 \(1\)，對於的是 \(x\) 正半軸的向量。
那怎麼計算其它的呢？
我們可以用歐拉公式來求：\(\omega_n^k=\cos k\dfrac{2\pi}{n}+i\sin k\dfrac{2\pi}{n}\)
（大概就是把模長那條線往 \(x\) 軸做垂線得到直角三角形，然後用三角函數得到兩條直角邊的長度，從而得到坐標，也就是複數）
當然顯而易見，單位根幅角就是 \(\frac{1}{n}\)
還有就是如果 \(z^n=1\)，那 \(z\) 就是 \(n\) 次單位根。
單位根的一些性質
\(\omega_{2n}^{2k}=\omega_n^k\)
證明：
\(\omega_{2n}^{2k}=\cos 2k\dfrac{2\pi}{2n}+i\sin 2k\dfrac{2\pi}{2n}\)
\(=\cos k\dfrac{2\pi}{n}+i\sin k\dfrac{2\pi}{n}=\omega_n^k\)
\(\omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^{k}\)
證明：
\(\omega_{n}^{\frac{n}{2}}=\cos \dfrac{n}{2}*\dfrac{2\pi}{n}+i\sin \dfrac{n}{2}*\dfrac{2\pi}{n}\)
\(=\cos \pi+i\sin\pi=-1\)
（\(\cos \pi=\cos 180^\circ=-1,\sin\pi=\sin 180^\circ=0\)）
那就有 \(\omega_{n}^{k+\frac{n}{2}}=\omega_{n}^{k}*\omega_{n}^{\frac{n}{2}}=\omega_{n}^{k}*(-1)=-\omega_{n}^{k}\)
OK，現在前置知識結束，開始正文。
快速傅里葉變換（FFT）
我們設 \(A(x)\) 的係數為 \((a_0,a_1,…,a_{n-1})\)（\(n\) 次多項式）
那就有 \(A(x)=a_0+a_1x+a_2{x^2}+a_3*{x^3}+a_4*{x^4}+a_5*{x^5}+ \dots+a_{n-2}*x^{n-2}+a_{n-1}*x^{n-1}\)
我們可以把它奇偶分開，\(A(x)=(a_0+a_2x^2+…+a_{n-2}x^{n-2})+(a_1+a_3x^3+…+a_{n-1}x^{n-1})\)
那我們可以設 \(A_1(x)=a_0+a_2x+…+a_{n-2}x^{\frac{n}{2}-1},A_2(x)=a_1+a_3x+…+a_{n-1}x^{\frac{n}{2}-1}\)
那我們容易看出 \(A(x)=A_1(x^2)+xA_2(x^2)\)。
把 \(\omega_n^k\) 帶入，\(A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})\)
把 \(\omega_n^{k+\frac{n}{2}}\) 帶入，\(A(\omega_n^k)=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k+n})\)，然後把這個式子化一下：
\(=A_1(\omega_n^{2k}*\omega_n^n)-\omega_n^{k}A_2(\omega_n^{2k}*\omega_n^n)\)
\(=A_1(\omega_n^{2k})-\omega_n^{k}A_2(\omega_n^{2k})\)
你會發現這兩個東西不同的只有一個常數項。
那我們可以枚舉得到第一個答案，然後通過第一個答案得到第二個。
那你會發現你只用枚舉前面的一半，那問題規模就減半。
那分出的兩個問題可以繼續分半解決，那我們可以用遞歸來搞。
那這個其實就是類似分治的感覺，是 \(O(nlogn)\) 的。
快速傅里葉逆變換（IFFT）
前面我們在弄的時候都是用點值表示法。
但一般我們（題目）給的多半都是係數表示法。
從係數表示法得到點值表示法我們已經學會了，但是怎麼轉回來呢？
這個時候就要用到 IFFT 了。
\((y_0,y_1,…,y_{n-1})\) 為 \((a_0,a_1,…,a_{n-1})\) 的傅里葉變換（就是點值表示）
那另外有一個 \((c_0,c_1,…,c_{n-1})\) 滿足：
\(c_k=\sum\limits_{i=0}^{n-1}y_i(\omega_n^{-k})^i\)
（其實就是把 \((y_0,y_1,…,y_{n-1})\) 當做多項式，求它在 \(\omega_n^0,\omega_n^{-1},…,\omega_n^{-(n-1)}\) 的點值表示）
然後來推推公式：
\(c_k=\sum\limits_{i=0}^{n-1}y_i(\omega_n^{-k})^i\)
\(=\sum\limits_{i=0}^{n-1}(\sum\limits_{j=0}^{n-1}a_j(\omega_n^i)^j)(\omega_n^{-k})^i\)
\(=\sum\limits_{i=0}^{n-1}(\sum\limits_{j=0}^{n-1}a_j(\omega_n^j)^i)(\omega_n^{-k})^i\)
\(=\sum\limits_{i=0}^{n-1}(\sum\limits_{j=0}^{n-1}a_j(\omega_n^j)^i(\omega_n^{-k})^i)\)
\(=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}a_j(\omega_n^{j-k})^i\)
\(=\sum\limits_{j=0}^{n-1}a_j\sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i\)
我們設 \(S(x)=\sum\limits_{i=0}^{n-1}x^i\)。
帶入 \(\omega_n^k\)：\(S(\omega_n^k)=\sum\limits_{i=0}^{n-1}(\omega_n^k)^i\)
然後分類討論：
當 \(k=0\) 時，\(\omega_n^k=\omega_n^0=1\)。
那式子就變成了 \(S(\omega_n^k)=\sum\limits_{i=0}^{n-1}1^i=n\)
當 \(k\neq0\) 時，我們考慮怎麼求。
觀察到沒項都乘了一個值，我們考慮用這麼一種方法：
乘 \(\omega_n^k\)：\(\omega_n^kS(\omega_n^k)=\sum\limits_{i=1}^{n}(\omega_n^k)^i\)
然後兩個相減：\(\omega_n^kS(\omega_n^k)-S(\omega_n^k)=(\omega_n^k)^n-(w_n^k)^0\)
\((\omega_n^k-1)S(\omega_n^k)=(\omega_n^k)^n-1\)
\((\omega_n^k-1)S(\omega_n^k)=(\omega_n^n)^k-1\)
\((\omega_n^k-1)S(\omega_n^k)=1^k-1\)
\(S(\omega_n^k)=\dfrac{1-1}{\omega_n^k-1}=0\)
OK，分析完了 \(S(x)\)，我們回到之前卡克的地方：
\(c_k=\sum\limits_{j=0}^{n-1}a_j\sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i\)
\(=\sum\limits_{j=0}^{n-1}a_jS(\omega_n^{j-k})\)
那隻會有一個 \(j=k\)，使得 \(j-k=0\)，貢獻是 \(n\)，其它的時候，貢獻都是 \(0\)。
那就是 \(c_k=a_kn\)
那你就會發現，\(a_k=\dfrac{c_k}{n}\)
那我們就可以通過求出 \((c_0,c_1,…,c_{n-1})\)，然後再根據這個一弄就可以得到 \(a_k\) 了。
（而且你會發現求 \((c_0,c_1,…,c_{n-1})\) 和求 FFT 的很像，只是由 \(\omega_n^i\) 變成了 \(\omega_n^{-i}\)，那我們只要再搞一個 \(op\)，記錄它單位的改變即可以了）
演算法總結
它的主要思想就是把係數表示法轉成點值表示法，然後用點值表示法的分支方法來快速得到答案，然後再將答案轉回係數表示法。
實現
妹想到把，還有。
別的地方都沒有問題，我們來講講遞歸的做法。
遞歸
（我一定不會告訴你這個是過不了的）
遞歸其實就是根據我們前面的奇偶分開做法，把一個序列弄成兩個，然後遞歸處理之後合併。
遞歸到只剩一個常數項的時候，就可以直接返回。
然後我們先來一個小小的優化。
它有一個很牛逼的名字——蝴蝶效應。
（其實就是記憶化）
因為向量的乘法耗比較多時間，我們可以先把要乘的乘出來放一個變數里，然後要用的時候直接就是這個變數。
（這樣如果這個乘法值要用多次的話就可以節省時間，原本要乘很多次，現在只用乘一次）
然而就算你加了這個優化，也是過不了的。
因為你遞歸，加上常數比較大，在 \(10^6\) 就爆了。
那我們就要找一個更優的方法來弄，那就是迭代實現。
迭代實現
我們考慮觀察一下原序列和翻轉之後的序列。
我們觀察一下它是怎麼變的：
\(0,1,2,3,4,5,6,7\)
\(0,2,4,6,1,3,5,7\)
\(0,4,2,6,1,5,3,7\)
似乎沒有什麼想法，轉成二進位：
\(000,100,010,110,001,101,011,111\)
再把一開始的也轉成二進位：
\(000,001,010,011,100,101,110,111\)
你會發現它就是把它們的二進位翻轉了。
那你可以用一個 \(O(n)\) 的方法得到對應關係——DP。
\(f_i=(f_{i/2}/2)|((i\&1)^{n-1})\)
大概就是把之前翻轉好的往左邊移一個，流出位置給原序列中右邊新出現的一位放。
接著就是怎麼通過這個翻轉對應得到我們迭代的序列。
現有一個顯然的東西，就是 \(f_{f_i}=i\)
那也就是說，它們之間是兩兩對應的。
那就需要把每對都只翻轉一次，就可以了。
那簡單，我們只要找一個在一對裡面只會發生一次的事，就比如 \(f_i>i\)。（當然你小於也可以，只要讓一對只有一個會發生就可以了）
然後我們來講講迭代要怎麼搞。
其實就想到與從下段直接網上合併。（因為你已經確定了最下的序列）
就枚舉合併的大小，然後枚舉區間，然後就合併。
小聲嗶嗶
終於寫完了，好累啊。
awa
程式碼
#include<cstdio>
#include<cmath>
#include<iostream>
using namespace std;
struct complex {
double x, y;
complex (double xx = 0, double yy = 0) {
x = xx;
y = yy;
}
}a[5000005], b[5000005];
int n, m, limit, l_size;
int an[5000005];
double Pi = acos(-1.0);
complex operator +(complex x, complex y) {//定義向量的加減乘
return complex(x.x + y.x, x.y + y.y);
}
complex operator -(complex x, complex y) {
return complex(x.x - y.x, x.y - y.y);
}
complex operator *(complex x, complex y) {
return complex(x.x * y.x - x.y * y.y, x.x * y.y + x.y * y.x);
}
//op=1則係數變點值，為-1則點值邊係數
//至於為啥看看求的兩個公式就知道了
void FFT(complex *now, int op) {
for (int i = 0; i < limit; i++)
if (i < an[i]) swap(now[i], now[an[i]]);//求出要迭代的序列
for (int mid = 1; mid < limit; mid <<= 1) {//枚舉合併序列的大小
complex Wn(cos(Pi / mid), op * sin(Pi / mid));//單位根
for (int R = mid << 1, j = 0; j < limit; j += R) {//R是區間右端點，j表示左端點
complex w(1, 0);//冪
for (int k = 0; k < mid; k++, w = w * Wn) {//枚舉區間的每個位置
complex x = now[j + k], y = w * now[j + mid + k];//先求出來減少時間
now[j + k] = x + y;//得到一個求出另外一個的
now[j + mid + k] = x - y;
}
}
}
}
int main() {
scanf("%d %d", &n, &m);
for (int i = 0; i <= n; i++) scanf("%lf", &a[i].x);
for (int i = 0; i <= m; i++) scanf("%lf", &b[i].x);
limit = 1;
while (limit <= n + m) {
limit <<= 1;
l_size++;
}
for (int i = 0; i < limit; i++)
an[i] = (an[i >> 1] >> 1) | ((i & 1) << (l_size - 1));
FFT(a, 1);
FFT(b, 1);
for (int i = 0; i <= limit; i++)
a[i] = a[i] * b[i];
//將點值表示法轉換為係數表示法再輸出
FFT(a, -1);
for (int i = 0; i <= n + m; i++)
printf("%d ", (int)(a[i].x / limit + 0.5));
return 0;
}
