问题
给定两个十进制数(105位),求它们的乘积。
不妨把它归为一个多项式乘积的问题:每一位都是一个系数,那么1234*2333
就变成了:
(x3+2x2+3x+4)∗(2x3+3x2+3x+3)
对于多项式乘积问题,显然跑O(n2)的朴素高精度乘法是过不了的。考虑我们计算a×b的方式:令X为答案,则有:
Xn=i=0∑nai×bn−i
因此,我们做的事情是直接计算答案的每一位。现在我们换一种思路。
点值表达
之前我们使用的(x5+233x3+x)这种表达方式,被称为系数表达
。因为它给出了系数向量:(1,0,233,0,1,0)。
但是考虑这样一个事实:**给定n个点,可以唯一确定一个n−1次多项式函数。**至于如何确定,有高斯消元和拉格朗日插值法。
因此我们拥有了一种全新的表达多项式的方式:点值表达。给出n+1个点,可以表达一个多项式。
例如:(0,0),(1,1),(2,4)是多项式(x2)的一种点值表达。一个多项式有无数组点值表达。
有什么用呢?
考虑多项式的加法A+B,生成的多项式的点值表达,可以由A和B的点值表达得到。在A和B上取相同的一些x,求出对应的点值表达:
A:(x1,ya1),(x2,ya2)⋯
B:(x1,yb1),(x2,yb2)⋯
则A+B:(x1,ya1+yb1),(x2,ya2+yb2)⋯
看图很好理解:
我们把从点值表达变成系数表达
的过程称作插值。
那么多项式的乘法也类似。大概是这样的步骤:
单位复数根
对着数学书脑补一番就解决了。
单位复数根:ωnk均匀分布在以(0,0)为中心的复平面圆上。n=8时长成这样:
ωnk=e2πink(注意其中的i是虚数单位)。欧拉告诉我们,eiu=cos(u)+isin(u),故有:ωnk=cos(2πk/n)+isin(2πk/n).
看图应该能脑补出来:尽管ωnk有n种取值,(ωnk)2只有$\frac{n}{2\dfrac。
在图中的意义是:第一象限点
的平方与第三象限点
的平方一致;第二象限点
的平方与第四象限点
的平方一致。因为(a+b)2=(−a−b)2。
由于这货的性质,我们选择单位复数根
作为x坐标进行求值。
FFT
关键步骤:
将A=i=0∑naixi拆分为:
A0=a0,a2x,a4x2,⋯,an−2x2n−1
A_1={a_1,a_3x,a_5x^2,…,a_{n-1}x^{\\dfracn}{2}-1}}
则有:A(x)=A0(x2)+xA1(x2)
由于我们使用单位复数根进行求值,则x2只有n/2种取值。我们把问题规模成功降低了一半!
那么如何插值回来呢?
**令ωnk=cos(2πk/n)−isin(2πk/n)**即可。(这里变成减号)
上面那段话并不清楚,因为我太弱了,根本讲不清楚。要完全理解上面的话,推荐看完代码,然后啃一啃《导论》。
所以我们就写出了代码:uoj #34.多项式乘法
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
| #include <cstdio> #include <cstdlib> #include <cmath> #include <complex> #include <iostream> using namespace std; typedef complex<double> cp;
void fft(cp *a,int n,int flag) { int i; cp a0[n/2+1],a1[n/2+1];
if(n==1) return; cp w_n(cos(2*M_PI/n),sin(flag*2*M_PI/n)); cp w(1,0);
for(i=0;i<n/2;i++) a0[i]=a[i*2],a1[i]=a[i*2+1];
fft(a0,n/2,flag); fft(a1,n/2,flag);
for(i=0;i<n/2;i++) { a[i]=a0[i]+w*a1[i]; a[i+n/2]=a0[i]-w*a1[i]; w=w*w_n; }
}
cp x[300005]={0},y[300005]={0}; int n,m;
void init() { int i; scanf("%d%d",&n,&m); for(i=0;i<=n;i++) cin>>x[i].real(); for(i=0;i<=m;i++) cin>>y[i].real(); m+=n; for(n=1;n<=m;n=n*2); }
int main(void) { int i; init(); fft(x,n,1); fft(y,n,1);
for(i=0;i<n;i++) x[i]=x[i]*y[i]; fft(x,n,-1);
for(i=0;i<=m;i++) printf("%d ",int((x[i].real())/n+0.5));
return 0; }
|
跑了4460ms。提交记录
递归版的比较慢(废话),迭代版的跑得比香港记者还快。
然而我并不会蝴蝶操作那套理论,请看riteme的博客:有关多项式的算法
相关资料
riteme 快速数论变换(NTT)
简要介绍了快速数论变换。
xlightgod UOJ34 多项式乘法
适合入门FFT。我就是在那里学习了一个。
iamzky 快速傅里叶变换
讲解很详细。