FFT/快速傅立叶变换
本文最后更新于:2017年8月18日 下午
好了,这图是用来搞笑
的~
具体内容
UOJ 34/LOJ 108 多项式乘法
时间限制:1s
空间限制:256MB
这是一道模板题。
给你两个多项式,请输出乘起来后的多项式。
输入格式
第一行两个整数n和m,分别表示两个多项式的次数。
第二行n+1个整数,分别表示第一个多项式的0 到n 次项前的系数。
第三行m+1个整数,分别表示第一个多项式的0到m次项前的系数。
输出格式
一行n+m+1个整数,分别表示乘起来后的多项式的0到n+m次项前的系数。
样例一
input
1 2
1 2
1 2 1
output
1 4 5 2
explanation
$(1+2x)⋅(1+2x+x^2)=1+4x+5x^2+2x^3$
限制与约定
$0≤n,m≤10^5$,保证输入中的系数大于等于0且小于等于9。
Code
//high* or many* UOJ 34
#include<iostream>
#include<cstdio>
#include<complex>
#include<cmath>
#include<cstring>
using namespace std;
typedef complex<double> E;
const int maxn=3e5+10;//max len
const double PI=acos(-1),eps=1e-8;
E a[maxn],b[maxn];
int n,m;
double x;
inline void Rader(E*,int);
void FFT(E*,int,int);
int main()
{
scanf("%d%d",&n,&m);
memset(a,0,sizeof(a)),memset(b,0,sizeof(b));
for(int i=0;i<=n;i++) scanf("%lf",&x),a[i].real(x);
for(int i=0;i<=m;i++) scanf("%lf",&x),b[i].real(x);
int len=1;
while(len<=n+m) len<<=1;
FFT(a,len,1),FFT(b,len,1);
for(int i=0;i<len;i++) a[i]*=b[i];
FFT(a,len,-1);
for(int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].real()+0.5));
return 0;
}
inline void Rader(E *y,int n)
{
for(int i=0,j=0;i<n;i++)
{
if(i<j) swap(y[i],y[j]);
int k=n;
while(j&(k>>=1)) j&=~k;//is & not &&
j|=k;
}
}
void FFT(E *y,int n,int inv)
{
Rader(y,n);
for(int h=2;h<=n;h<<=1)
{
E Wn(cos(2*PI/h),inv*sin(2*PI/h));//e^(2*PI/h)=cos(2*PI/h)+i*sin(2*PI/h)
for(int i=0;i<n;i+=h)
{
E W(1,0);
for(int j=i;j<i+h/2;j++)
{
E u=y[j],t=W*y[j+h/2];
y[j]=u+t;
y[j+h/2]=u-t;
W*=Wn;//旋转因子
}
}
}
if(inv==-1) for(int i=0;i<n;i++) y[i]/=n;
}
hdu1402 A * B Problem Plus
给出若干组自然数A,B(A,B的位数不大于50000),求$A \times B$.
Code
#include<iostream>
#include<cstdio>
#include<complex>
#include<cmath>
#include<cstring>
#define max(a,b) ((a)>(b)?(a):(b))
using namespace std;
typedef complex<double> E;
const int maxn=3e5+10;
const double PI=acos(-1),eps=1e-8;
string str1,str2;
E a[maxn],b[maxn];
int ans[maxn];
bool flag=0;
void Rader(E*,int);
void FFT(E*,int,int);
int main()
{
// freopen("data1.in","r",stdin);
while(getline(cin,str1)&&getline(cin,str2))
{
flag=0;
if(str1[0]=='-') str1.erase(0,1),flag^=1;
if(str2[0]=='-') str2.erase(0,1),flag^=1;
int len=1,len1=str1.size(),len2=str2.size();
while(len<(max(len1,len2)<<1)) len<<=1;
memset(a,0,sizeof(a));
memset(b,0,sizeof(b));
for(int i=0;i<len1;i++) a[i].real(str1[len1-i-1]-'0');
for(int i=0;i<len2;i++) b[i].real(str2[len2-i-1]-'0');
for(int i=len1;i<len;i++) a[i].real(0.0);
for(int i=len2;i<len;i++) b[i].real(0.0);
FFT(a,len,1),FFT(b,len,1);
for(int i=0;i<len;i++) a[i]*=b[i];
FFT(a,len,-1);
for(int i=0;i<len;i++) ans[i]=(int)(a[i].real()+0.5);
for(int i=0;i<len;i++)
{
ans[i+1]+=ans[i]/10;
ans[i]%=10;
}
while(--len&&ans[len]==0);//不能把0全部去掉,因为有可能乘数为0.
if(flag) putchar('-');
for(int i=len;i>=0;--i) putchar(ans[i]+'0');
putchar('\n');
}
return 0;
}
void Rader(E *y,int n)
{
for(int i=0,j=0;i<n;i++)
{
if(i<j) swap(y[i],y[j]);
int k=n;
while(j&(k>>=1)) j&=~k;
j|=k;
}
}
void FFT(E *y,int n,int inv)
{
Rader(y,n);
for(int h=2;h<=n;h<<=1)
{
E Wn(cos(2*PI/h),inv*sin(2*PI/h));
for(int i=0;i<n;i+=h)
{
E W(1,0);
for(int j=i;j<i+h/2;j++)
{
E u=y[j],t=W*y[j+h/2];
y[j]=u+t;
y[j+h/2]=u-t;
W*=Wn;
}
}
}
if(inv==-1) for(int i=0;i<n;i++) y[i]/=n;
}
本博客所有文章除特別声明外,均采用
CC BY 4.0
许可协议。转载请注明来源
SpaceSkyNet
!