用途
- 让我们看一下下面这一个问题:
- 对于
A(x)=a0+a1x+a2x2+...+an−1xn−1,
B(x)=b0+b1x+b2x2+...+bm−1xm−1
- 求
C(x)=A(x)∗B(x)
暴力
- 我们显然可以枚举A的每一次项,B的每一次项。
- 然后一一乘起来,就是答案了!
-
ci=j=0∑j=iaj∗bi−j
- 但是这样子做的时间复杂度
O(n∗m)的,显然很慢!
- 所以,我们要想办法优化
点值与插值
- 在介绍更优的算法之前,就让我们先了解另一个多项式的表示方式
- 在通常情况下,我们是用“系数表示法”,也就是刚才的形如:
-
A(x)=a0+a1x+a2x2+...+an−1xn−1的形式
- 但其实,还有另一种表示法,叫做“点值表示法”。
- 对于一个最高次数为
n−1的多项式,我们可以选
n个不同的数
(x0,x1,...,xn−1)代入次多项式,得出
n个结果
(y0,y1,...,yn−1)
- 可以证明,只要给出这
2∗n个数,
x与
y,就可以得出有且仅有一个多项式
- 证明,考虑变成矩阵乘法:
- 只要证明“x矩阵”是可逆的,就可以用
a=y∗x−1来得出解。
- 我们把x矩阵变成范德蒙德行列式,可以通过证明得出右图的等式(数学归纳法)
- 我们知道这是不为0的,因为这些数两两不相同。
- 所以x矩阵是可逆的,a也就是唯一的
- 而我们把点值变回原多项式的操作叫做插值
瓶颈
- 既然我们只要知道了点值与插值,我们就可以把
A(x)与
B(x)用相同的
x带入,然后得出两组不同的
y,再把这两组
y一一乘起来,最后用插值求出
C(x)就行了。
- 但很可惜的是,朴素的点值和插值不管怎么样,似乎都是
O(n2)的,这样子还不如直接暴力呢!!!(说那么多有什么用?)
- 但不要就此放弃,就在这时,傅里叶站了出来,他找到了突破口。
突破口
- 他提出了一个关键的思想:如果把点值时取的
x特殊化,会不会有什么新发现呢?
- 顺着这个思想,他找到了n次单位复根。
n次单位复根
- 这个名词听起来很高深,实际也挺高深的。
- 说白了,就是如果
xn=1,它就是n次单位复根。
- 在实数数域内,可能符合条件的只有1和-1,而且-1得在偶数次方的情况下才成立。
- 但如果推广到复数,满足条件的就有n个数了。
- 首先来简单地讲一下复数是什么。
- 其实,实数已经包含了所有存在的数,但就在这时,人们发现又有一类数,它们虽然不存在,但在解决一些数学问题中有重要的作用,这类数的核心就是——
−1
,也就是
i。
- 在很多情况下,我们把复数数域看成是一个平面直角坐标系,x轴表示实数部分,y轴表示虚数部分。例如一个数
3+5−1
=3+5i,它在这个坐标系上对应的点就是
(3,5)。
- 复数作为一种数,自然也有乘法
-
(a+bi)∗(c+di)=(ac−bd)+(ad+bc)i,显然
- 但如果把它放在坐标系上,乘积也是有规律的。就是模长相乘,幅角相加。
- 模长指的是这个复数所在的点与原点的距离,而幅角就是这一条连线与x轴正半轴的夹角。
- 具体证明这里就不赘述了。
- 那么,对于一个复数,它的n次方的模长是原数的n次方,而幅角就是原数的n倍。
- 好的,下面是关键部分:
- 我们知道,如何一个非零数的0次方都为1,所以
(n1
)0=1(这里的次方根是在复数数域内的)
- 而
(n1
)n=1也是成立的。
- 也就是说在乘了
n次
n1
之后,模长没有变(1),幅角则是转了360度之后又回到了原位。
- 所以,
wn1(我们一般把
(x1
)y表示为
wxy)在平面坐标系上模长为1,幅角为
n2π。
- 那么,剩下的就是一些简单的三角函数了,这里就不赘述了,最终
wn1的坐标就是
(cos(n2π),sin(n2π))。
有关的定理
群的性质
- 因为
wn0=wnn=1
- 所以
wni∗wnj=wni+j=wn(i+j)%n
- 于是有
wn−1=wnn−1,
wnk=wnk+n,
wnk=wn2k
消去引理
-
wdndi=wni
折半引理
- 对于任意正整数n,2*n次单位根的平方的集合与n次单位根集合相同。
- 也就是
{(w2n)2}={wn}
求和引理
-
j=0∑n−1(wnk)j=0
- 证明,通过等比数列求和可得到原式等于
wnk−1(wnk)n−1
- 其分子=
wnnk−1=
1−1=
0
- 得证。
FFT
- 好了,铺垫了那么久,终于来到了最终的正题。
- 我们把这特殊的n次单位根带入,就会有特殊的结果。
- 这特殊的结果就是多项式A的离散傅里叶变换(DFT)
- 紧接着,我们考虑分治。
- 首先,我们把n补成形如2^x(x为正整数)的数,原本没有的项就把系数补为0。
- 然后,我们考虑这样的一个多项式
A(x):
-
a0+a1x+a2x2+...+an−1xn−1
- 我们把奇数次项和偶数次项分开:
-
(a0+a2x2+...+an−2xn−2)+(a1x+a3x3+...+an−1xn−1)
- 于是,我们把这个多项式分成了两个多项式:
-
A0(x2)+x∗A1(x2)
- 在这里,项数也就变为了
n/2,一直到变为1为止,这时候,只剩下1个常数项,就返回1个代入
w10(因为有多少项就要有多少个点值),也就是当前的这个系数。
- 假设我们已经处理好了
A0与
A1两个
wn/20~
wn/2n/2−1的点值,考虑如何合并。
- 对于代入
wnk,有以下两种情况:
-
-
k<2n,那么
A(wnk)=A0(wn/2k)+wnk∗A1(wn/2k),直接相加即可(非常简单)
- 2.否则,
A(wnk+n/2)=A0(wn/2k+n/2)+wnk+n/2∗A1(wn/2k+n/2)=A0(wn/2k)−wnk∗A1(wn/2k)
- 然后,就直接合并就好了(是不是比想象中的要简单呢?)
- 然后就这样合并
log2n次,就可以得到点值了!!!!!
- 这样点积的时间复杂度就是
O(nlog2n)
插值
- 点积已经被我们KO了,那么插积怎么搞呢?
- 实际上,DFT有一个非常巧妙的结论。
- 对于一个多项式A,我们求出它的DFT,然后把得出来的
(y0,y1,...,yn−1)当做是一个新的多项式B的系数,然后用n次单位复根的的倒数作为代入值再对B求一遍DFT,得出一个
z
- 可证
ai=nzi
- 证明,我们先把
zk用式子求出来
-
zk=i=0∑n−1yi∗(wn−k)i=i=0∑n−1(j=0∑n−1aj∗(wni)j)∗(wn−k)i
-
=i=0∑n−1j=0∑n−1aj∗(wn−k)i∗(wni)j=i=0∑n−1j=0∑n−1aj∗wn(j−k)∗i
-
=j=0∑n−1aj∗i=0∑n−1wn(j−k)∗i
- 若
j=k,
∑i=0n−1wn(j−k)∗i显然为
n,不然,就是
∑i=0n−1(wni)j−k,通过求和引理我们可以知道这个式子为0.
- 所以
ak=nzk,得证
- 这样子的话,就相当于把插值变成了点值,时间复杂度为
O(nlog2n)
一些关于实现的东西
- 直接把所有分治时点值存下来是十分不现实的。
- 所以我们最好开滚动数组,这样子合并起来就很方便,但注意复数的运算要打好点,不要让常数太大(毕竟是复数嘛)。
代码(多项式乘法)(洛谷3803 模板题)
#include<cstdio>
#include<cstring>
#include<cmath>
#define pi M_PI
using namespace std;
struct xushu{
double x,y;
xushu operator +(const xushu &p){
return {x+p.x,y+p.y};
}
xushu operator -(const xushu &p){
return {x-p.x,y-p.y};
}
xushu operator *(const xushu &p){
return {x*p.x-y*p.y,x*p.y+y*p.x};
}
xushu operator /(const int &p){
return {x/p,y/p};
}
};
xushu a[2100000][2],b[2100000];
int read(){
int x=0;char ch=getchar();
while(ch<'0'||ch>'9') ch=getchar();
while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
return x;
}
int mymax(int x,int y) {return x>y?x:y;}
int abs(int x) {return x>0?x:-x;}
inline xushu mi(int a,int b){
if(a<0) a+=b;
xushu t={cos(2*pi/b*a),sin(2*pi/b*a)};
return t;
}
int g[2100000],t;
void find(int len,int bk){
for(int i=0;i<len;i++){
if(g[i]>i)
{xushu t=a[i][0];a[i][0]=a[g[i]][0];a[g[i]][0]=t;}
}
int d=1;t=0;
while(1){
len/=2;d*=2;
if(!len) break;t=1-t;
for(int i=1;i<=len;i++){
xushu one=mi(bk,d),temp={1,0};
for(int j=0;j<d/2;j++){
int t2=(i-1)*d+j,t3=t2+d/2;
a[t2][t]=a[t2][1-t]+temp*a[t3][1-t];
a[t3][t]=a[t2][1-t]-temp*a[t3][1-t];
temp=temp*one;
}
}
}
}
xushu c[2100000],d[2100000];
xushu ans[2100000];
int main()
{
int n,m,xx;scanf("%d%d\n",&n,&m);n++;m++;
for(int i=0;i<n;i++) a[i][0].x=read();
for(int i=0;i<m;i++) b[i].x=read();
int len=n+m-1,d1=1,d2=0;
while(d1<len) d1*=2,d2++;len=d1;
for(int i=0;i<len;i++){
int e=i,h=0;
for(int j=1;j<=d2;j++) h=h*2+e%2,e/=2;
g[i]=h;
}
find(len,1);
for(int i=0;i<len;i++) c[i]=a[i][t],a[i][0]=b[i];
find(len,1);
for(int i=0;i<len;i++) d[i]=a[i][t];
for(int i=0;i<len;i++) c[i]=c[i]*d[i];
for(int i=0;i<len;i++) a[i][0]=c[i];
find(len,-1);
for(int i=0;i<len;i++) ans[i]=a[i][t]/len;
for(int i=0;i<n+m-1;i++) printf("%.0lf ",ans[i].x+0.001);
printf("\n");
return 0;
}