2179: FFT快速傅立叶
Time Limit: 10 Sec Memory Limit: 259 MB
Submit: 3425 Solved: 1795
[Submit][Status][Discuss]
2179: FFT快速傅立叶
Time Limit: 10 Sec Memory Limit: 259 MBSubmit: 3425 Solved: 1795
[Submit][Status][Discuss]
Description
Input
Output
Sample Input
3
4
Sample Output
数据范围:
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;
}
/*
1
3
4
*/
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;
}