点值表示法
我们正常表示一个多项式的方式,形如 A ( x ) = a 0 + a 1 x + a 2 x 2 + . . . + a n x n A(x)=a_0+a_1x+a_2x^2+...+a_nx^n A(x)=a0+a1x+a2x2+...+anxn,这是正常人容易看懂的,但是,我们还有一种表示法。
我们知道, n + 1 n+1 n+1个点可以表示出一个 n n n 次的多项式。
于是,我们任意地取 n + 1 n+1 n+1 个不同的值,表示 x x x ,求出的值与 x x x 对应,形成 n + 1 n+1 n+1个点,这就可以表示。
复数
一种表示坐标的方法,对于坐标 ( x , y ) (x,y) (x,y),可写作 x + i y x+iy x+iy,其中 x x x为实部, y y y为虚部。
C++中有复数的模板,complex
,可以直接作为变量类型使用。
运算规则,自然不用多说了,也就是直接拿式子算即可。
如果你不会,可以看看百度百科。
我们将点至原点的距离称为模长,将其与原点相连之后与 x x x 轴形成的一个夹角称为辐角,不过呢,对于第三第四象限的点,自然要加上一个 180 ° 180° 180°了。
我的表述自然不够专业,希望可以表述出这个意思吧。
复数的乘法可以理解为,模长相乘,辐角相加。
Tips:
此部分的证明来自 cjx 犇犇。
为啥是这样呢?证明如下:
( a + b i ) ( c + d i ) = a c − b d + i ( b c + a d ) (a+bi)(c+di)=ac-bd+i(bc+ad) (a+bi)(c+di)=ac−bd+i(bc+ad)
那么这个点就是 ( a c − b d , b c + a d ) (ac-bd,bc+ad) (ac−bd,bc+ad),其模长:
( a c − b d ) 2 + ( b c + a d ) 2 = ( a c ) 2 − 2 a b c d + ( b d ) 2 + ( b c ) 2 + 2 a b c d + ( a d ) 2 = ( a c ) 2 + ( b d ) 2 + ( b c ) 2 + ( a d ) 2 = ( a 2 + b 2 ) ( c 2 + d 2 ) \sqrt{(ac-bd)^2+(bc+ad)^2}\\ =\sqrt{(ac)^2-2abcd+(bd)^2+(bc)^2+2abcd+(ad)^2}\\ =\sqrt{(ac)^2+(bd)^2+(bc)^2+(ad)^2}\\ =\sqrt{(a^2+b^2)(c^2+d^2)}\\ (ac−bd)2+(bc+ad)2=(ac)2−2abcd+(bd)2+(bc)2+2abcd+(ad)2=(ac)2+(bd)2+(bc)2+(ad)2=(a2+b2)(c2+d2)
那么我们应该可以看出来这个模长相乘了。
接下来是辐角相加,我们设原来两个辐角为 θ 1 , θ 2 \theta_1,\theta_2 θ1,θ2,而模长为 t 1 , t 2 t_1,t_2 t1,t2。
cos ( θ 1 + θ 2 ) = cos ( θ 1 ) cos ( θ 2 ) − sin ( θ 1 ) sin ( θ 2 ) = a t 1 × c t 2 − b t 1 × d t 2 = a c − b d t 1 t 2 \cos(\theta_1+\theta_2)=\cos(\theta_1)\cos(\theta_2)-\sin(\theta_1)\sin(\theta_2)\\ =\frac {a}{t_1}\times\frac {c}{t_2}-\frac {b}{t_1}\times\frac {d}{t_2}\\ =\frac {ac-bd}{t_1t_2} cos(θ1+θ2)=cos(θ1)cos(θ2)−sin(θ1)sin(θ2)=t1a×t2c−t1b×t2d=t1t2ac−bd
我们知道分母是乘积的模长,分子是横坐标,所以这个式子恰好就是乘积辐角对应的 cos \cos cos 值。
那么显然,辐角的值就是 θ 1 + θ 2 \theta_1+\theta_2 θ1+θ2了。
把圆均分
如图是坐标轴上一个以原点为圆心的半径为 1 1 1 的圆。
我们定义 ω n k \omega_n^k ωnk表示从 ( 1 , 0 ) (1,0) (1,0)开始,把圆均分为 n n n份的第 k k k个点的复数,其中 ( 1 , 0 ) (1,0) (1,0)为 ω n 0 \omega_n^0 ωn0。
那么,不难发现, ω n k \omega_n^k ωnk表示的点为 ( cos ( 2 π k n ) , sin ( 2 π n k ) ) (\cos(2\pi\frac k n),\sin(2\pi\frac n k)) (cos(2πnk),sin(2πkn))。
这是为什么呢?我们考虑它的辐角,由于其平分了一整个圆,所以其辐角为 360 k n ° 360\frac k n° 360nk°,转换为弧度后则为 2 π k n 2\pi\frac k n 2πnk,且模长为 1 1 1,利用三角函数易得其坐标了。
简单推推式子不难发现 ( ω n 1 ) k = ω n k (\omega_n^1)^k=\omega_n^k (ωn1)k=ωnk,这是利用了模长相乘,辐角相加,因为模长是 1 1 1 ,怎么乘都是 1 1 1,于是辐角不断叠加,从定义上看是这样的。
- 性质1: ω d n d k = ω n k \omega _{dn}^{dk}=\omega_n^k ωdndk=ωnk,根据定义可证。
- 性质2: ω n k + n 2 = − ω n k \omega _{n}^{k+\frac n 2}=-\omega_n^k ωnk+2n=−ωnk,两点对称。
这个东西有什么用呢?
离散傅里叶变换
离散傅里叶变换(Discrete Fourier Transform,简称DFT)的思想是利用 ω n k \omega_n^k ωnk将一个多项式转为点值表示法。
对于一个多项式 A ( x ) = a 0 + a 1 x + a 2 x 2 + . . . + a n − 1 x n − 1 A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1} A(x)=a0+a1x+a2x2+...+an−1xn−1,我们按照前文所云,将所有的 ω n k \omega_n^k ωnk作为 x x x 代入。
于是我们得到了 n − 1 n-1 n−1 个点,使用复数形式表示,成为一个数组 ( y 0 , y 1 , y 2 , . . . , y n − 1 ) (y_0,y_1,y_2,...,y_{n-1}) (y0,y1,y2,...,yn−1)的。
这被称为 A ( x ) A(x) A(x) 的傅里叶变换。
傅里叶逆变换
我们再将其作为一个多项式 B ( x ) = y 0 + y 1 x + . . . . + y n − 1 x n − 1 B(x)=y_0+y_1x+....+y_{n-1}x^{n-1} B(x)=y0+y1x+....+yn−1xn−1。
对于这个多项式 B ( x ) B(x) B(x),代入所有的 ω n − k \omega_n^{-k} ωn−k,也就是 ω n k \omega_n^k ωnk的倒数,得到 ( z 0 , z 1 , z 2 , . . . , z n − 1 ) (z_0,z_1,z_2,...,z_{n-1}) (z0,z1,z2,...,zn−1)。
易得:
z k = ∑ i = 0 n − 1 y i ( ω n − k ) i = ∑ i = 0 n − 1 ( ∑ j = 0 n − 1 a j ( ω n i ) j ) ( ω n − k ) i = ∑ j = 0 n − 1 a j ( ω n i ) j ∑ i = 0 n − 1 ( ω n − k ) i = ∑ j = 0 n − 1 a j ( ω n j ) i ∑ i = 0 n − 1 ( ω n − k ) i = ∑ j = 0 n − 1 a j ∑ i = 0 n − 1 ( ω n j − k ) i z_k=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i \\ =\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^i)^j)(\omega_n^{-k})^i\\ =\sum_{j=0}^{n-1}a_j(\omega_n^i)^j\sum_{i=0}^{n-1}(\omega_n^{-k})^i\\ =\sum_{j=0}^{n-1}a_j(\omega_n^j)^i\sum_{i=0}^{n-1}(\omega_n^{-k})^i\\ =\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j-k})^i\\ zk=i=0∑n−1yi(ωn−k)i=i=0∑n−1(j=0∑n−1aj(ωni)j)(ωn−k)i=j=0∑n−1aj(ωni)ji=0∑n−1(ωn−k)i=j=0∑n−1aj(ωnj)ii=0∑n−1(ωn−k)i=j=0∑n−1aji=0∑n−1(ωnj−k)i
关于后面那个等比数列,若 j = k j=k j=k,可得 1 1 1,否则用等比数列式子可知为 0 0 0。
因此:
z k = n a k a k = z k n z_k=na_k\\ a_k=\frac {z_k} n zk=nakak=nzk
所以我们可以求出原来的多项式了。
有几点需要注意:
- 我们提取 a i a_i ai 只提取实部,因为虚部是虚数,无法经过我们的转换后变成实数。
- 数字有一定误差(毕竟你使用了三角函数等东西,小数是会有误差的),所以要四舍五入。
for(int i=0;i<=len-1;i++)
{
ans[i]+=floor(a[i].real()/len+0.5);
}
快速傅里叶变换
快速傅里叶变换(Fast Fourier Transform,简称FFT),是在 DFT的基础上我们发现时间复杂度依然需要 O ( n 2 ) O(n^2) O(n2),没有含金量,所以我们要给他含金量!
我们可以使用分治的思想,使得时间复杂度降至 O ( n log n ) O(n\log n) O(nlogn)。
对于 A ( x ) = a 0 + a 1 x + a 2 x 2 + . . . + a n − 1 x n − 1 A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1} A(x)=a0+a1x+a2x2+...+an−1xn−1,我们可以:
A 1 ( x ) = a 0 + a 2 x + . . . , A 2 ( x ) = a 1 + a 3 x + . . . A ( x ) = A 1 ( x 2 ) + x A 2 ( x 2 ) A ( ω n k ) = A 1 ( ω n 2 k ) + ω n k A 2 ( ω n 2 k ) A ( ω n k ) = A 1 ( ω n 2 k ) + ω n k A 2 ( ω n 2 k ) A_1(x)=a_0+a_2x+...,A_2(x)=a_1+a_3x+...\\ A(x)=A_1(x^2)+xA_2(x^2)\\ A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})\\ A(\omega_n^k)=A_1(\omega_{\frac n 2}^{k})+\omega_n^kA_2(\omega_{\frac n 2}^{k})\\ A1(x)=a0+a2x+...,A2(x)=a1+a3x+...A(x)=A1(x2)+xA2(x2)A(ωnk)=A1(ωn2k)+ωnkA2(ωn2k)A(ωnk)=A1(ω2nk)+ωnkA2(ω2nk)
同理:
A ( ω n k + n 2 ) = A 1 ( ω n 2 k + n ) + ω n k + n 2 A 2 ( ω n 2 k + n ) A ( ω n k + n 2 ) = A 1 ( ω n 2 k ) − ω n k A 2 ( ω n 2 k ) A(\omega_n^{k+\frac n 2})=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac n 2}A_2(\omega_n^{2k+n})\\ A(\omega_n^{k+\frac n 2})=A_1(\omega_{\frac n 2}^{k})-\omega_n^{k}A_2(\omega_{\frac n 2}^{k})\\ A(ωnk+2n)=A1(ωn2k+n)+ωnk+2nA2(ωn2k+n)A(ωnk+2n)=A1(ω2nk)−ωnkA2(ω2nk)
不错!利用这两个式子,我们可以在 O ( n log n ) O(n \log n) O(nlogn) 的时间复杂度求出 A ( x ) A(x) A(x)。
不过注意这里一直有一个除二操作,为了方便,我们需要把多项式补成一个次数为 2 x − 1 2^x-1 2x−1 的多项式。
可以写一个递归来求解。
注意这个取倒可以利用共轭复数,对于 a + i b a+ib a+ib,其共轭复数为 a − i b a-ib a−ib。
1 a + i b = a − i b ( a + i b ) ( a − i b ) = a − i b a 2 + b 2 \frac {1}{a+ib}=\frac {a-ib}{(a+ib)(a-ib)}=\frac {a-ib}{a^2+b^2} a+ib1=(a+ib)(a−ib)a−ib=a2+b2a−ib
由于此处 a 2 + b 2 = 1 a^2+b^2=1 a2+b2=1 ,因此其共轭复数为其倒数。
void FFT(cp *a,LL n,bool inv)//inv 表示omega是否取倒
{
if(n==1)return;
static cp buf[N];
LL m=n/2;
for(int i=0;i<=m-1;i++)buf[i]=a[2*i],buf[i+m]=a[2*i+1];//奇偶分开
for(int i=0;i<=n-1;i++)a[i]=buf[i];
FFT(a,m,inv),FFT(a+m,m,inv);
for(int i=0;i<=m-1;i++)
{
cp x=omega(n,i);
if(inv)x=conj(x);//conj(x)求解共轭复数
buf[i]=a[i]+x*a[i+m],buf[i+m]=a[i]-x*a[i+m];
}
for(int i=0;i<=n-1;i++)a[i]=buf[i];
}
利用FFT求解多项式的乘积
这个还是十分简单的,直接把两个多项式转化为两个长度相同的次数为 2 x − 1 2^x-1 2x−1的多项式。
求解出其傅里叶变换形式之后,对于每个点对应的复数,相乘即可。
为什么呢?
这里有一个误区大家要注意, a i a_i ai 不是一个点,而是单纯因为 ω n i \omega_n^i ωni是一个复数,所以 a i a_i ai 才变成了一个复数, a i a_i ai只是表示 A ( ω n i ) A(\omega_n^i) A(ωni)的值。
首先我们知道当前的 a i a_i ai表示的是 A ( ω n i ) A(\omega_n^i) A(ωni), b i b_i bi表示的是 B ( ω n i ) B(\omega_n^i) B(ωni),那么我们将两个值直接相乘。
因此多项式相乘以后,我们希望 a i = A ( ω n i ) B ( ω n i ) a_i=A(\omega_n^i)B(\omega_n^i) ai=A(ωni)B(ωni),那么就是直接相乘喽。
给一个简单的实现:
nclude<bits/stdc++.h>
#define LL int
#define cp complex<double>
using namespace std;
const double PI=acos(-1.0000);
const int N=5e6+5;
cp omega(LL n, LL k)
{
return cp(cos(2*PI*k/n),sin(2*PI*k/n));
}
LL n,x,len,ans[N];
cp a[N],b[N];
//上文的 FFT 实现省去
int main()
{
scanf("%d",&n);
len=1;
while(len<2*n)len*=2;
for(int i=n-1;i>=0;i--)scanf("%1d",&x),a[i].real(x);
for(int i=n-1;i>=0;i--)scanf("%1d",&x),b[i].real(x);
FFT(a,len,0),FFT(b,len,0);
for(int i=0;i<=len-1;i++)a[i]*=b[i];
FFT(a,len,1);
for(int i=0;i<=len-1;i++)//进位
{
ans[i]+=floor(a[i].real()/len+0.5);
ans[i+1]+=ans[i]/10;
ans[i]%=10;
}
int i=len;
for(i;i>=0&&ans[i]==0;i--);//前导零
if(i==-1)len=0;
for(;i>=0;i--)printf("%d",ans[i]);
}
非递归FFT
这里有一个优化,我们发现每次递归有一个把 a i a_i ai 奇偶分开的过程,本质来看,就是将二进制末尾为 0 0 0 的数字与二进制末尾为 1 1 1 的数字分开。
我们不妨想一下,对于一个数
x
x
x,其位置可以根据其二进制确定,就是其二进制倒过来的数字。
我们先将每个 a i a_i ai 放置在对应的位置,然后向上逐渐合并。
void FFT(cp *a,bool inv)
{
LL lim=0;
while((1<<lim)<len)lim++;
for(int i=0;i<=len-1;i++)
{
LL t=0;
for(int j=0;j<lim;j++)
if((i>>j)&1)t|=(1<<(lim-j-1));//处理其翻转后的值
if(i<t)swap(a[i],a[t]);
}
static cp buf[N];
for(int l=2;l<=len;l*=2)
{
LL m=l/2;
for(LL j=0;j<=len-1;j+=l)
{
for(LL i=0;i<=m-1;i++)
{
cp x=omega(l,i+j);
if(inv)x=conj(x);
buf[i+j]=a[i+j]+x*a[i+j+m];
buf[i+j+m]=a[i+j]-x*a[i+j+m];
}
}
for(int j=0;j<=len-1;j++)a[j]=buf[j];
}
}
蝴蝶操作
这个东西其实就是想了个办法使得把工具人数组 buf
除掉了。
调调顺序即可。
void FFT(cp *a,bool inv)
{
LL lim=0;
while((1<<lim)<len)lim++;
for(int i=0;i<=len-1;i++)
{
LL t=0;
for(int j=0;j<lim;j++)
if((i>>j)&1)t|=(1<<(lim-j-1));
if(i<t)swap(a[i],a[t]);
}
for(int l=2;l<=len;l*=2)
{
LL m=l/2;
for(LL j=0;j<=len-1;j+=l)
{
for(LL i=0;i<=m-1;i++)
{
cp x=omega(l,i+j);
if(inv)x=conj(x);
x*=a[i+j+m];
a[i+j+m]=a[i+j]-x;
a[i+j]=a[i+j]+x;
}
}
}
}
一些小优化
对于 i i i 的二进制翻转可以先预处理出来,我们用 R i R_i Ri 表示 i i i 的二进制翻转。
然后 ω n k \omega_n^k ωnk 可以利用性质累乘,会快很多,因为复数的各种操作常数巨大,直接累乘常数较小。
最后代码就长这样了:
#include<bits/stdc++.h>
#define LL int
#define cp complex<double>
using namespace std;
const double PI=acos(-1.0000);
const int N=5e6+5;
cp omega(LL n, LL k)
{
return cp(cos(2*PI*k/n),sin(2*PI*k/n));
}
LL n,len,lim,x,ans[N],r[N];
cp a[N],b[N];
void FFT(cp *a,bool inv)
{
for(int i=0;i<=len-1;i++)
{
LL t=r[i];
if(i<t)swap(a[i],a[t]);
}
for(int l=2;l<=len;l*=2)
{
LL m=l/2;
cp omg=omega(l,1);
for(LL j=0;j<=len-1;j+=l)
{
cp x(1,0);
for(LL i=0;i<=m-1;i++)
{
cp t=x;
if(inv)t=conj(t);
t*=a[i+j+m];
a[i+j+m]=a[i+j]-t,a[i+j]=a[i+j]+t;
x*=omg;
}
}
}
}
int main()
{
scanf("%d",&n);
len=1;
while(len<2*n)len*=2,lim++;
for(int i=0;i<=len-1;i++)
{
LL t=0;
for(int j=0;j<lim;j++)if((i>>j)&1)t|=(1<<(lim-j-1));
r[i]=t;
}
for(int i=n-1;i>=0;i--)scanf("%1d",&x),a[i].real(x);
for(int i=n-1;i>=0;i--)scanf("%1d",&x),b[i].real(x);
FFT(a,0),FFT(b,0);
for(int i=0;i<=len-1;i++)a[i]*=b[i];
FFT(a,1);
for(int i=0;i<=len-1;i++)
{
ans[i]+=floor(a[i].real()/len+0.5);
ans[i+1]+=ans[i]/10;
ans[i]%=10;
}
int i=len;
for(i;i>=0&&ans[i]==0;i--);
if(i==-1)len=0;
for(;i>=0;i--)printf("%d",ans[i]);
}
参考
胡小兔-小学生都能看懂的FFT!!!