FFT
今天看了算法导论，对FFT感受颇深。
感觉我就在抄算法导论。
回归正题。
多项式的表示
系数表示法
系数表示法其实非常常见，其实就是：
$$A(x)=\sum\limits_{j=0}^{n-1}a_jx^j$$
这是一个n次多项式，每个项的次数为0,1,...,n-1
如果用系数表示法来做多项式加法，那么时间复杂度是$O(N)$的
但是用系数表示法来做多项式乘法，那么时间复杂度是$O(N^2)$的，这样非常慢，所以这里介绍的FFT可以把多项式乘法的时间复杂度优化到$O(Nlog_2N)$
点值表示法
一个n次多项式$A(x)$的点值表示就是n个点值对所形成的集合：
$$\{(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1})\}$$
其中所有的$x_k$各不相同，并且当k=0,1,...,n-1时有
$$y_k=A(x_k)$$
一个多项式可以有很多种点值表示，只要选取的n个$x_k$互异即可。
对于许多关于多项式的操作来说，点值表示法是很方便的。
如果用点值表示法来做多项式加法，那么时间复杂度是$O(N)$的：
如果$C(x)=A(x)+B(x)$，则对于任意的$x_k$，$C(x_k)=A(x_k)+B(x_k)$。准确地说，如果已知$A(x)$的点集表示：
$$\{(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1})\}$$
和$B(x)$的点集表示
$$\{(x_0,y^{'}_0),(x_1,y^{'}_1),...,(x_{n-1},y^{'}_{n-1})\}$$
（注意，A和B对相同的n个点求值）则$C(x)$的点集表示为：
$$\{(x_0,y_0+y^{'}_0),(x_1,y_1+y^{'}_1),...,(x_{n-1},y_{n-1}+y^{'}_{n-1})\}$$
如果点值表示法来做多项式乘法，那么时间复杂度也是$O(N)$的：
如果$C(x)=A(x)B(x)$，则对于任意的$x_k$，$C(x_k)=A(x_k)B(x_k)$。准确地说，如果已知$A(x)$的点集表示：
$$\{(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1})\}$$
和$B(x)$的点集表示
$$\{(x_0,y^{'}_0),(x_1,y^{'}_1),...,(x_{n-1},y^{'}_{n-1})\}$$
（注意，A和B对相同的n个点求值）则$C(x)$的点集表示为：
$$\{(x_0,y_0y^{'}_0),(x_1,y_1y^{'}_1),...,(x_{n-1},y_{n-1}y^{'}_{n-1})\}$$
系数表示法转化成点值表示法称为求值；
点值表示法转化成系数表示法称为插值。
我们发现，对于多项式乘法，点值表示法比系数表示法有优势。
如果我们任选n个点求值（把系数表示法转化成点值表示法），那么时间复杂度是$O(N^2)$，但是稍后我们可以看到，如果我们巧妙地选取$x_k$，那么时间复杂度可以做到$O(Nlog_2N)$，并且插值（把点值表示法转化成系数表示法）也可以做到$O(Nlog_2N)$
于是我们得到一个思路：
单位复根
这里为稍后我们讨论奠定了一下基础。
n次单位复根是指满足$w^n=1$的复数$w$。n次单位复根有n个，分别是$e^{2\pi ik/n},k=0,1,..,n-1$
我们可以从欧拉公式理解：
$$e^{ix}=cosx+isinx$$
如图我们可以看到，n个单位复根平均分布在以原点为圆心的复平面上，其中
$$w_n=e^{2\pi i/n}$$
称为主n次单位复根，任何一个n次单位复根都是主n次单位复根$w_n$的幂
所以n个n次单位复根就是$w_n^0,w_n^1,...,w_n^{n-1}$
所以其实$w_n^k=e^{2\pi ik/n},k=0,1,..,n-1$
然后然后n次单位复根有一些性质：
引理一（相消引理）：对于任意整数$n\geq 0$，$k\geq 0$，$d>0$，有$w_{dn}^{dk}=w_n^k$
证明：$w_{dn}^{dk}=e^{2\pi idk/dn}=e^{2\pi ik/n}=w_n^k$
引理二：对于偶数$n>0$，有$w_{n}^{n/2}=w_2=-1$
引理三（折半引理）：对于偶数$n>0$，n个n次单位复根的平方等于n/2个n/2次单位复根。
引理四（求和引理）：对于任意整数$n\geq 1$和不能被n整除的非负整数k，有：
$$\sum\limits_{j=0}^{n-1}(w_{n}^{k})^j=0$$
证明：等比数列的求和公式也可以用于复数：
$\sum\limits_{j=0}^{n-1}(w_{n}^{k})^j=\frac{(w_{n}^{k})^{n}-1}{w_{n}^{k}-1}=\frac{(w_{n}^{n})^{k}-1}{w_{n}^{k}-1}=\frac{(1)^{k}-1}{w_{n}^{k}-1}=0$
从系数表示法到点值表示法
不失一般性，我们设n是2的幂，如果n不是2的幂，可以适当增大n。
前面已经提到，其实我们是对n个n次单位复根$w_n^0,w_n^1,...,w_n^{n-1}$进行求值，即：
如果
$$A(x)=\sum\limits_{j=0}^{n-1}a_jx^j$$
那么
$$y_k=A(w_n^{k})=\sum\limits_{j=0}^{n-1}a_jw_n^{kj},其中k=0,1,...,n-1$$
我们运用分治策略，用$A(x)$中偶数下标的系数和奇数下标的系数，分别定义了两个新的n/2次多项式$A^{[0]}(x)$和$A^{[1]}(x)$：
$A^{[0]}(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{n/2-1}$
$A^{[1]}(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{n/2-1}$
那么有如下式：
$$A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)$$
我们发现，$A^{[0]}(x)$和$A^{[1]}(x)$是一个子问题，我们递归处理。
现在我们已经知道了$A^{[0]}(x)$在$w_{n/2}^0,w_{n/2}^1,...,w_{n/2}^{n/2-1}$处的取值：
$$\{(w_{n/2}^0,y^{'}_{0}),(w_{n/2}^1,y^{'}_{1}),...,(w_{n/2-1}^{n/2-1},y^{'}_{n/2-1})\}$$
$$其中y^{'}_{k}=A^{[0]}(w_{n/2}^k)=A^{[0]}(w_{n}^{2k})，k=0,1,...,n/2-1$$
和$A^{[1]}(x)$在$w_{n/2}^0,w_{n/2}^1,...,w_{n/2}^{n/2-1}$处的取值：
$$\{(w_{n/2}^0,y^{''}_{0}),(w_{n/2}^1,y^{''}_{1}),...,(w_{n/2-1}^{n/2-1},y^{''}_{n/2-1})\}$$
$$其中y^{''}_{k}=A^{[1]}(w_{n/2}^k)=A^{[1]}(w_{n}^{2k})，k=0,1,...,n/2-1$$
我们想知道$A(x)$在$w_n^0,w_n^1,...,w_n^{n-1}$处的取值的求值：
$$\{(w_{n/2}^0,y_{0}),(w_{n/2}^1,y_{1}),...,(w_{n-1}^{n-1},y_{n-1})\}$$
运算一下：
$若k=0,1,...,n/2-1$
$y_{k}$
$=A(w_n^{k})$
$=A^{[0]}((w_n^{k})^2)+w_n^{k}A^{[1]}((w_n^{k})^2)$
$=A^{[0]}(w_n^{2k})+w_n^{k}A^{[1]}(w_n^{2k})$
$=y^{'}_{k}+w_n^{k}y^{''}_{k}$
$y_{k+n/2}$
$=A(w_n^{k+n/2})$
$=A^{[0]}((w_n^{k+n/2})^2)+w_n^{k+n/2}A^{[1]}((w_n^{k+n/2})^2)$
$=A^{[0]}(w_n^{2k+n})+w_n^{k+n/2}A^{[1]}(w_n^{2k+n})$
$=A^{[0]}(w_n^{2k})-w_n^{k}A^{[1]}(w_n^{2k})$
$=y^{'}_{k}-w_n^{k}y^{''}_{k}$
（其实就是上几行的等式，不懂可以往回看几行）
我们发现我们已经知道了$y^{'}_{k}$和$y^{''}_{k}$了！！！！！
所以我们完全可以算出$y_{0},y_{1},...y_{n-1}$了。
好了，到现在为止可以在$O(Nlog_2N)$的时间内完成从系数表示法到点值表示法的转化。
伪代码是这样的：
从点值表示法到系数表示法
我们通过矩阵来理解：
从左到右3个矩阵分别记为$Y$，$V$和$A$，其中$Y_{j,1}=y_j(j=0,1,...,n-1)$，$V_{i,j}=w_n^{ij}(i=0,1,...,n-1，j=0,1,...,n-1)$
我们可以表示为矩阵积：$Y=VA$
现在我们已经知道了$Y$和$V$，想求$A$
如果我们知道$V$的逆$V^{-1}$，那么我们可以求$A$：
$Y=VA$
$V^{-1}Y=V^{-1}VA$
$V^{-1}Y=A$
我们先给出一个结论：$V^{-1}_{i,j}=w_n^{-ij}/n(i=0,1,...,n-1，j=0,1,...,n-1)$
证明：
记$X=V^{-1}V$
$X_{i,j}=\sum\limits_{k=0}^{n-1}V^{-1}_{i,k}V_{k,j}=\sum\limits_{k=0}^{n-1}\frac{w_n^{-ik}w_n^{kj}}{n}=\sum\limits_{k=0}^{n-1}\frac{(w_n^{j-i})^k}{n}$
若$i=j$，那么$w_n^{j-i}=1$，所以$X_{i,j}=1$
若$i\neq j$，根据上面单位复根的求和引理，可知$X_{i,j}=0$
$X$是一个单位矩阵
所以$V^{-1}$是$V$的逆。
根据$V^{-1}Y=A$我们可以得到：
$a_k=\sum\limits_{j=0}^{n-1}V^{-1}_{k,j}y_j=\sum\limits_{j=0}^{n-1}w_n^{-kj}y_j/n=\frac{1}{n}\sum\limits_{j=0}^{n-1}y_jw_n^{-kj}$
即：
$$a_k=\frac{1}{n}\sum\limits_{j=0}^{n-1}y_jw_n^{-kj},其中k=0,1,...,n-1$$
我们来看一下”从系数表示法到点值表示法“时我们是怎么求值的：
$$y_k=\sum\limits_{j=0}^{n-1}a_jw_n^{kj},其中k=0,1,...,n-1$$
我们发现这两条式子惊人地相似，就是$y_k$换成了$a_k$，$a_j$换成了$y_j$，$w_n^{kj}$换成了$w_n^{-kj}$，还多了了一个$\frac{1}{n}$，不过这个很好处理。
我们完全可以沿用求值的函数就可以，将($y_0,y_1,...,y_{n-1})$传递进函数，然后$w_n^{kj}$变成$w_n^{-kj}$即可；最后还要除以$n$
多项式乘法的程序是这样的：
#include<cstdio> #include<cstdlib> #include<iostream> #include<fstream> #include<algorithm> #include<cstring> #include<string> #include<cmath> #include<queue> #include<stack> #include<map> #include<utility> #include<set> #include<bitset> #include<vector> #include<functional> #include<deque> #include<cctype> #include<climits> #include<complex> using namespace std; typedef long long LL; typedef double DB; typedef pair<int,int> PII; typedef complex<DB> CP; #define mmst(a,v) memset(a,v,sizeof(a)) #define mmcy(a,b) memcpy(a,b,sizeof(a)) #define re(i,a,b) for(i=a;i<=b;i++) #define red(i,a,b) for(i=a;i>=b;i--) #define fi first #define se second #define m_p(a,b) make_pair(a,b) #define SF scanf #define PF printf #define two(k) (1<<(k)) template<class T>inline T sqr(T x){return x*x;} template<class T>inline void upmin(T &t,T tmp){if(t>tmp)t=tmp;} template<class T>inline void upmax(T &t,T tmp){if(t<tmp)t=tmp;} const DB EPS=1e-9; inline int sgn(DB x){if(abs(x)<EPS)return 0;return(x>0)?1:-1;} const DB Pi=acos(-1.0); inline int gint() { int res=0;bool neg=0;char z; for(z=getchar();z!=EOF && z!='-' && !isdigit(z);z=getchar()); if(z==EOF)return 0; if(z=='-'){neg=1;z=getchar();} for(;z!=EOF && isdigit(z);res=res*10+z-'0',z=getchar()); return (neg)?-res:res; } inline LL gll() { LL res=0;bool neg=0;char z; for(z=getchar();z!=EOF && z!='-' && !isdigit(z);z=getchar()); if(z==EOF)return 0; if(z=='-'){neg=1;z=getchar();} for(;z!=EOF && isdigit(z);res=res*10+z-'0',z=getchar()); return (neg)?-res:res; } /* 注意事项 1.数组要开大 保证高位有足够的0 2.n必须是2的幂次，补0即可 3.最后所以的数要除n 4.n次多项式与m次多项式的乘积为n+m次多项式 长度为n的多项式与长度为m的多项式的乘积为长度为n+m-1的多项式 */ const int maxN=100000*4; int alen,blen,clen,N; CP a[maxN+1000],b[maxN+1000],c[maxN+1000]; CP t[maxN+1000]; inline void FFT(CP *a,int N,int l,int r,int flag) { if(l==r) return; int i,len=r-l+1,mid=(l+r)/2,ln=l,rn=mid+1; for(i=l;i<=r;i+=2)t[ln++]=a[i],t[rn++]=a[i+1]; re(i,l,r)a[i]=t[i]; FFT(a,N,l,mid,flag);FFT(a,N,mid+1,r,flag); CP wn(cos(2*Pi/DB(len)),sin(DB(flag)*2*Pi/DB(len))),w(1.0,0);//注意没有=号，常错 re(i,l,mid){CP x=a[i],t=w*a[i+len/2];a[i]=x+y;a[i+len/2]=x-y;w=w*wn;} if(flag==-1 && len==N)re(i,l,r)a[i]/=DB(N); } int main() { freopen("UOJ#34.in","r",stdin); freopen("UOJ#34.out","w",stdout); int i; alen=gint()+1;blen=gint()+1; re(i,0,alen-1)a[i]=DB(gint()); re(i,0,blen-1)b[i]=DB(gint()); clen=alen+blen-1; for(N=1;N<clen;N<<=1); FFT(a,N,0,N-1,1); FFT(b,N,0,N-1,1); re(i,0,N-1)c[i]=a[i]*b[i]; FFT(c,N,0,N-1,-1); re(i,0,clen-1)printf("%d ",int(c[i].real()+0.5));printf("\n"); return 0; }
最后我们得到一个定理：
模运算下的FFT
上面我们谈论FFT都是用实数的。
但有时候我们要将系数$a_i$模一个数$P$，准确地说：
$已知A(x)和B(x):$
$$A(x)=\sum\limits_{j=0}^{n-1}a_jx^j，其中0\leq a_j<P$$
$$B(x)=\sum\limits_{j=0}^{n-1}b_jx^j，其中0\leq b_j<P$$
$A(x)和B(x)的积为C(x)=A(x)B(x),那么$
$$c_j=\sum\limits_{k=0}^{j}a_kb_{j-k}\%P，其中0\leq j \leq2n-2$$
那么怎么做呢？
用原根（如果不知道原根，可以看我另一篇博客）
假设$P$是质数，那么$P$一定有原根（不妨设为$G$）。
设$P=2^{a}\times b+1,其中b为奇数$
我们知道：
$G^{\varphi (P)}\equiv 1(\% P)$
由费马小定理得：
$G^{P-1}\equiv 1(\% P)$
又因为$P=2^{a}\times b+1$，所以$P-1=2^{a}\times b$，即：
$G^{2^{a}\times b}\equiv 1(\% P)$
$(G^{b})^{2^{a}}\equiv 1(\% P)$
类似于复数的做法，我们重新定义$w_n=(G^{b})^{2^{a}/n}\%P，其中n是2的幂，且1 \leq n\leq 2^{a}$
那么$w_n^{k}=w_n=(G^{b})^{k*2^{a}/n}\%P，其中k=0,1,...,n-1$
我们发现此时任然满足n次单位复根的性质：
引理一（相消引理）：对于任意整数$n\geq 0$，$k\geq 0$，$d>0$，有$w_{dn}^{dk}=w_n^k$
引理二：对于偶数$n>0$，有$w_{n}^{n/2}=w_2=-1$
证明：
我们先考虑一个问题：
$求解x^2\equiv 1(\%P)，其中P是质数，且0\leq x< P$
变形得：
$(x^2-1)\%P=0$
$(x+1)(x-1)\%P=0$
因为$P$是质数
所以$x-1=0$或$x+1=P$，即$x=1或P-1$
所以$w_2=G^{b*2^{a-1}}$要么是$1$，要么是$P-1$
又因为$G$是原根
所以$G^0,G^1,...,G^{P-2}$和$1,2,...,P-1$是一一对应的。
显然$G^0$对应了$1$
所以$G^{b*2^{a-1}}$不能对应$1$，只能对应$P-1$
所以$w_2=G^{b*2^{a-1}}=P-1$
又因为在模$P$意义下，$-1$和$P-1$是一样的
所以$w_{n}^{n/2}=w_2=-1$
引理三（折半引理）：对于偶数$n>0$，$n$个$w_n^{k}(k=0,1,...,n-1)$的平方等于$n/2$个$w_{n/2}^{k}(k=0,1,...,n/2-1)$
引理四（求和引理）：对于任意整数$n\geq 1$和不能被n整除的非负整数k，有：
$$\sum\limits_{j=0}^{n-1}(w_{n}^{k})^j=0$$
我们发现居然有和n次单位复根的性质一模一样的性质。
我们只需要改动一下即可。
因此，若满足$2^{k}|\varphi(p)且2^{k}\geq n$,我们可以用$g^{ \frac{\varphi(p)}{2^k}}$作为单位根，其中$g$是$p$的原根
如果$p$没有原根呢？我们可以将$p$分解质因数，取若干个便于FFT的大质数进行取模，然后再用中国剩余定理还原系数即可。
#include<cstdio> #include<cstdlib> #include<iostream> #include<fstream> #include<algorithm> #include<cstring> #include<string> #include<cmath> #include<queue> #include<stack> #include<map> #include<utility> #include<set> #include<bitset> #include<vector> #include<functional> #include<deque> #include<cctype> #include<climits> #include<complex> //#include<bits/stdc++.h>适用于CF,UOJ，但不适用于poj using namespace std; typedef long long LL; typedef double DB; typedef pair<int,int> PII; typedef complex<DB> CP; #define mmst(a,v) memset(a,v,sizeof(a)) #define mmcy(a,b) memcpy(a,b,sizeof(a)) #define fill(a,l,r,v) fill(a+l,a+r+1,v) #define re(i,a,b) for(i=(a);i<=(b);i++) #define red(i,a,b) for(i=(a);i>=(b);i--) #define ire(i,x) for(typedef(x.begin()) i=x.begin();i!=x.end();i++) #define fi first #define se second #define m_p(a,b) make_pair(a,b) #define p_b(a) push_back(a) #define SF scanf #define PF printf #define two(k) (1<<(k)) template<class T> T sqr(T x){return x*x;} template<class T> void upmin(T &t,T tmp){if(t>tmp)t=tmp;} template<class T> void upmax(T &t,T tmp){if(t<tmp)t=tmp;} const DB EPS=1e-9; int sgn(DB x){if(abs(x)<EPS)return 0;return(x>0)?1:-1;} const DB Pi=acos(-1.0); int gint() { int res=0;bool neg=0;char z; for(z=getchar();z!=EOF && z!='-' && !isdigit(z);z=getchar()); if(z==EOF)return 0; if(z=='-'){neg=1;z=getchar();} for(;z!=EOF && isdigit(z);res=res*10+z-'0',z=getchar()); return (neg)?-res:res; } LL gll() { LL res=0;bool neg=0;char z; for(z=getchar();z!=EOF && z!='-' && !isdigit(z);z=getchar()); if(z==EOF)return 0; if(z=='-'){neg=1;z=getchar();} for(;z!=EOF && isdigit(z);res=res*10+z-'0',z=getchar()); return (neg)?-res:res; } const int maxm=8000; const LL MOD=1004535809; const LL G=3; int n,m,X,cnt,d; LL power(LL a,LL k,LL Mod){LL x=1;while(k){if(k&1)x=x*a%Mod;a=a*a%Mod;k>>=1;}return x;} struct GF { LL a[4*maxm+200]; LL& operator [](int i){return a[i];} GF(){mmst(a,0);} void FFT(LL *a,int n,int type) { if(n==1)return; static LL t[4*maxm+200]; int i,ln=0,rn=n>>1; for(i=0;i<=n-1;i+=2)t[ln++]=a[i],t[rn++]=a[i+1]; re(i,0,n-1)a[i]=t[i]; LL *l=a,*r=a+(n>>1); FFT(l,n>>1,type); FFT(r,n>>1,type); LL wn=power(G,(LL)type*(MOD-1)/n%(MOD-1),MOD),w=1; for(i=0;i<n>>1;i++,w=w*wn%MOD)t[i]=l[i]+w*r[i],t[i+(n>>1)]=l[i]-w*r[i]; re(i,0,n-1)a[i]=(t[i]%MOD+MOD)%MOD; } GF& operator *=(GF &f) { static LL b[4*maxm+200]; int i; re(i,0,d-1)b[i]=f[i]; FFT(a,d,1); FFT(b,d,1); re(i,0,d-1)a[i]=a[i]*b[i]%MOD; FFT(a,d,MOD-2); re(i,m-1,(m-2)<<1)(a[i-(m-1)]+=a[i])%=MOD,a[i]=0; LL inv=power(d,MOD-2,MOD); re(i,0,m-2)a[i]=a[i]*inv%MOD; return *this; } };
