FFT/快速傅立叶变换


本文最后更新于:2017年8月18日 下午

FFT
好了,这图是用来搞笑的~

具体内容

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;
}

文章作者: SpaceSkyNet
版权声明: 本博客所有文章除特別声明外,均采用 CC BY 4.0 许可协议。转载请注明来源 SpaceSkyNet !
  目录
评论