一、 多项式的表示法:
1.系数表示法:
一个(n-1)的(n)项多项式(f(x))可以表示为(sumlimits_{i=1}^{n-1}a_i*x^i)。
于是可以用每一项的系数表示:(f(x)={a_0,a_1,...a_{n-1}})。
2.点值表示法:
将多项式看成一个函数,那么找(n)个不同的数(x)代入,可以得到(n)个不同的(y)。
显然这(n)对((x,y))可以唯一确定一个多项式,因为这相当于(n)个方程的线性方程组。
写出来就是:(f(x)={(x_0,f(x_0)),(x_1,f(x_1)),...,(x_{n-1},f(x_{n-1}))})。
二、多项式乘法:
1.朴素算法:
我们现在要求(h(x)=f(x)*g(x))。
我们分别考虑上两种表示法做乘法的复杂度:
<1>系数表示法:我们要枚举(f(x))和(g(x))的每一位系数相乘,复杂度显然为(O(n^2))。
<2>点值表示法:两个点值表示法的多项式相乘是(O(n))的!!!
举个例子:
(f(x)={(x_0,f(x_0)),(x_1,f(x_1)),...,(x_{n-1},f(x_{n-1}))})。
(g(x)={(x_0,g(x_0)),(x_1,g(x_1)),...,(x_{n-1},g(x_{n-1}))})。
那么我们只需要(O(n))枚举每一位,便可得到:
(h(x)={(x_0,f(x_0)*g(x_0)),(x_1,f(x_1)*g(x_1)),...,(x_{n-1},f(x_{n-1})*g(x_{n-1}))})。
但是我们从系数表示转点值表示是(O(n^2))的,从点值转系数还要(O(n^3))消元。
系数表示法似乎无法优化,但是点值表示似乎有优化的空间。
2.复数
见OI-WIKI,讲的很详细。
3.单位复根
1.定义:
以下的(n)都是(2)的整次幂。
我们在系数表示转点值表示时会找(n)个不同的(x_i)代入,对于每个(x_i),我们都要算出(x_i^0,x_i^1,...,x_i^{n-1}),这显然是(O(n^2))的。
我们考虑取一些特殊的(x)优化这个过程,比如我们可以找(+1,-1,+i,-i),这些数的若干次方都是(1),但是只找这(4)个显然是不够的。
现在我们在实轴和虚轴形成的坐标系中画一个单位圆,那么这个圆上的所有点都满足我们的要求,即若干次方后为(1)。
我们再将这个圆(n)等分,从((1,0))开始标号,分别标为(0,1,2,...,n-1)。
比如(n=8):
设(omega_n^k)表示(n)等分时标号为(k)的点代表的值,根据复数的运算(模长相乘,角相加)可以得到:(omega_n^k=(omega_n^1)^k)。
我们将(omega_n^1)称为(n)次单位根。
根据三角函数的公式可以推出(omega_n^k)的公式:(omega_n^k=cosfrac{k}{n}2pi+isinfrac{k}{n}2pi)。
于是(omega_n^0,omega_n^1,..,omega_n^{n-1})就是我们要代入的(x_0,x_1,..,x_{n-1})。
2.性质:
<1>(omega_n^k=omega_{2n}^{2k})
证明只需代入公式。
<2>(omega_n^{k+frac{n}{2}}=-omega_n^k)
证明:
(w_n^{frac{n}{2}}=cosfrac{n}{2}frac{2pi}{n}+isinfrac{n}{2}frac{2pi}{n})
(=cospi+i*sinpi)
(=-1)
<3>(omega_n^n=omega_n^0=1)
现在我们找到了一组特殊的(x),但是我们依旧要将每个(x_i)暴力代入算,所以复杂度仍然为(O(n^2))。
三、FFT
上面那种方法叫做离散傅里叶变换((DFT)),我们考虑分治。
假设我们有个多项式为(f(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1})。
我们按照奇偶分类并将奇数类都提出一个(x):
(f(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+(a_1x^1+a_3x^3...+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}))
此时我们定义另外两个多项式:
(f1(x)=a_0+a_2x^1+a_4x^2+...+a_{n-2}x^{frac{n}{2}-1})
(f2(x)=a_1+a_3x^1+...+a_{n-1}x^{frac{n}{2}-1})
我们发现(f(x)=f1(x^2)+xf2(x^2))!!!!
设(k<frac{n}{2}):
将(omega_n^k)代入:
(f(omega_n^k)=f1((omega_n^k)^2)+omega_n^kf2((omega_n^k)^2))
(=f1(omega_n^{2k})+omega_n^kf2(omega_n^{2k}))
(=f1(omega_{frac{n}{2}}^k)+omega_n^kf2(omega_{frac{n}{2}}^{k}))
将(omega_n^{k+frac{n}{2}})代入:
(f(omega_n^{k+frac{n}{2}})=f1(omega_n^{2k+n})+omega_n^{k+frac{n}{2}}f2(omega_n^{2k+n}))
(=f1(omega_n^{2k}omega_n^n)-omega_n^kf2(omega_n^{2k}omega_n^n))
(=f1(omega_n^{2k})+omega_n^kf2(omega_n^{2k}))
(=f1(omega_{frac{n}{2}}^k)-omega_n^kf2(omega_{frac{n}{2}}^{k}))
上面两个式子的结果除了符号不同以外完全相同,于是我们只要知道(f1(omega_{frac{n}{2}}^k),f2(omega_{frac{n}{2}}^{k})),我们就能求出
(f(omega_n^k),f(omega_n^{k+frac{n}{2}}))。
于是就可以递归分治了,复杂度显然是(O(nlog_2n))的。
四、IFFT
我们用(FFT)求出了两个多项式乘积的点值表示((x_i,y_i)),我们还要将其转为系数表示。
考虑怎么将点值表示快速转为系数表示:
这里有个结论:
一个多项式在分治的过程中乘上单位根的共轭复数,分治完的每一项除以n即为原多项式的每一项系数。
我们设(c_k=sumlimits_{i=0}^{n-1}y_i(omega_n^{-k})^i)。
那么((c_i,w_n^{-i}))就是多项式(H(x)=y_0+y_1x+...+y_{n-1}x^{n-1})的点值表示。
先推一波式子:
(c_k=sumlimits_{i=0}^{n-1}y_i(omega_n^{-k})^i)
(=sumlimits_{i=0}^{n-1}(sumlimits_{j=0}^{n-1}a_j(omega_n^i)^j)(omega_n^{-k})^i)
(=sumlimits_{i=0}^{n-1}(sumlimits_{j=0}^{n-1}a_j(omega_n^j)^i)(omega_n^{-k})^i)
(=sumlimits_{i=0}^{n-1}(sumlimits_{j=0}^{n-1}a_j(omega_n^j)^i(omega_n^{-k})^i))
(=sumlimits_{i=0}^{n-1}(sumlimits_{j=0}^{n-1}a_j(omega_n^{j-k})^i))
设(S(x)=sumlimits_{i=0}^{n-1}x^i)。
当(k!=0)时:
代入(omega_n^k)可得:
(S(omega_n^k)=1+omega_n^k+(omega_n^k)^2+...+(omega_n^k)^{n-1})
两边同乘(omega_n^k)得:
(omega_n^kS(omega_n^k)=omega_n^k+(omega_n^k)^2+...+(omega_n^k)^{n})
两式相减可得:
((omega_n^k-1)S(omega_n^k)=(omega_n^k)^n-1)
于是可得:
(S(omega_n^k)=frac{(omega_n^k)^n-1}{omega_n^k-1})
(=frac{(omega_n^n)^k-1}{omega_n^k-1})
(=frac{1-1}{omega_n^k-1})
(=0)
当(k=0)时:
显然(S(omega_n^k)=n)
回到之前的式子:
(=sumlimits_{i=0}^{n-1}(sumlimits_{j=0}^{n-1}a_j(omega_n^{j-k})^i))
当(j!=k)时中间那个式子都是(0)
于是:
(c_k=na_k)
(a_k=frac{c_k}{n})
证完了。
五、迭代FFT
递归的常数太大,我们需要非递归的方法。
先放张经典图:
我们发现:
每个位置分治后的最终位置为其二进制翻转后得到的位置。
于是我们预先处理好每个数的最终位置,之后向上合并即可。
模板题
code:
#include<bits/stdc++.h>
using namespace std;
const int maxn=1e7+10;
const double Pi=acos(-1.0);
int n,m,lim=1,len;
int pos[maxn];
struct cplx{double x,y;}a[maxn],b[maxn];
cplx operator+(cplx a,cplx b){return (cplx){a.x+b.x,a.y+b.y};}
cplx operator-(cplx a,cplx b){return (cplx){a.x-b.x,a.y-b.y};}
cplx operator*(cplx a,cplx b){return (cplx){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
inline int read()
{
char c=getchar();int res=0,f=1;
while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
while(c>='0'&&c<='9')res=res*10+c-'0',c=getchar();
return res*f;
}
void fft(cplx* a,int op)
{
for(int i=0;i<lim;i++)if(i<pos[i])swap(a[i],a[pos[i]]);
for(int mid=1;mid<lim;mid<<=1)
{
cplx wn=(cplx){cos(Pi/mid),op*sin(Pi/mid)};
for(int i=0,r=mid<<1;i<lim;i+=r)
{
cplx w=(cplx){1,0};
for(int j=0;j<mid;j++,w=w*wn)
{
cplx x=a[i+j],y=w*a[i+mid+j];
a[i+j]=x+y;a[i+mid+j]=x-y;
}
}
}
}
int main()
{
n=read(),m=read();
for(int i=0;i<=n;i++)a[i].x=read();
for(int i=0;i<=m;i++)b[i].x=read();
while(lim<=n+m)lim<<=1,len++;
for(int i=0;i<lim;i++)pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
fft(a,1);fft(b,1);
for(int i=0;i<lim;i++)a[i]=a[i]*b[i];
fft(a,-1);
for(int i=0;i<=n+m;i++)printf("%d ",(int)(a[i].x/lim+0.5));
return 0;
}