1. 程式人生 > >快速傅立葉變換(FFT)和數論變換(NTT)模板

快速傅立葉變換(FFT)和數論變換(NTT)模板

#include<bits/stdc++.h>
#define ll long long
#define dob complex<double>
const int N=120010;
const double pi=acos(-1.0);
using namespace std;
inline int rd()
{
    int x=0,res=1;char ch=getchar();
    while(ch>'9'||ch<'0'){if(ch=='-')res*=-1;ch=getchar();}
    while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
    return x*res;
}
int n,m,L,R[N];
dob a[N],b[N];
void FFT(dob *A,int n,int f)
{
    for(int i=0;i<n;i++)if(i<R[i])swap(A[i],A[R[i]]);
    for(int i=1;i<n;i<<=1)
    {
        dob wn(cos(pi/i),sin(f*pi/i)),x,y;
        for(int j=0;j<n;j+=(i<<1))
        {
            dob w(1,0);
            for(int k=0;k<i;k++,w*=wn)
            {
                x=A[j+k];y=w*A[j+i+k];
                A[j+k]=x+y;
                A[j+i+k]=x-y;
            }
        }
    }
}
int main()
{
    while(~scanf("%d%d",&n,&m))
    {
        for(int i=0;i<n;i++)a[i]=rd();//複數一定要用讀入掛讀入
        for(int i=0;i<m;i++)b[i]=rd();
        int len=1,s=n+m;
        L=0;
        for(;len<s;len<<=1)++L;
        for(int i=n;i<len;i++)a[i]=0;
        for(int i=n;i<len;i++)b[i]=0;
        for(int i=0;i<len;i++)R[i]=(R[i>>1]>>1)|((i&1)<<(L-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<=s;i++)printf("%d ",int(a[i].real()/len+0.5));puts("");//idft以後要除len
    }
}
#include<bits/stdc++.h>
#define ll long long
const int N=262144;
const ll MOD=50000000001507329LL;//998244353 1004535809
using namespace std;
int n,m;
ll a[N],b[N],x[N],y[N];
ll wn[25];
ll Mul(ll x,ll y)//乘法超ll用快速乘,主函式也需要用
{
    ll ans=(x*y-(ll)((long double)x/MOD*y+1e-8)*MOD);
    return ans<0?ans+MOD:ans;
}
ll Qpow(ll a,ll b,ll M)
{
    ll ans=1;a%=M;
    while(b)
    {
        if(b&1) ans=Mul(ans,a);
        a=Mul(a,a);
        b>>=1;
    }
    return ans;
}
void Getwn()//主函式預處理getwn()
{
    for(int i=0;i<25;i++)
    {
        wn[i]=Qpow(3,(MOD-1)/(1<<i),MOD);
    }
}
void NTT(ll *x,int n,int rev)
{
    int i,j,k,ds;
    ll w,u,v;
    for(i=1,j=n>>1,k=n>>1;i<n-1;i++,k=n>>1)
    {
        if(i<j) swap(x[i],x[j]);
        while(j>=k) j-=k,k>>=1;
        if(j<k) j+=k;
    }
    for(i=2,ds=1;i<=n;i<<=1,ds++)
    {
        for(j=0;j<n;j+=i)
        {
            w=1;
            for(k=j;k<j+i/2;k++)
            {
                u=x[k];
                v=Mul(w,x[k+i/2]);
                x[k]=(u+v)%MOD;
                x[k+i/2]=(u-v+MOD)%MOD;
                w=Mul(w,wn[ds]);
            }
        }
    }
    if(rev==-1)
    {
        for(i=1;i<n/2;i++) swap(x[i],x[n-i]);
        w=Qpow(n,MOD-2,MOD);
        for(i=0;i<n;i++) x[i]=Mul(x[i],w);
    }
}
int main()
{
    Getwn();
    while(~scanf("%d%d",&n,&m))
    {
        for(int i=0;i<n;i++)scanf("%lld",&a[i]);
        for(int i=0;i<m;i++)scanf("%lld",&b[i]);
        int len=1,s=n+m;
        while(len<s)len<<=1;
        for(int i=n;i<len;i++)a[i]=0;
        for(int i=m;i<len;i++)b[i]=0;
        NTT(a,len,1);NTT(b,len,1);
        for(int i=0;i<len;i++)a[i]=Mul(a[i],b[i]);
        NTT(a,len,-1);
        for(int i=0;i<=s;i++)printf("%lld ",a[i]);puts("");
    }
}