BZOJ 2179: FFT快速傅立叶

2179: FFT快速傅立叶

Time Limit: 10 Sec  Memory Limit: 259 MB
Submit: 3425  Solved: 1795
[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.（这一步叫插值） 
如下图， 
ai表示多项式的系数； 
ωkn就是我们带入多项式的值： 
第一个多项式x=ω0n，第二个是x=ω1n，第三个是x=ω2n….知道第n个是x=ωn−1n； 
yi就是带入不同的x的求出的多项式的值。 

那么，通过上一步已经求出了两个多项式（两个乘数）的yi，现在我们把对应的yi乘起来就是乘积的多项式的Yi，我们要求的是这个多项式的系数即ai，因此我们只要求出ωkn矩阵的逆矩阵即可。

最后推出ai的表达式 
 
和第二步要求的式子几乎一样~再来一次FFT即可解决~ 
时间复杂度依然是O(nlogn)

转载出处：http://blog.csdn.net/regina8023/article/details/44908753

附代码模板

#include<cmath>#include<ctime>#include<cstdio>#include<cstdlib>#include<cstring>#include<complex>#include<iostream>#include<algorithm>#include<iomanip>#include<vector>#include<string>#include<queue>#include<set>#include<map>#define pi acos(-1)using namespace std;inline int read(){int x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}while(ch<='9'&&ch>='0'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}return f*x;}const int N=1<<17;complex<double> a[N],b[N],t[N];char s[N];int c[N],r[N];void FFT(complex<double> x[],int n,int p){for(int i=0;i<n;i++)if(i<r[i])swap(x[i],x[r[i]]);for(int m=2;m<=n;m<<=1){complex<double> wn(cos(2*p*pi/m),sin(2*p*pi/m));for(int i=0;i<n;i+=m){complex<double> w(1,0),temp;int k=m>>1;//不知道为什么这个在这里定义跑得快for(int j=0;j<k;j++,w*=wn){temp=x[i+j+k]*w;x[i+j+k]=x[i+j]-temp;x[i+j]=x[i+j]+temp;}}}}int main(){int n=read();scanf("%s",s+1);for(int i=0;i<n;i++)a[i]=s[n-i]-'0';scanf("%s",s+1);for(int i=0;i<n;i++)b[i]=s[n-i]-'0';for(int i=1,j=n;(i>>2)<j;i<<=1)n=i;for(int i=0;i<n;i++) r[i]=(i&1)*(n>>1)+(r[i>>1]>>1); FFT(a,n,1);FFT(b,n,1);for(int i=0;i<n;i++)t[i]=a[i]*b[i];FFT(t,n,-1);for(int i=0;i<n;i++)c[i]=t[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]);puts("");    return 0;}/*134*/

2017.9.3更新模板

采用重载复数，快很很多

#include<cmath>#include<ctime>#include<cstdio>#include<cstring>#include<cstdlib>#include<complex>#include<iostream>#include<algorithm>#include<iomanip>#include<vector>#include<string>#include<bitset>#include<queue>#include<map>#include<set>using namespace std;typedef long long ll;typedef double db;inline int read(){int x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}while(ch<='9'&&ch>='0'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}return x*f;}void print(int x){if(x<0)putchar('-'),x=-x;if(x>=10)print(x/10);putchar(x%10+'0');}#define pi acos(-1)const int N=200100;struct cp{db real,image;friend cp operator -(const cp &x,const cp &y){return (cp){x.real-y.real,x.image-y.image};}friend cp operator +(const cp &x,const cp &y){return (cp){x.real+y.real,x.image+y.image};}friend cp operator *(const cp &x,const cp &y){return (cp){x.real*y.real-x.image*y.image,x.real*y.image+x.image*y.real};}}a[N],b[N],c[N];int r[N],ans[N];void fft(cp x[],int n,int p){register int i,j,k,m;for(i=0;i<n;++i)if(r[i]<i)swap(x[i],x[r[i]]);for(m=2;m<=n;m<<=1){cp wn=(cp){cos(2*pi*p/m),sin(2*pi*p/m)};for(i=0;i<n;i+=m){cp w=(cp){1,0},tmp;k=m>>1;for(j=0;j<k;++j,w=w*wn){tmp=x[i+j+k]*w;x[i+j+k]=x[i+j]-tmp;x[i+j]=x[i+j]+tmp;}}}}char s[N];int main(){int n=read();register int i,j,len=0;scanf("%s",s);for(i=0;i<n;++i)a[i].real=s[n-i-1]-'0';scanf("%s",s);for(i=0;i<n;++i)b[i].real=s[n-i-1]-'0';for(i=1,j=n;(i>>2)<j;i<<=1)n=i;for(i=0;i<n;++i)r[i]=((i&1)*(n>>1))+(r[i>>1]>>1);fft(a,n,1);fft(b,n,1);for(i=0;i<n;++i)c[i]=a[i]*b[i];fft(c,n,-1);for(i=0;i<n;++i)ans[i]=c[i].real/n+0.1;for(i=0;i<n;++i)if(ans[i]){len=i;ans[i+1]+=ans[i]/10;ans[i]%=10;}for(i=len;~i;i--)putchar(ans[i]+'0');puts("");return 0;}