什么是FFT#
既然打開了這篇博客，大家肯定都已經對FFT(Fast Fourier Transformation)有一點點了解了吧
FFT即為快速傅里葉變換，可以快速求卷積（當然不止這一些應用，但是我不會）
系數表示法與點值表示法#
我們通常表示一個\(n-1\)次多項式是利用系數表示法like this:\(f(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}\)
點值表示法即為將多項式用坐標系上的若干個點表示
我們對這個多項式代入不同的值{\(x_1,x_2,...,x_n\)}
我們就可以得到\(n\)個點\((x_1,f(x_1)),(x_2,f(x_2)),...,(x_n,f(x_n))\)
實際上只要保證代入的\(n\)個數互不相同，那么這\(n\)個點就對應了唯一的\(f(x)\)（類似\(n\)元一次方程組？）
點值表示法它的優勢在哪呢？
我們如果想要求出兩個多項式相乘，系數表示法就很麻煩，而點值表示法卻只需要將相同的\(x\)對應的點值相乘就行了
DFT與IDFT#
DFT##
我們若將一個多項式強行轉為點值表示法則時間復雜度為\(O(n^2)\)
自然有人表示：太慢啦，就不能快點嗎！
於是我們開始思考，如果我們對於\((x_1,f(x_1)),(x_2,f(x_2)),...,(x_n,f(x_n))\)取一系列特殊的有關系的值，是不是能加速呢？
當然可以，要不然就不會有這篇博客了
在DFT中，我們可以將復平面中的單位園n等分，然后每一等分的頂點\(\omega _n^i\)作為\(x_i\)，\(\omega _n^i\)就是單位根
（如圖為\(8\)等分）
其中\(\omega _n^i=\cos (\frac{{2\pi i}}{n}) + i\sin (\frac{{2\pi i}}{n})\)
當然DFT只是讓\((x_1,f(x_1)),(x_2,f(x_2)),...,(x_n,f(x_n))\)變得特殊了一些，並沒有優化復雜度的效果（要不然哪來的FFT）
IDFT##
IDFT(Inverse Discrete Fourier Transform)顧名思義就是DFT的逆變換（當然我不會，這個東西，可以會，但沒必要）
關於單位根#
首先，顯然\(\omega _{2n}^{2i}=\omega _n^i\)（結合圖像可得）\((n=2^k)\)
\(\omega _n^{i+\frac{n}{2}}=-\omega _n^i\)為什么呢，因為這兩個單位根在坐標中終點關於原點對稱
最后一個沒有多大用的東西，但是后面會用到\(\omega _n^0=\omega _n^n=1\)，這個顯然吧，因為這個頂點是重合的，且虛部為0
FFT#
遞歸解法##
FFT(Fast Fourier Transformation)快速傅里葉變換，有了上面這些東西，我們就可以考慮對原多項式進行變形
\(f(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}\)
\(=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+x(a_1+a_3x^2+...+a_{n-1}x^{n-2})\)
是不是發現這兩個東西很像啊
我們設\(f_1(x)=(a_0+a_2x+...+a_{n-2}x^{\frac{n}{2}-1}),f_2(x)=(a_1+a_3x+a_{n-1}x^{\frac{n}{2}-1})\)
所以\(f(x)=f_1(x^2)+xf_2(x^2)\)
我們發現，\(f_1(x),f_2(x)\)和\(f(x)\)的形式一樣誒！
也就是說我們現在這個FFT已經可以用遞歸解決了！
然而，遞歸自帶sb大常數，時間上並不能通過
於是，厲害的人們又發現了新的操作，蝴蝶變換
蝴蝶變換##
首先我們令\(i<\frac{n}{2}\)
\(f(\omega _n^i)=f_1({(\omega _n^i)}^2)+\omega _n^if_2({(\omega _n^i)}^2)\)
\(=f_1(\omega _n^{2i})+\omega _n^if_2(\omega _n^{2i})\)
\(=f_1(\omega _{\frac{n}{2}}^i)+\omega _n^if_2(\omega _{\frac{n}{2}}^i)\)
\(f(\omega _n^{i+\frac{n}{2}})=f_1({(\omega _n^i)}^2\omega _n^n)+\omega _n^{i+\frac{n}{2}}f_2({(\omega _n^i)}^2\omega _n^n)\)
\(=f_1(\omega _n^{2i})-\omega _n^if_2(\omega _n^{2i})\)
\(=f_1(\omega _{\frac{n}{2}}^i)-\omega _n^if_2(\omega _{\frac{n}{2}}^i)\)
是不是很神奇，是不是！
就差了一個符號，是不是很蝴蝶啊(親愛的，你慢慢飛，小心...)
有了蝴蝶變換，我們發現我們只要知道\(f_1(\omega _{\frac{n}{2}}^i),f_2(\omega _{\frac{n}{2}}^i)\)就可以求出\(n\)等分下兩個不同單位根帶入后的點值了，是不是很厲害！
rev數組##
但是我們又發現了一個新的問題，每次我們合並的時候，都是把要求的序列分為奇數位和偶數位做
可是我們原數組並不符合這個條件，直接做的話肯定死
我們可以研究一下每個數字和它的二進制表示
以8個數為例：
但是我們合並的順序應該長成這樣：
（讀者可以自己模擬一下，段長每次為\(2^k\)，然后合並）
我們發現，其實如果把每個數的二進制位反過來一下，就變成順序的了！
把二進制反過來其實也不復雜，就一行代碼
rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1))
其實這段代碼的意義大概是把除了最低位之外的位用i>>1翻轉，然后再將最低位移到最高位
也就是說，做FFT的時候只要將系數重新排列，然后每段進行合並更新就好啦（這段看一下代碼就會很清楚）
IFFT#
IFFT(Inverse Discrete Fourier Transform)快速傅立葉逆變換
實際上我們的傅立葉變換可以寫成矩陣的形式：
\(\left[ {\begin{array}{*{20}{c}} 1&1&1&{...}&1\\ 1&{\omega _n^1}&{\omega _n^2}&{...}&{\omega _n^{n - 1}}\\ 1&{\omega _n^2}&{\omega _n^4}&{...}&{\omega _n^{2n - 2}}\\ {...}&{...}&{...}&{...}&{...}\\ 1&{\omega _n^{n - 1}}&{\omega _n^{2n - 2}}&{...}&{\omega _n^{(n - 1) \times (n - 1)}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{a_0}}\\ {{a_1}}\\ {{a_2}}\\ {...}\\ {{a_{n - 1}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {f(\omega _n^0)}\\ {f(\omega _n^1)}\\ {f(\omega _n^2)}\\ {...}\\ {f(\omega _n^{n - 1})} \end{array}} \right]\)
當然，看到這個形式，我們第一個想法就是對
\(\left[ {\begin{array}{*{20}{c}} 1&1&1&{...}&1\\ 1&{\omega _n^1}&{\omega _n^2}&{...}&{\omega _n^{n - 1}}\\ 1&{\omega _n^2}&{\omega _n^4}&{...}&{\omega _n^{2n - 2}}\\ {...}&{...}&{...}&{...}&{...}\\ 1&{\omega _n^{n - 1}}&{\omega _n^{2n - 2}}&{...}&{\omega _n^{(n - 1) \times (n - 1)}} \end{array}} \right]\)
這個矩陣求逆，但是這樣的話，復雜度就吃不消了，我們好不容易把\(O(n^2)\)的DFT降為\(O(nlogn)\)的FFT，結果逆不回去，這不是非常尷尬嗎
但是我們發現，我們對這個矩陣求一下逆，可以得到這個形式是固定的
\(\left[ {\begin{array}{*{20}{c}} {\frac{1}{n}}&{\frac{1}{n}}&{\frac{1}{n}}&{...}&{\frac{1}{n}}\\ {\frac{1}{n}}&{\frac{1}{n}\omega _n^{ - 1}}&{\frac{1}{n}\omega _n^{ - 2}}&{...}&{\frac{1}{n}\omega _n^{1 - n}}\\ {\frac{1}{n}}&{\frac{1}{n}\omega _n^{ - 2}}&{\frac{1}{n}\omega _n^{ - 4}}&{...}&{\frac{1}{n}\omega _n^{2 - 2n}}\\ {...}&{...}&{...}&{...}&{...}\\ {\frac{1}{n}}&{\frac{1}{n}\omega _n^{1 - n}}&{\frac{1}{n}\omega _n^{2 - 2n}}&{...}&{\frac{1}{n}\omega _n^{ - (n - 1) \times (n - 1)}} \end{array}} \right]\)
然后我們可以把\(\frac{1}{n}\)提出來，那么我們就可以先利用快速傅立葉變換，做出系數的\(n\)倍，然后再除掉
坑終於填完啦
FFT_Code#
這是一個高精乘的模板
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
char s[60050],s1[60050];
struct r{
double real,image;
}a[131090],b[131090];
double pi=acos(-1);
r operator +(r a,r b){return (r){a.real+b.real,b.image+a.image};}
r operator -(r a,r b){return (r){a.real-b.real,a.image-b.image};}
r operator *(r a,r b){return (r){a.real*b.real-b.image*a.image,a.real*b.image+a.image*b.real};}
int ans[131090],rev[131090],bit,n,l;
void fft(r *a,int n,int op){
for (int i=0;i<n;i++)if (i<rev[i])std::swap(a[rev[i]],a[i]);
for (int i=1;i<n;i<<=1){
r wn=(r){cos(pi/i),op*sin(pi/i)};//將單位圓分成i*2個部分
for (int j=0;j<n;j+=i<<1){
r wnk=(r){1,0};
for (int k=j;k<i+j;k++,wnk=wnk*wn){
r x=a[k],y=wnk*a[i+k];
a[k]=x+y;
a[k+i]=x-y;//這里很蝴蝶
}
}
}
if (op==-1)for (int i=0;i<n;i++)a[i].real/=n;
}
int main(){
scanf("%d",&l);
scanf("%s",s);
scanf("%s",s1);
int n=2;
for (bit=1;(1<<bit)<(l<<1)-1;bit++)n<<=1;
for (int i=0;i<l;i++)a[i]=(r){(double)(s[l-i-1]-'0'),0};
for (int i=0;i<l;i++)b[i]=(r){(double)(s1[l-i-1]-'0'),0};
for (int i=0;i<n;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));//處理rev數組
fft(a,n,1);
fft(b,n,1);
for (int i=0;i<n;i++)a[i]=a[i]*b[i];
fft(a,n,-1); //IFFT
for (int i=0;i<n;i++){
ans[i]+=(a[i].real+0.5);
ans[i+1]=ans[i]/10;
ans[i]%=10;
}
int i;
for (i=n;i>=0&&!ans[i];i--);
if (i<0)putchar('0');
else for (;i>=0;i--)printf("%d",ans[i]);
}
什么是NTT#
NTT(Number Theoretic Transform)快速數論變換（不要問我快速是哪來的，我不知道）
可以將類似FFT的做法用整數來做
與FFT區別##
1、可以取模（必須取模）
2、比FFT快這是一個遙遠的傳說，究竟哪個快還是和數據有關
原根#
對於一個正整數m，存在一個整數g，滿足\(g^{\varphi (m)}\equiv 1(mod\) \(m)\)，則g為m的一個原根
一個數存在原根的充要條件是該數可以表示為\(2,4,p^k,2*p^k\)(\(p\)為奇質數,\(k \geq 1\))
NTT#
有了原根我們就可以考慮用原根來代替單位根，因為原根也有類似的性質（一般此類題目模數為質數，所以以下用\(p-1\)代替\(\varphi(p)\)）
那么我們考慮用\(G_n^i\)=\(g^{i(p-1)/n}\)來代替\(\omega _n^i\)
為什么這樣代替呢，因為這樣代替就滿足：
\(G_n^0=G_n^n=1\)，這個顯然吧？
\(G_n^i=-G_n^{i+n/2}\)，這個通過原根的定義稍微推一下就能得到
\(G _{2n}^{2i}=G _n^i\)，這個應該問題也不大吧？
然后我們就發現NTT它鍋了，哪里鍋了呢？
性質當然是沒有問題的，鍋就鍋在\(n \nmid (p-1)\)的時候，我們沒辦法做浮點數的快速冪
那我們該怎么辦，難道拋棄原根去尋找新方法？
Of course not.我們就規定NTT的模數只能取\(k*2^m+1\)並且\(n \leq 2^m\)就行啦（最常用的可以NTT的模數就是998244353，如果出現新的大家自己檢驗一下啦）
（順便送給各位，998244353的一個原根是3是不是非常良心）
NTT_Code#
這還是一個高精乘的模板
#include<cstdio>
#include<algorithm>
int g=3,mo=998244353,a[131090],b[131090],rev[131090],l,bit;
char s[60050],s1[60050];
int ksm(int p,int k){
int ret=1;
while (k){
if ((k&1))ret=(ret*1ll*p)%mo;
p=(p*1ll*p)%mo;k>>=1;
}
return ret;
}
void ntt(int *a,int n,int f){
for (int i=0;i<n;i++)if (i<rev[i])std::swap(a[i],a[rev[i]]);
for (int i=1;i<n;i<<=1){
int gn=ksm(g,(mo-1)/(i<<1));if (f==-1)gn=ksm(gn,mo-2);
for(int j=0;j<n;j+=i<<1){
int gqn=1;
for (int k=j;k<j+i;k++){
int x=a[k],y=gqn*1ll*a[k+i]%mo;
a[k]=(x+y)%mo;
a[k+i]=(x-y+mo)%mo;
gqn=(gqn*1ll*gn)%mo;
}
}
}
}
int main(){
scanf("%d",&l);
scanf("%s",s);
scanf("%s",s1);
int n=2;
for (bit=1;(1<<bit)<(l<<1)-1;bit++)n<<=1;
for (int i=0;i<l;i++)a[i]=s[l-i-1]-'0';
for (int i=0;i<l;i++)b[i]=s1[l-i-1]-'0';
for (int i=0;i<n;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
ntt(a,n,1);
ntt(b,n,1);
for (int i=0;i<n;i++)a[i]=(a[i]*1ll*b[i])%mo;
ntt(a,n,-1);
int ny=ksm(n,mo-2);
for (int i=0;i<n;i++)a[i]=(a[i]*1ll*ny)%mo;
for (int i=0;i<n;i++)a[i+1]+=a[i]/10,a[i]%=10;
int i;
for (i=n;i>=0&&!a[i];i--);
if (i==-1)printf("0");
else for (;i>=0;i--)printf("%d",a[i]);
}
一點小應用#
多項式求逆##
多項式求逆是個啥呢，就是我們有多項式\(f(x)\)我們希望求出多項式\(g(x)\)滿足
\(f(x)g(x)\equiv 1(mod\) \(x^n)\)
假設我們現在已經得到了\(p(x)\)滿足\(f(x)p(x)\equiv 1(mod\) \(x^{\left \lceil \frac{n}{2} \right \rceil})\)
因為\(g(x)\)滿足\(f(x)g(x)\equiv 1(mod\) \(x^n)\)，則必然滿足\(f(x)g(x)\equiv 1(mod\) \(x^{\left \lceil \frac{n}{2} \right \rceil})\)
所以\(f(x)(g(x)-p(x))\equiv 0(mod\) \(x^{\left \lceil \frac{n}{2} \right \rceil})\)
因為\(f(x)\not\equiv 0(mod\) \(x^{\left \lceil \frac{n}{2} \right \rceil})\)（否則不存在\(g(x)\)，可以證明這樣的多項式不存在逆多項式）
所以\(g(x)-p(x)\equiv 0(mod\) \(x^{\left \lceil \frac{n}{2} \right \rceil})\)
將其平方一下得到：
\(g^2(x)-2g(x)p(x)+p^2(x)\equiv 0(mod\) \(x^n)\)（為什么后面mod的東西也平方了？因為\(x\)次多項式和\(y\)次多項式相乘后得到的是\(x+y\)次的多項式）
\(g^2(x)\equiv 2g(x)p(x)-p^2(x)(mod\) \(x^n)\)
\(g(x)\equiv 2p(x)-\frac{p^2(x)}{g(x)}(mod\) \(x^n)\)
當然\(f(x)\)與\(g(x)\)是互逆的，所以\(g(x)\equiv 2p(x)-f(x)p^2(x)(mod\) \(x^n)\)
所以\(g(x)\equiv p(x)(2-f(x)p(x))(mod\) \(x^n)\)
到這里，我們發現，我們只要能夠求出\(p(x)\)就可以求出\(g(x)\)，那遞歸求解就好啦
void inv(int *b,int n){
if (n==1){a[0]=ksm(b[0],mo-2);return;}
inv(b,(n+1)>>1);
int bit,len=2;
for (bit=1;(1<<bit)<(n<<1);bit++)len<<=1;
for (int i=0;i<n;i++)d[i]=b[i];
for (int i=n;i<len;i++)d[i]=0;
for (int i=0;i<len;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
ntt(a,len,1);ntt(d,len,1);
for (int i=0;i<len;i++)a[i]=(2ll-d[i]*1ll*a[i]%mo+mo)%mo*1ll*a[i]%mo;
ntt(a,len,-1);
for (int i=n;i<len;i++)a[i]=0;
}
多項式取ln##
相信大家已經熟練掌握了FFT（可能我還沒？）我們就來學習一下多項式取ln的藝術吧
對於一個多項式\(f(x)\)，我們想求出\(ln(f(x))\)
令\(g(x)=ln(f(x))\)，兩邊同時求導得：
\(g'(x)=\frac{f'(x)}{f(x)}\)
那么\(g(x)=\int \frac{f'(x)}{f(x)}\)
那我們就只要做一下多項式求導，多項式求逆，多項式積分就好啦