【快速傅里葉變換】FFT
參考：從多項式乘法到快速傅里葉變換 by miskcoo
FFT 學習筆記 by Menci
（一）多項式的表示法
系數表示法：f(x)=a[n-1]*x^(n-1)+...+a[0]，稱為n-1次多項式。
點值表示法：一個n-1次多項式在復數域中有n個根，即n個(x,y)可以唯一確定一個n-1次多項式。
對於一個多項式，從其系數表示法到其點值表示法的變換稱為離散傅里葉變換（DFT），反之稱為傅里葉逆變換（IDFT）。
朴素的離散傅里葉變換，枚舉實現的復雜度為O(n^2)。
快速傅里葉變換是指以O(n log n)的復雜度實現IDF和IDFT的算法，常用Cooley-Tukey算法。
（二）復數
復數是形如a+bi的數，當b=0時為實數。
定義一個平面為復平面，那么平面內的每個點(a,b)唯一對應一個復數a+bi，i可以理解為y軸上的單位長度，正如1是x軸上的單位長度。
i的本質是在數軸上定義旋轉變換，i是1逆時針旋轉90°，那么i^2=-1。
復數相加，遵循平行四邊形定則。
復數相乘，模長相乘，幅角相加。
（三）單位根
以圓點為起點，以復平面單位圓的n等分點為終點，作n個向量，設所得幅角為正且最小的向量對應的復數為ω(1,n)，即n次單位根。（括號左為上標，右為下標）。
圖片來源：其中B點是單位根ω(1,n)，逆時針依次為ω(2,n)，ω(3,n)...，ω(n,n)=ω(n,0)=1。
計算公式：ω(k,n)=cos ( 2kπ/n ) + i*sin ( 2kπ/n )
單位根的性質：
(1)消去：ω(2n,2k)=ω(n,k)
(2)折半：ω(n,k+n/2)=-ω(n,k)
將ω(n,0)~ω(n,n-1)這n個單位根作為代表n-1次多項式的n個點的橫坐標，可以得到很好的性質。
（四）快速傅里葉變換（FFT解決DFA）
這部分因為不會操作數學公式，直接粘貼Menci博客QAQ。
將n-1次多項式A(x)的系數奇偶分成兩個多項式A1(x)和A2(x)，則A(x)=A1(x^2)+x*A2(x^2)。
對於k<n/2，有A(ω(n,k))=A1(ω(n/2,k)) + ω(n,k)*A2(ω(n/2,k))
同時，有A(ω(n,k+n/2))=A1(ω(n/2,k)) - ω(n,k)*A2(ω(n/2,k))
對於一個k次多項式，通過奇偶分項得到兩個k/2次多項式，分別計算后再調用其值解決k次多項式，即分治解決。
（五）傅里葉逆變換（IDFA）
對於n-1次多項式，其n-1維系數向量{a0,a1...an-1}通過DFA得到點值向量{b0,b1...bn-1}，反之操作稱為IDFA。
將點值向量作為系數，以單位根的倒數進行FFT，得到的每個數除以n，就是IDFA的結果。
（六）迭代實現FFT
對於多項式A(x)，已知系數向量a[]，求橫坐標為ω(n,0)~ω(n,n-1)的點值向量b[]。
將多項式奇偶分項后，對於k<n/2，有A(ω(n,k))=A1(ω(n/2,k)) + ω(n,k)*A2(ω(n/2,k))，同時有A(ω(n,k+n/2))=A1(ω(n/2,k)) - ω(n,k)*A2(ω(n/2,k))，分治邊界是a[i]*ω(0,0)即a[i]。
邊界元素：FFT遞歸邊界的數組排布恰好是原數組每個位置二進制反轉后的數字，例如：
原：00 01 10 11
終：00 10 01 11
蝴蝶操作：為了在合並時不引入新數組，進行一下操作。ω(l,k)=ω(n,n/l*k)，預處理以n為底的ω[]，IDFT時預處理倒數。
t=ω(n/l*k)*a[l/2+k]
a[l/2+k]=a[k]-t
a[k]+=t
（七）多項式乘法
多項式的點值表示法易於進行乘法，因為對於fc(x)=fa(x)*fb(x)，每個點x在多項式A,B中的點值相乘即可得到在多項式C中的點值。
將n-1次多項式A和m-1次多項式B均視為n+m-2次多項式（高位補0），進行DFT后相乘再通過IDFT即可得到多項式C。
#include<cstdio> #include<cstring> #include<cctype> #include<cmath> #include<queue> #include<stack> #include<set> #include<vector> #include<algorithm> #define ll long long #define lowbit(x) x&-x using namespace std; int read(){ char c;int s=0,t=1; while(!isdigit(c=getchar()))if(c=='-')t=-1; do{s=s*10+c-'0';}while(isdigit(c=getchar())); return s*t; } int min(int a,int b){return a<b?a:b;} int max(int a,int b){return a<b?b:a;} int ab(int x){return x>0?x:-x;} //int MO(int x){return x>=MOD?x-MOD:x;} //void insert(int u,int v){tot++;e[tot].v=v;e[tot].from=first[u];first[u]=tot;} /*------------------------------------------------------------*/ const int inf=0x3f3f3f3f; const int maxn=300010;//2^18!!! const double PI=acos(-1); int n,m; struct cp{ double x,y; cp(double a,double b){x=a;y=b;} cp(){x=y=0;}; cp operator + (cp a){return cp(x+a.x,y+a.y);} cp operator - (cp a){return cp(x-a.x,y-a.y);} cp operator * (cp a){return cp(x*a.x-y*a.y,x*a.y+y*a.x);} }a[maxn],b[maxn]; void fft(cp *a,int n,int f){ int k=0; for(int i=0;i<n;i++){ if(i>k)swap(a[i],a[k]); for(int j=n>>1;(k^=j)<j;j>>=1); } for(int l=2;l<=n;l<<=1){ int m=l/2; cp wn(cos(2*PI*f/l),sin(2*PI*f/l)); for(cp *p=a;p!=a+n;p+=l){ cp w(1,0); for(int i=0;i<l/2;i++){ cp t=w*p[i+m]; p[i+m]=p[i]-t; p[i]=p[i]+t; w=w*wn; } } } if(f==-1){for(int i=0;i<n;i++)a[i].x/=n;}// } int main(){ scanf("%d%d",&n,&m);n++;m++; for(int i=0;i<n;i++){int u;scanf("%d",&u);a[i]=cp(u,0);} for(int i=0;i<m;i++){int u;scanf("%d",&u);b[i]=cp(u,0);} int N=1;while(N<n+m)N*=2; 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); for(int i=0;i<n+m-1;i++)printf("%d ",(int)(a[i].x+0.1)); return 0; }
注意：
1.數組空間必須是大於兩個卷積數組長度和的2的冪（記為N），當n為1e5時數組空間為300000。
2.DFT后的操作一定要以N為單位，例如點值相乘。這個問題在多重fft的題目中很容易寫錯，多重fft必須把上一個的結果后半部分清零再繼續。
3.代碼需要手寫復數模板來減小常數，NTT還要預處理omega。枚舉反二進制的方法是從高位開始模擬進位。
（八）卷積
對於函數f(n)和g(n)，定義其卷積為函數(f⊗g)
(fg)(n)=Σf(i)g(n-i)，i=0~n。——形式冪級數
卷積的形式和多項式乘法類似，(fg)的生成函數是f(n)和g(n)的生成函數的乘積。
卷積是和為定值的形式，若差為定值，將其中一個數組反轉后即可卷積。卷積中多余的部分數組為0，不影響答案。
高精度乘法：將數字從低位到高位編號0~len-1，每一位代入多項式系數，那么數字乘法就是多項式乘法，最后從低到高處理進位。
（九）模數任意——fft合並
將每個數字拆成a*32768+b，然后四次dft后4次idft合並。
#include<cstdio> #include<cstring> #include<cctype> #include<cmath> #include<queue> #include<stack> #include<set> #include<vector> #include<complex> #include<algorithm> #define ll long long #define lowbit(x) x&-x using namespace std; int read(){ char c;int s=0,t=1; while(!isdigit(c=getchar()))if(c=='-')t=-1; do{s=s*10+c-'0';}while(isdigit(c=getchar())); return s*t; } int min(int a,int b){return a<b?a:b;} int max(int a,int b){return a<b?b:a;} int ab(int x){return x>0?x:-x;} //int MO(int x){return x>=MOD?x-MOD:x;} //void insert(int u,int v){tot++;e[tot].v=v;e[tot].from=first[u];first[u]=tot;} /*------------------------------------------------------------*/ const int inf=0x3f3f3f3f,MOD=1e9+7,maxn=300010;// const double PI=acos(-1); namespace fft{ complex<double>o[maxn],oi[maxn]; void init(int n){ for(int k=0;k<n;k++){o[k]=complex<double>(cos(2*PI*k/n),sin(2*PI*k/n));oi[k]=conj(o[k]);} } void transform(complex<double>*a,int n,complex<double>*o){ int k=0; while((1<<k)<n)k++; for(int i=0;i<n;i++){ int t=0; for(int j=0;j<k;j++)if(i&(1<<j))t|=(1<<(k-j-1)); if(i<t)swap(a[i],a[t]); } for(int l=2;l<=n;l*=2){ int m=l/2; for(complex<double>*p=a;p!=a+n;p+=l){ for(int i=0;i<m;i++){ complex<double>t=o[n/l*i]*p[i+m]; p[i+m]=p[i]-t; p[i]+=t; } } } } void dft(complex<double>*a,int n){transform(a,n,o);} void idft(complex<double>*a,int n){ transform(a,n,oi); for(int i=0;i<n;i++)a[i]/=n; } } int n,N,m,kind,f[maxn][30],h[30],F[maxn],g[2][maxn]; complex<double>a[maxn],b[maxn],c[maxn],d[maxn],v[maxn]; ll ans1[maxn],ans2[maxn],ans3[maxn],ans4[maxn];// void multply(complex<double>*x,complex<double>*y,ll *z){ for(int i=0;i<N;i++)v[i]=x[i]*y[i]; fft::idft(v,N); for(int i=0;i<n;i++)z[i]=(ll)(v[i].real()+0.5); } void MTT(int *x,int *y,int *z){ for(int i=0;i<N;i++)a[i]=b[i]=c[i]=d[i]=complex<double>(0,0);// for(int i=0;i<n;i++)a[i].real(x[i]>>15); for(int i=0;i<n;i++)b[i].real(x[i]&32767); for(int i=0;i<n;i++)c[i].real(y[i]>>15); for(int i=0;i<n;i++)d[i].real(y[i]&32767); fft::dft(a,N);fft::dft(b,N);fft::dft(c,N);fft::dft(d,N); multply(a,c,ans1);multply(a,d,ans2); multply(b,c,ans3); multply(b,d,ans4); for(int i=0;i<n;i++)z[i]=(ans1[i]*32768%MOD*32768%MOD+(ans2[i]+ans3[i])*32768%MOD+ans4[i])%MOD; } int main(){ n=read();m=read();kind=read(); f[0][0]=1;h[0]=1; int mx=0; for(mx=1;h[mx-1]<n;mx++)h[mx]=h[mx-1]*m;mx--; for(int i=1;i<=n;i++){ for(int j=0;j<=mx;j++)if(i-h[j]>=0){ for(int k=0;k<=j;k++)f[i][j]=(f[i][j]+f[i-h[j]][k])%MOD; } } for(int i=0;i<=n;i++)for(int j=0;j<=mx;j++)F[i]=(F[i]+f[i][j])%MOD; int x=0; for(int i=0;i<=n;i++)g[x][i]=F[i]; n++; N=1;// while(N<n+n)N*=2;fft::init(N); for(int k=1;k<kind;k++)MTT(g[x],F,g[1-x]),x=1-x; int sum=0; for(int i=0;i<=n;i++)sum=(sum+g[x][i])%MOD; printf("%d",sum); return 0; }
注意清空complex的時候實部和虛部(imag)要一起清空。
好看的模板：L_0_Forever_LF
【快速數論變換】NTT
只記錄重要概念，證明略去。
（一）原根
當(a,m)=1時，對於滿足a^x=1(%m)的最小正整數x，稱x為a模m的階。
根據歐拉定理a^φ(m)=1(%m)，當x=φ(m)時，稱a為m的原根。
以下只討論m為素數的情況，則當a為m的原根時，a^0~a^(p-2)取遍1~p-1所有值。
模m有原根的充要條件：m=2,4,p^e,2*p^3，p是奇素數。（也就是說，m為素數時一定有原根）
求m的原根：p-1= p1^a1 * p2^a2 * pk^ak，g是p的原根當且僅當對於所有的pi滿足g^[ (p-1)/pi ] ≠ 1 (%p)
（二）快速數論變換
當模數為形如p=r*2^k+1的素數（費馬素數）時，則有n|p-1，可以進行NTT。
先找到p的原根g（p=998244353 || 1004535809，g=3）
在原來FFT的基礎上，omega[i]=g^[ (p-1)/n*i ] % p，倒數為逆元。
IDFT時，除以n改為乘n的逆元。
（三）模數任意的NTT
找到三個費馬素數滿足相乘結果>n*(m-1)^2，分別進行NTT后用CRT合並。
p=998244353,1004535809,469762049，g=3。
我寫的是FFT合並。
（四）離散對數
當(a,p)=1時，若滿足a^x=b (%p)，則稱在模p意義下，x是b以a為底的離散對數，即logab=x（單個快速求解可用BSGS算法）。
1.對於x*y=z(%p)，有log x+log y=log z(%p-1)，因此離散對數常用於乘法轉加法（生成函數）。
2.對於x^y=z(%p)，有y*log x=log z(%p-1)。
其中log以p的原根g為底。
例題：【BZOJ】3992: [SDOI2015]序列統計 NTT+生成函數
【生成函數】母函數
生成函數的三大要素：①選擇項，②大小，③元素個數。一般最終要求某個“大小”的元素個數。
對於一類組合對象構成的集合A：
1.每個元素a∈A都定義了一個非負整數的”大小“，記為|a|。
2.大小為n的元素個數記為$A_n$。
那么A的一般生成函數是
$$A(x)=\sum_{i=0}^{n}A_ix^i$$
在這里每個元素都抽象成”大小“，元素a可以理解為有|a|的單位元素的元素。
組合對象集合D為A和B的笛卡爾積，即D中的每個元素都是A中某個元素a和B中某個元素b組成的有序二元組(a,b)，那么顯然有D(x)=A(x)B(x)。
若干一般生成函數的乘積中，第n項的含義是：每個選擇項取一個元素，大小相加為n的元素個數。
每個生成函數本質上是一個集合，那么若干生成函數的乘積就是★每個集合取一個元素的組合。例如生成函數A,B,C，A*B*C的每個元素就是有序三元組(a,b,c)。
指數型生成函數是
$$A(x)=\sum_{i=0}^{n}A_i\frac{x^i}{i!}$$
這樣就會有：
$$D_n=\sum_{i+j=n}A_iB_j\frac{(i+j)!}{i!j!}=\sum_{i+j=n}A_iB_j\binom{i+j}{i}$$
這里乘(i+j)!是因為這只是系數，后面要除以(i+j)!。
理解為每個元素內部有序就可以了，這樣元素內部是排列。
生成函數都是處理對於n個選擇項各選一個組成對應”大小“的元素個數，而一般生成函數元素內部是組合，指數型生成函數元素內部是排列。
一般生成函數還有個化簡公式，令x∈[-1,1]時套等比數列公式即可收斂：
$$\sum_{i=0}^{n}x^i=\frac{1}{1-x}$$
指數型生成函數也有個化簡公式——泰勒展開：
$$\sum_{i=0}^{n}\frac{x^i}{i!}=e^x$$
例題：
1.熱身：蘋果只能取偶數個，橘子只能取1~4個，求拿n個水果的方案數
題解：定義每種水果為一個集合（每個集合選一個），“大小”為水果個數，最后求總集合”大小“為n的元素個數。
f(x)=(1+x^2+x^4+...)*(1+x+x^2+x^3+x^4)。
如果覺得集合很難理解，不妨用”選擇項“這個詞。
題意：給定n個物品，價值為ai，物品價格互不相同，求選一個或兩個或三個的價值為x的方案數，輸出所有存在的x和對應方案數。ai<=40000。
題解：要求什么就定義什么為”大小“，所以定義”大小“為價值，[第一個物品][第二個物品][第三個物品]為三個選擇項。
那么每個選擇項的每個系數記錄對應價值的物品數量（1個）。
這樣拼起來就好了嗎？不是，物品不能重復取，所以拼起來之后再容斥掉選相同的。
我們可以直接寫出選一個物品的集合的生成函數f，兩個相同物品的g和三個相同物品的h。
考慮有AAB,ABA,BAA,AAA四種不合法情況，答案就是f^3-3f*g+2h。最后這個求得排列數，需要/3!。選1個或2個的隨便推推也一樣。
3.【BZOJ】3992: [SDOI2015]序列統計 NTT+生成函數
題意：給定一個[0,m-1]范圍內的數字集合S，從中選擇n個數字（可重復）構成序列。給定x，求序列所有數字乘積%m后為x的序列方案數%1004535809。1<=n<=10^9，3<=m<=8000，m為素數，1<=x<=m-1。
題解：要求乘積，定義”大小“為數字的乘積。但是我們不能加減”大小“啊？
換成離散對數就可以了，然后定義每個數字為選擇項，答案就是f^n。
題意：有n種物品，並且知道每種物品的數量。要求從中選出m件物品的排列數。例如有兩種物品A,B，並且數量都是1，從中選2件物品，則排列有”AB”,”BA”兩種。
題解：定義”大小“為物品數量，選擇項為每種物品，那么組合數就是一般生成函數（元素內部有序，嚴格按物品編號排），排列數就是指數型生成函數（元素內部帶標號，可以打亂）。
這里指數型生成函數，指的是方案就是元素內部帶標號的方案數，這樣再計算過程中那個公式自動/i!后進行計算再乘回(i+j)!的。
暴力枚舉求解這個生成函數。
另一部分知識：生成函數公式的化簡
參考：什么是生成函數？ by M67
首先根據二項式定理：
$$(1+x)^n=\sum_{i=0}^{n}\binom{n}{i}x^i$$
擴展到負數和實數，即廣義二項式定理，只要將組合數表示成下降冪即可：
$$(1+x)^{n}=\sum_{i=0}^{\infty}\frac{n^{\underline{i}}}{i!}x^i$$
這里有一個特殊的變換，當n>0時：
$$(1-x)^{-n}=\sum_{i=0}^{\infty}\frac{(-n)^{\underline{i}}}{i!}(-x)^i$$
$$(1-x)^{-n}=\sum_{i=0}^{\infty}\frac{(-1)^i*n^{\overline{i}}}{i!}*(-1)^i*x^i$$
$$(1-x)^{-n}=\sum_{i=0}^{\infty}\binom{n+i-1}{n-1}*x^i$$
所以這是將i個相同的數分割n個非空部分的方案數的生成函數。
另外還常用等比數列遞推公式來收斂等比數列：
$$\sum_{i=0}^{\infty}(x^q)^i=\frac{1}{1-x^q}$$
另外對於有限的等比數列還常用類似錯位相減的方法，左邊乘上(1-x)就會變成右邊分子，即：
$$\sum_{i=0}^{n}x^i=\frac{1-x^{n+1}}{1-x}$$
例題：
1.求：
$$g(x)=(1+x^2+x^4+...)(1+x^5+x^{10}+...)(1+x+x^2+x^3+x^4)(1+x)$$
用上面提到的技巧即可得到：
$$g(x)=\frac{1}{1-x^2}*\frac{1}{1-x^5}*\frac{1-x^5}{1-x}*(1+x)$$
不斷地約分，最后用平方差公式化簡可得：
$$g(x)=\frac{1}{(1-x)^2}=(1-x)^{-2}=\sum_{i=0}^{\infty}(i+1)x^i$$
再代入上面那個結論就可以得到多項式了。
2.食物：用上面的技巧化簡約分最后剩個小組合數，n=10^500用讀入取模。
【多項式求逆】
核心原理是$\%x^{\frac{n}{2}}$的多項式平方后可以轉化為$\%x^n$。
已知多項式f(x)，求多項式g(x)，滿足f(x)g(x)=1(%x^n)。
$$f(x)g(x)=1(\%x^{\frac{n}{2}})$$
$$f(x)h(x)=1(\%x^n)$$
現在已知g(x)，求h(x)。
$$f(x)g(x)-1=0(\%x^{\frac{n}{2}})$$
將1移到左邊后就可以平方了。
$$f(x)^2g(x)^2-2f(x)g(x)+1=0(\%x^n)$$
將1換成f(x)*g(x)，從而將h(x)代入，提出f(x)消去。
$$f(x)g(x)^2-2g(x)+h(x)=0(\%x^n)$$
最終得到：
$$h(x)=g(x)(2-f(x)g(x))(\%x^n)$$
然后就可以遞歸求解，邊界條件:當n=1時，f(x)g(x)=1(%x)，g(0)是f(0)在%MOD意義下的數論逆元。
注意：每次n不同都要重新預處理Omega[]。
例題：【BZOJ】4555: [Tjoi2016&Heoi2016]求和 排列組合+多項式求逆 或 斯特林數+NTT
【拉格朗日插值法】
參考：拉格朗日插值法(圖文詳解) by Angel_Kitty
對於一個n次多項式，如果已知n+1個點，可以構造拉格朗日多項式L(x)：
$$L(x)=\sum_{i=0}^{n}y_il_i(x)$$
其中$l_i(x)$為插值基函數：
$$l_i(x)=\prod_{j=0,j\neq i}^{n}\frac{x-x_j}{x_i-x_j}$$
通過代入需要的x即可得到答案。
每次插值的復雜度為O(n^2)。
例題：【BZOJ】4559: [JLoi2016]成績比較 計數DP+排列組合+拉格朗日插值
當橫坐標連續時，上式可以表示為：
$$f(x)=\sum_{i=0}^{n}f(i)*\prod_{j=0,j\neq i}^{n}\frac{x-j}{i-j}$$
預處理階乘的逆元$\frac{1}{i!}$，預處理$v=\prod x-j$，每次乘上v/(i-j)，分母是兩段階乘，再根據負數的個數判斷正負性。
復雜度O(n)。
例題：【BZOJ】3453: tyvj 1858 XLkxc 拉格朗日插值（自然數冪和）