标签：
今天看了算法导论，对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; }
最后我们得到一个定理：
标签：
原文地址：http://www.cnblogs.com/maijing/p/4787535.html
