多项式的表示
1 系数表示法
设是一个关于的次多项式,则:
2 点值表达法
我们可以把次多项式看作一个函数,那么它可以用平面直角坐标系上的个点来确定。我们把这个点代入,我们就可以得到一个元一次方程组,然后通过高斯消元就可以确定这个多项式。
系数表达法的多项式乘法时间复杂度显然是的,但是点值表达法的多项式的乘法的时间复杂度却是的(两个多项式的每一个点的横坐标都相等)。那我们就会希望可以利用点值表达法的这一优势来快速地进行多项式乘法。但是我们平时使用的都是系数表示法,想要利用这个优势,就要快速地将系数表达法转化为点值表达法。
复数
复数的定义
我们把形如(,均为实数)的数称为复数,其中称为实部,称为虚部,称为虚数单位。当虚部等于零时,这个复数可以视为实数;当的虚部不等于零时,实部等于零时,常称为纯虚数。复数域是实数域的代数闭包,也即任何复系数多项式在复数域中总有根。 复数是由意大利米兰学者卡当在十六世纪首次引入,经过达朗贝尔、棣莫弗、欧拉、高斯等人的工作,此概念逐渐为数学家所接受。(摘自百度百科)
复数的运算法则
1 复数的加法
设,那么:
2 复数的乘法
设,那么:
3 (一个奇奇怪怪的东西)
单位根
定义
如果,那么我们称为次单位根。
单位根的值
事实上,在这里:
由欧拉恒等式得到。
消去引理
证明:
折半引理
如果是大于的一个偶数,那么前个次单位根的平方的集合等于后个次单位根的平方的集合。
证明:对于任意的,有:
因此,
DFT 离散傅里叶变换
终于切入正题了。
DFT就是将一个多项式的系数表达法转化为点值表达法。
具体操作:
对于次多项式,我们有:
令
则
把代入可得到:
由折半引理我们会发现
易证
所以
我们会发现这两个式子只有常数项不同
那么当我们在求第一个式子的时候,我们可以直接在的时间内同时得到第二个式子的值
然后就这样递归下去,就可以得到的点值表达式了
IDFT 离散傅里叶逆变换
我们用点值表达法求出了两个多项式的乘积,但是我们还需要将这个点值表达法转换回系数表达法。那这个又该怎么操作呢?
首先,我们可以吧从系数表达法转换为点值表达法的过程用矩阵来表示:
事实上,我们会发现,把中间那个矩阵的所有元素全部取一个倒数,再乘上一个,就可以得到原本的矩阵。
这样,我们就可以把点值表达法转换回系数表达法。
递归版FFT
从上面我们可以很容易得到FFT的递归写法。
代码实现:
void init(){
int limit=1;
while(limit<=n+m)
limit<<=1;
}
void FFT(int limit,cmplx *a,int type){
if(limit==1)
return;
cmplx a1[limit>>1],a2[limit>>1];
for(int i=0;i<=limit;i+=2){
a1[i>>1]=a[i];
a2[i>>1]=a[i+1];
}
FFT(limit>>1,a1,type);
FFT(limit>>1,a2,type);
const double x=2.0*pai/(const int)limit;
cmplx wn=cmplx(cos(x),type*sin(x)),w=cmplx(1,0);
for(int i=0;i<(limit>>1);i++,w*=wn){
const cmplx y=w*a2[i];
a[i]=a1[i]+y;
a[i+(limit>>1)]=a1[i]-y;
}
return;
}
但是递归版的FFT的常数非常大,一不小心就会被卡飞掉。所以我们需要一个更加高效的做法。
迭代版FFT
在递归版FFT中,我们要将奇数项和偶数项分离。
以时为例,我们可以观察递归版的FFT每一次进行分离之后的序列:
0 1 2 3 4 5 6 7
0 2 4 6 1 3 5 7
0 4 2 6 1 5 3 7
0 4 2 6 1 5 3 7
再把它们转化为二进制看看:
000 100 010 110 001 101 011 111
再把它们的二进制反转后:
000 001 010 011 100 101 110 111
神奇的事情发生了!得到的这个序列的二进制反转后与原序列是一样的。
所以我们就可以通过这个很轻松地求出这个序列。
然后我们就可以迭代实现FFT了!
代码如下:
void init(){
int limit=1,l=0;
while(limit<=n+m){
limit<<=1;
l++;
}
for(int i=0;i<limit;i++)
r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
return;
}
void FFT(int limit,cmplx *a,int type){
for(int i=0;i<limit;i++)
if(i<r[i])
swap(a[i],a[r[i]]);
for(int mid=1;mid<limit;mid<<=1){
const double x=pai/(const int)mid;
cmplx wn=cmplx(cos(x),type*sin(x));
for(int R=mid<<1,j=0;j<limit;j+=R){
cmplx w=cmplx(1,0);
for(int k=0;k<mid;k++,w*=wn){
cmplx y=a[j+k],z=w*a[j+k+mid];
a[j+k]=y+z;
a[j+k+mid]=y-z;
}
}
}
return;
}
给大家几道题:
【BZOJ3160】万径人踪灭
【BZOJ3771】Triple
【BZOJ4827】礼物
如果有误在评论区吼一声哦!