「题解」快速傅立叶之二

「题解」快速傅立叶之二

前言

这篇博客是我从晴歌。大佬这里扒下来的。

写得非常好，如果大佬发现了，也请见谅！

题目描述

点这里

给出序列 \(a[0],a[1],...,a[n-1]\) 和 \(b[0],b[1],...,b[n-1]\) .

\[c[k]=\sum_{i=k}^{n-1}a[i]b[i-k]\] .

求序列 \(c[]\) .

题解

这题就是BZOJ_3527[ZJOI2014]力(FFT+卷积)的后半段...

我们来重新分析一下.

首先我们要知道卷积的标准形式:

\[c[i]=\sum_{j=0}^ia[j]b[i-j]\]

很明显这道题并不是卷积的形式,因为卷积是和一定,二这道题却是差一定.

我们其实可以画画图(我脑洞大)...

然后可以发现差一定的时候就是你+1,我也+1,你-1,我也-1.

但是如果我们把其中一个序列倒过来,就变成了你+1,我-1,你-1,我+1,就变成和一定的了!这一点灰常重要!

然后上次我推的那个太不自然,我们这次好好分析一下.

1.把a倒置.

把a倒置之前原式为(我们这里令 \(n=n-1\) ,序列就是 \(0~n\) ,方便一些)

\[\sum_{j=k}^na[i]b[i-k]\]

我们设倒置之后的序列为 \(a'[]\) ,则有

\[原式\Longleftrightarrow\sum_{i=k}^na'[n-i]b[i-k]\]

换元,得到:

\[\sum_{i=0}^{n-k}a'[n-(i+k)]b[(i+k)-k]\]

即

\[\sum_{i=0}^{n-k}a'[n-i-k]b[i]\]

也就是:

\[c[k]=\sum_{i=0}^{n-k}a'[n-i-k]b[i]\]

如果我们设 \(A[k]=\sum_{i=1}^ka'[k-i]b[i]\) ,那么就有

\[c[k]=A[n-k]\]

这样我们求个卷积,然后倒过来输出就好了.

2.把b倒置

在网上看到好几篇题解都说是倒置b,但是自己推了好长时间都没有推出来,最后发现其中有一点奥妙...

倒置之前原式

\[\sum_{i=k}^na[i]b[i-k]\]

我们设倒置之后的序列为 \(b'[]\) ,则有

\[原式\Longleftrightarrow\sum_{i=k}^na[i]b'[n-i+k]\]

换元,得到

\[\sum_{i=0}^{n-k}a[i+k]b'[n_(i+k)+k]\]

也就是

\[\sum_{i=0}^{n-k}a[i+k]b'[n-i]\].

可以发现和是定值 \(n+k\) ,但是循环上界却只有 \(n-k\) .

我们想要得到的应该是

\[\sum_{i=0}^{n+k}a[i+k]b'[n-i]\].

我们得到的式子少了一部分.但是观察可以发现,我们得到的式子的循环上界是 \(n-k\) ,对应 \(a[n]b'[k]\) .

继续向上的 \(a[i]\) 都为 \(0\) ,而且都后的 \(b[i]\) 会越界 \(b[负数]\) .

所以这个就可以表示一个卷积了.

\[c[k]=\sum_{i=0}^{n+k}a[n+k-i]b'[i]\]

这个式子是根据原式表示一个卷积二构造出来的等价的式子,只是看起来比较方便而已.

我们设 \(B[i]=\sum_{i=0}^ka[i]b[k-i]\) .

这样就可以得到：

\[c[k]=B[n+k]\]

接下来是代码：

#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;

#define rep(i,__l,__r) for(register int i=__l,i##_end_=__r;i<=i##_end_;++i)
#define fep(i,__l,__r) for(register int i=__l,i##_end_=__r;i>=i##_end_;--i)
#define writc(a,b) fwrit(a),putchar(b)
#define mp(a,b) make_pair(a,b)
#define ft first
#define sd second
#define LL long long
#define ull unsigned long long
#define pii pair<int,int>
#define Endl putchar('\n')
// #define FILEOI
// #define int long long

#ifdef FILEOI
    #define MAXBUFFERSIZE 500000
    inline char fgetc(){
        static char buf[MAXBUFFERSIZE+5],*p1=buf,*p2=buf;
        return p1==p2&&(p2=(p1=buf)+fread(buf,1,MAXBUFFERSIZE,stdin),p1==p2)?EOF:*p1++;
    }
    #undef MAXBUFFERSIZE
    #define cg (c=fgetc())
#else
    #define cg (c=getchar())
#endif
template<class T>inline void qread(T& x){
    char c;bool f=0;
    while(cg<'0'||'9'<c)f|=(c=='-');
    for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
    if(f)x=-x;
}
inline int qread(){
    int x=0;char c;bool f=0;
    while(cg<'0'||'9'<c)f|=(c=='-');
    for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
    return f?-x:x;
}
template<class T,class... Args>inline void qread(T& x,Args&... args){qread(x),qread(args...);}
template<class T>inline T Max(const T x,const T y){return x>y?x:y;}
template<class T>inline T Min(const T x,const T y){return x<y?x:y;}
template<class T>inline T fab(const T x){return x>0?x:-x;}
inline int gcd(const int a,const int b){return b?gcd(b,a%b):a;}
inline void getInv(int inv[],const int lim,const int MOD){
    inv[0]=inv[1]=1;for(int i=2;i<=lim;++i)inv[i]=1ll*inv[MOD%i]*(MOD-MOD/i)%MOD;
}
template<class T>void fwrit(const T x){
    if(x<0)return (void)(putchar('-'),fwrit(-x));
    if(x>9)fwrit(x/10);putchar(x%10^48);
}
inline LL mulMod(const LL a,const LL b,const LL mod){//long long multiplie_mod
    return ((a*b-(LL)((long double)a/mod*b+1e-8)*mod)%mod+mod)%mod;
}

class fft_task{
private:
    static const int MAXN=1<<20;
    const double Pi=acos(-1.0);
    struct cplx{
        double vr,vi;//实部和虚部
        cplx(const double R=0,const double I=0):vr(R),vi(I){}//构造函数
        //------------------overload----------------//
        cplx operator + (const cplx a)const{return cplx(vr+a.vr,vi+a.vi);}//重载加法
        cplx operator - (const cplx a)const{return cplx(vr-a.vr,vi-a.vi);}
        cplx operator * (const cplx a)const{return cplx(vr*a.vr-vi*a.vi,vr*a.vi+a.vr*vi);}
        cplx operator / (const double var)const{return cplx(vr/var,vi/var);}
    };

    cplx a[MAXN+5],b[MAXN+5];
    int revi[MAXN+5],n,m;

    inline void init(){
        for(int i=0;i<n;++i)revi[i]=(revi[i>>1]>>1)|((i&1)?n>>1:0);
    }

    inline void fft(cplx* f,const short opt=1){
        for(int i=0;i<n;++i)if(i<revi[i])
            swap(f[i],f[revi[i]]);
        for(int p=2;p<=n;p<<=1){
            int len=p/2;
            cplx tmp(cos(Pi/len),opt*sin(Pi/len));
            for(int k=0;k<n;k+=p){
                cplx buf(1,0);
                for(int l=k;l<k+len;++l,buf=buf*tmp){
                    cplx tt=buf*f[len+l];
                    f[len+l]=f[l]-tt;
                    f[l]=f[l]+tt;
                }
            }
        }
        if(opt==-1)for(int i=0;i<n;++i)f[i]=f[i]/n;
    }
public:
    inline void launch(){
        n=m=qread()-1;
        int t=n;
        for(int i=0;i<=n;++i)scanf("%lf %lf",&a[i].vr,&b[n-i].vr);
        for(m+=n,n=1;n<=m;n<<=1);
        init();
        fft(a),fft(b);
        for(int i=0;i<n;++i)a[i]=a[i]*b[i];
        fft(a,-1);
        for(int i=t;i<=t*2;++i)printf("%d\n",(int)(a[i].vr+0.5));
    }
}This;

signed main(){
#ifdef FILEOI
    freopen("1.in","r",stdin);
    freopen("file.out","w",stdout);
#endif
    This.launch();
    return 0;
}