2179: FFT快速傅立叶
Time Limit: 10 Sec Memory Limit: 259 MB
Submit: 1595 Solved: 792
[Submit][Status][Discuss]
Description
给出两个n位10进制整数x和y,你需要计算x*y。
Input
第一行一个正整数n。 第二行描述一个位数为n的正整数x。 第三行描述一个位数为n的正整数y。
Output
输出一行,即x*y的结果。
Sample Input
1
3
4
Sample Output
12
数据范围:
n<=60000
FFT模板题。
直接进行高精度乘法是
O(n2)
的,于是我们采用FFT来
O(nlogn)
实现:
(注明:以下着重介绍算法流程和算法思想,具体细节参考《算法导论》)
1.我们把乘数的每一位看作多项式的系数,得到多项式 A(x) (因为高精度乘法的本质就是多项式乘法)
2.首先求出
A(ωkn)
,其中
k∈[0,n−1]
,
ωn
是n次单位复根。
由于n次单位复根的一些奇妙性质:
相消引理
折半引理
我们可以采用分治
O(nlogn)
的时间求出这
n
项的值,但是递归实现常数较大,我们采用蝴蝶算法来迭代实现。
如图,把原来顺次排列的数列变成叶子中的顺序就可以迭代了~
(叶子中的顺序就是原序列的二进制逆序)
3.(这一步叫插值)
如下图,
ωkn
就是我们带入多项式的值:
第一个多项式
x=ω0n
,第二个是
x=ω1n
,第三个是
x=ω2n
….知道第n个是
x=ωn−1n
;
yi
就是带入不同的
x
的求出的多项式的值。
那么,通过上一步已经求出了两个多项式(两个乘数)的
最后推出
ai
的表达式
和第二步要求的式子几乎一样~再来一次FFT即可解决~
时间复杂度依然是
O(nlogn)
详见代码注释。
#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <algorithm>
#include <cstring>
#include <complex>
#define pi acos(-1)
#define N 200005
using namespace std;
complex<double> a[N],b[N],p[N];
int n,c[N];
char s[N];
void FFT(complex<double> x[],int n,int p)
{
//把原来依次排列的数变成叶子中的顺序
for (int i=0,t=0;i<n;i++)
{
if (i>t) swap(x[i],x[t]);
for (int j=n>>1;(t^=j)<j;j>>=1);
}
for (int m=2;m<=n;m<<=1) //枚举每一层
{
complex<double> wn(cos(p*2*pi/m),sin(p*2*pi/m));
for (int i=0;i<n;i+=m)
{
complex<double> w(1,0),u;
int k=m>>1;
for (int j=0;j<k;j++,w*=wn)
{
//蝴蝶操作
u=x[i+j+k]*w;
x[i+j+k]=x[i+j]-u;
x[i+j]=x[i+j]+u;
}
}
}
}
int main()
{
cin>>n;
scanf("%s",s);
for (int i=0;i<n;i++)
a[i]=s[n-i-1]-'0';
scanf("%s",s);
for (int i=0;i<n;i++)
b[i]=s[n-i-1]-'0';
//把长度变为2的幂次,方便FFT中的迭代
for (int j=n,i=1;(i>>2)<j;i<<=1)
n=i;
FFT(a,n,1),FFT(b,n,1);
for (int i=0;i<n;i++)
p[i]=a[i]*b[i];
//插值
FFT(p,n,-1);
for (int i=0;i<n;i++)
c[i]=p[i].real()/n+0.1;
int len=0;
//进位
for (int i=0;i<n;i++)
if (c[i])
len=i,c[i+1]+=c[i]/10,c[i]%=10;
for (int i=len;i>=0;i--)
printf("%d",c[i]);
return 0;
}