又是一篇咕咕咕了好久的文章。
什么是多项式
一个形如(a_0+a_1x_1+a_2x_2+...+a_nx^n)的东西叫做关于(x)的(n)次多项式。
其中(a_i)叫做这个多项式的(i)次项系数。特别地,(a_0)是这个多项式的常数项。
多项式乘法
一个(n)次多项式和一个(m)次多项式的乘积是一个(n+m)次多项式。
设第一个多项式的系数构成的数列为(a),另一个为(b)。
则这两个多项式的乘积的系数构成的数列(c)是(a)和(b)的卷积。
即:
(c_x=sum_{i=0}^xa_i imes b_{x-i})
显然(c)有(n+m+1)项。
如何实现多项式乘法
首先考虑暴力。
直接按照上述式子求出(c)即可。时间复杂度(O(nm))。
能不能更快呢?
点值表示法
可以证明,对于一个(n)次多项式(f(x)),我们只需要它在(n+1)个不同的(x)的取值时对应的(f(x))的之,就可以唯一确定这个多项式。
于是就有了点值表示法。
一个(n)次多项式的点值表示法形如一个二元组集合(P),(|P|=n+1)。
(P={(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),...,(x_n,f(x_n))})
有什么用?
如果我们知道(n)次多项式(f(x))和(m)次多项式(g(x))在互不相同的(x_0,x_1,x_2,...,x_{n+m})上的取值,我们就可以知道(f(x) imes g(x))的点值表示
(P_{f(x) imes g(x)}={(x_0,f(x_0) imes g(x_0)),(x_1,f(x_1) imes g(x_1)),(x_2,f(x_2) imes g(x_2)),...,(x_{n+m},f(x_{n+m}) imes g(x_{n+m}))})
也就是说,点值表示下的多项式乘法是(O(n))的!
于是我们兴奋地尝试去把这两个多项式转成点值表示。
我们枚举(1)到(n+m)作为选定的(x),对于每个(x)可以(O(n+m))求得(f(x),g(x))的值。
然后我们(O(n+m))把它乘起来。
于是我们发现这样做的复杂度是(O((n+m)^2))的。
比暴力还要慢QWQ
不过我们的直觉告诉我们,这种奇怪的方法虽然慢,但是有优化的余地。
于是有一个叫傅里叶的神仙,怀着这样的信念,于是就有了——
快速傅里叶变换(FFT)
我们之前的点值表示都局限于实数。
既然实数不行,我们就需要请出......
复数
复数被表示为(a+bi),其中(a,b)为实数,(i=sqrt{-1})。
我们把复数表示在复平面上,就像我们把实数表示在数轴上。
复平面上有横轴和纵轴,即实轴和虚轴。复平面上的点((x,y))表示复数(x+yi)。
(个人喜欢用数对的形式表示复数)
然后介绍一下复数的四则运算。
类似实数,把(i)当成未知数,注意(i^2=-1)即可。
也就是
((a,b)+(c,d)=(a+b,c+d))
((a,b)-(c,d)=(a-c,b-d))
((a,b) imes(c,d)=ac+adi+bci+bdi^2=(ac-bd,ad+bc))
(frac{(a,b)}{(c,d)}=frac{(a,b) imes(c,-d)}{(c,d) imes(c,-d)}=frac{(ac+bd,bc-ad)}{c^2+d^2}=(frac{ac+bd}{c^2+d^2},frac{bc-ad}{c^2+d^2}))
哦对了,对于复数((a,b)),((a,-b))是它的共轭复数。复数除法的过程相当于将分子分母同乘分母的共轭复数。
我们来重点看一看复数乘法。
复数乘法有一个很有用的性质。如果我们将复数((a,b))对应向量((a,b)),那么复数的乘法可以表达为:
模长相乘,幅角相加。
模长:对于向量((a,b)),其模长为点((a,b))与原点((0,0))的距离;幅角:对于向量((a,b)),其幅角为它与横轴正半轴的夹角。
即对于复数乘法((a,b) imes(c,d)),其结果所对应的向量的模长为向量((a,b))的模长乘向量((c,d))的模长,幅角为向量((a,b))的幅角加向量((c,d))的幅角。
有啥用?
先别急,我们在复数下定义一个东西。
单位复根
满足((a,b)^n=1)的复数((a,b))叫(n)次单位复根。
那么它的分布如何呢?
首先,单位复根的模长必然为1。因为当模长小于1,它乘(n)次以后模长还是小于1;大于1同理。
然后,我们发现如果复数((a,b))是一个单位复根,那么它的幅角( heta)一定满足
(frac{n heta}{2pi}in Z)
不难理解吧,也就是一个向量从((1,0))开始,每次旋转( heta)的角度,旋转(n)次以后还落在((1,0))上,那么它一定旋转了整数圈。
于是我们发现(n)次单位根正好是模长为1,幅角为(frac{2kpi}{n})的向量对应的复数。
我们记(omega_n^k)表示模长为1,幅角为(frac{2kpi}{n})的向量对应的(n)次单位复根,称为第(k)个(n)次单位复根。
我们还能发现,(omega_n^k=omega_n^{k mod n}),所以一般情况下,我们认为(n)次单位复根有(n)个,即(omega_n^0,omega_n^1,omega_n^2,...,omega_n^{n-1})。
然后我们发现它们有一些美妙的性质。
1.(n)次单位复根对应的向量(n)等分单位圆。
显然吧。因为两个相邻(n)次单位复根(认为(omega_n^{n-1})与(omega_n^0)相邻)对应的向量的夹角是相等的。
2.(omega_n^k=omega_n^{k mod n})
这个之前有提到,因为在弧度制下,任意弧度( heta)与( heta+2kpi(kin Z))表示一个相同的角。
3.((omega_n^k)^p=omega_n^{kp})
第(k)个(n)次单位复根对应向量的幅角变为原来的(p)倍,自然是对应(omega_n^{kp})对应的向量啦。
4.(omega_n^{2k}=omega_{frac{n}{2}}^k)
前提是(n)是偶数。
不难理解,(n)次单位复根只是在(frac{n}{2})次单位复根的基础上增加了(k)为奇数的(omega_n^{k}),(k)为偶数的部分的集合刚好是(frac{n}{2})次单位复根。
5.(omega_n^{k+frac{n}{2}}=-omega_n^k)
前提是(n)是偶数。
相当于一个复数对应的向量进行一次中心对称,复数((a,b))变为((-a,-b))。
DFT
终于开始讲FFT了。
不知是不是收到了上天的旨意,傅里叶突发奇想将(n)次单位复根带入了(n-1)次多项式(f(x)),希望得到意想不到的效果。
假定(n)是偶数。
设(f(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1})
我们将这(n)个单项式按照次数的奇偶性分类。
(f(x)=(a_0+a_2x^2+a_4x^4+...+a_{n-2}x^{n-2})+(a_1x^1+a_3x^3+a_5x^5+...+a_{n-1}x^{n-1}))
(f(x)=(a_0+a_2x^2+a_4x^4+...+a_{n-2}x^{n-2})+x(a_1+a_3x^2+a_5x^4+...+a_{n-1}x^{n-2}))
再将(omega_n^k)带入
(f(omega_n^k)=(a_0+a_2(omega_n^k)^2+a_4(omega_n^k)^4+...+a_{n-2}(omega_n^k)^{n-2})+omega_n^k(a_1+a_3(omega_n^k)^2+a_5(omega_n^k)^4+...+a_{n-1}(omega_n^k)^{n-2}) (*))
(f(omega_n^k)=(a_0+a_2omega_n^{2k}+a_4(omega_n^{2k})^2+...+a_{n-2}(omega_n^{2k})^{frac{n-2}{2}})+omega_n^k(a_1+a_3omega_n^{2k}+a_5(omega_n^{2k})^2+...+a_{n-1}(omega_n^{2k})^{frac{n-2}{2}}))
(f(omega_n^k)=(a_0+a_2omega_frac{n}{2}^k+a_4(omega_frac{n}{2}^k)^2+...+a_{n-2}(omega_frac{n}{2}^k)^{frac{n-2}{2}})+omega_n^k(a_1+a_3omega_frac{n}{2}^k+a_5(omega_frac{n}{2}^k)^2+...+a_{n-1}(omega_frac{n}{2}^k)^{frac{n-2}{2}}))
发现左右两边明显是两个递归子问题。
但是我们发现(frac{n}{2})次单位复根只有(frac{n}{2})个,所以上面的式子只适用于(k<frac{n}{2})的情况。
但是我们需要(n)个值来确定这个多项式。
那么怎么求(omega_n^{k+frac{n}{2}}(0le k<frac{n}{2}))呢?
我们知道(omega_n^{k+frac{n}{2}}=-omega_n^k)
然后((*))式中括号内每一项的次数都是偶数,所以符号不变,只有括号外的(omega_n^k)变为了(-omega_n^k)。
所以
(f(omega_n^{k+frac{n}{2}})=-f(omega_n^k)=(a_0+a_2omega_frac{n}{2}^k+a_4(omega_frac{n}{2}^k)^2+...+a_{n-2}(omega_frac{n}{2}^k)^{frac{n-2}{2}})-omega_n^k(a_1+a_3omega_frac{n}{2}^k+a_5(omega_frac{n}{2}^k)^2+...+a_{n-1}(omega_frac{n}{2}^k)^{frac{n-2}{2}}))
然后就是实现。
等一下!
我们之前一直要求(n)是偶数。
那万一不是呢?
其实很简单,我们把原来的(n)增加到大于等于(n)的2的整数次幂即可。
我们定义(FFT(l,r))表示对区间(a_{l,r})进行(DFT)。设(r-l+1=len),(ln=frac{len}{2})。
当(len=1),什么都不用做直接返回即可。
首先,奇偶分类。那么我们先将(a_{l,r})复制到(cpy_{l,r}),然后从(l)到(l+ln-1)枚举(i),计数器(tmp)一开始等于(l),每次令(a_i=cpy_{tmp},a_{i+ln}=cpy_{tmp+1}),并让(tmp)增加2。
然后,递归到(FFT(l,l+ln-1),FFT(l+ln,r))。
之后定义(rt)为第1个(len)次单位复根,即(omega_n^1),由它所对应的向量幅角为(frac{2pi}{n}),模长为1得到(omega_n^1=(cos(frac{2pi}{n}),sin(frac{2pi}{n})))。
定义当前(n)次单位复根为(nw),一开始为(omega_n^0),即1。
然后将已经经过递归变换的(a_{l,r})复制到(cpy_{l,r}),从(l)到(l+ln-1)枚举(i),令(a_i=cpy_i+nw imes cpy_{i+ln},a_{i+ln}=cpy_i-nw imes cpy_{i+ln})
结束。
时间复杂度(O(nlog n))
代码会和IDFT一起放出。
IDFT
有了DFT,我们就可以轻松地将点值表示的多项式相乘。
可我们总不能输出点值表示的多项式吧。
我们需要将它转回系数表达的多项式。
其实,IDFT的过程,就是把DFT中的(omega_n^1)换成了它的共轭复数,即((cos(frac{2pi}{n}),-sin(frac{2pi}{n}))),得到的系数再除以(frac{1}{n})即可。
证明比较长,略。
其实我自己也是一知半解
那么终于到了放代码的时刻。
code:
const double pi=acos(-1);//pi
struct com{//复数结构体
double re,in;
}f[300010],f1[300010],f2[300010],cpy[300010];
com operator+(const com &x,const com &y){
return (com){x.re+y.re,x.in+y.in};
}
com operator-(const com &x,const com &y){
return (com){x.re-y.re,x.in-y.in};
}
com operator*(const com &x,const com &y){
return (com){x.re*y.re-x.in*y.in,x.re*y.in+x.in*y.re};
}
void FFT(com *f,int l,int len,int op){
if(len<2)return;
int nl=len>>1,tmp=l-1;
for(int i=l;i<l+len;++i)cpy[i]=f[i];//复制原f数组到cpy
for(int i=l;i<l+nl;++i)f[i]=cpy[++tmp],f[i+nl]=cpy[++tmp];//按照奇偶分类
FFT(f,l,nl,op),FFT(f,l+nl,nl,op);//递归处理
com w=(com){cos(pi/(1.0*nl)),sin(pi/(1.0*nl)*op)},nw=(com){1,0},t;//w是第一个单位复根,nw是当前单位复根
for(int i=l;i<l+nl;++i,nw=nw*w)t=nw*f[i+nl],cpy[i]=f[i]+t,cpy[i+nl]=f[i]-t;
for(int i=l;i<l+len;++i)f[i]=cpy[i];
}
例题
1.洛谷模板
code:
#include<bits/stdc++.h>
using namespace std;
const double pi=2*acos(-1);
struct comp{
double re,in;
}tmp[2100010],rt,pw,t,a[2100010],b[2100010];
comp operator+(const comp &x,const comp &y){
return (comp){x.re+y.re,x.in+y.in};
}
comp operator-(const comp &x,const comp &y){
return (comp){x.re-y.re,x.in-y.in};
}
comp operator*(const comp &x,const comp &y){
return (comp){x.re*y.re-x.in*y.in,x.re*y.in+x.in*y.re};
}
int n,m,w,tg,mid,nw,sz;
void Work(comp *p,const int &l,const int &r){
for(int i=l;i<=r;i++)tmp[i]=p[i];
nw=l,sz=l-1;
while(nw<=r)p[++sz]=tmp[nw],nw+=2;
nw=l+1;
while(nw<=r)p[++sz]=tmp[nw],nw+=2;
}
void FFT(comp *p,const int &l,const int &nl,const int &op){
if(nl<2)return;
Work(p,l,l+nl-1);
FFT(p,l,nl>>1,op),FFT(p,l+(nl>>1),nl>>1,op);
rt=(comp){cos(pi/nl),sin(pi/nl)*op},pw=(comp){1,0};
mid=(l+l+nl-1)/2;
for(int i=l;i<l+nl;i++)tmp[i]=p[i];
for(int i=l;i<=mid;i++,pw=pw*rt)t=pw*tmp[mid+i-l+1],p[i]=tmp[i]+t,p[mid+i-l+1]=tmp[i]-t;
}
void scan(double &x){
x=0;
char c=getchar();
while('0'>c||c>'9')c=getchar();
while('0'<=c&&c<='9')x=x*10+c-'0',c=getchar();
}
void print(int x){
x>9?print(x/10),0:0,putchar(x%10+'0');
}
int main(){
scanf("%d%d",&n,&m);
n++,m++;
for(int i=0;i<n;i++)scan(a[i].re);
for(int i=0;i<m;i++)scan(b[i].re);
n+=m-1,m=1;
while(m<n)m<<=1;
FFT(a,0,m,1),FFT(b,0,m,1);
for(int i=0;i<m;i++)a[i]=a[i]*b[i];
FFT(a,0,m,-1);
for(int i=0;i<n;i++)print((int)(fabs(a[i].re)/m+0.1)),putchar(' ');
return 0;
}
2.高精度乘法洛谷模板
code:
#include<bits/stdc++.h>
using namespace std;
const double pi=2*acos(-1);
struct comp{
double re,in;
}tmp[150010],rt,pw,t,a[150010],b[150010];
comp operator+(const comp &x,const comp &y){
return (comp){x.re+y.re,x.in+y.in};
}
comp operator-(const comp &x,const comp &y){
return (comp){x.re-y.re,x.in-y.in};
}
comp operator*(const comp &x,const comp &y){
return (comp){x.re*y.re-x.in*y.in,x.re*y.in+x.in*y.re};
}
int n,m,w,tg,mid,nw,sz;
void Work(comp *p,const int &l,const int &r){
for(int i=l;i<=r;i++)tmp[i]=p[i];
nw=l,sz=l-1;
while(nw<=r)p[++sz]=tmp[nw],nw+=2;
nw=l+1;
while(nw<=r)p[++sz]=tmp[nw],nw+=2;
}
void FFT(comp *p,const int &l,const int &nl,const int &op){
if(nl<2)return;
Work(p,l,l+nl-1);
FFT(p,l,nl>>1,op),FFT(p,l+(nl>>1),nl>>1,op);
rt=(comp){cos(pi/nl),sin(pi/nl)*op},pw=(comp){1,0};
mid=(l+l+nl-1)/2;
for(int i=l;i<l+nl;i++)tmp[i]=p[i];
for(int i=l;i<=mid;i++,pw=pw*rt)t=pw*tmp[mid+i-l+1],p[i]=tmp[i]+t,p[mid+i-l+1]=tmp[i]-t;
}
void scan(double &x){
x=0;
char c=getchar();
while('0'>c||c>'9')c=getchar();
x=c-'0';
}
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++)scan(a[n-i].re);
for(int i=1;i<=n;i++)scan(b[n-i].re);
n=n*2-1,m=1;
while(m<n)m<<=1;
FFT(a,0,m,1),FFT(b,0,m,1);
for(int i=0;i<m;i++)a[i]=a[i]*b[i];
FFT(a,0,m,-1);
for(int i=0;i<n;i++)a[i].re=fabs(a[i].re)/m,a[i].re=(int)(w+a[i].re+0.1),w=(int)a[i].re/10,a[i].re-=10*w,(w&&i==n-1)?n++:0;
for(int i=n-1;i>=0;i--)tg|=a[i].re>0,tg?putchar((int)a[i].re+'0'):0;
return 0;
}