「题解」快速傅立叶之二
「题解」快速傅立叶之二
前言
这篇博客是我从晴歌。大佬这里扒下来的。
写得非常好,如果大佬发现了,也请见谅!
题目描述
给出序列 \(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+卷积)](http://www.cnblogs.com/Sunnie69/p/5574222.html])的后半段...
我们来重新分析一下.
首先我们要知道卷积的标准形式:
$$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;
}