bzoj3527: [Zjoi2014]力 卷積+FFT
阿新 • • 發佈:2018-11-19
ans second 既然 pair -i 一個 fine eal can
先寫個簡要題解:本來去桂林前就想速成一下FFT的,結果一直沒有速成成功,然後這幾天斷斷續續看了下,感覺可以寫一個簡單一點的題了,於是就拿這個題來寫,之前式子看著別人的題解都不太推的對,然後早上6點多推了一個多小時終於發現了一個很巧妙的方法,首先問題的關鍵在於後半個式子,因為顯然前半個式子很容易想到卷積的形式,那麽直接FFT就好了,但是後半部分不好考慮,一般肯定是通過類似換元的做法得出結論,所以到中間有一步就有點難度,那個地方我一直卡。後來突然想到,既然前半部分i<j時那麽好處理,那麽i>j的情況我把i和j分別用(n-i)和(n-j)代入不就轉化過去了,一步就過去了,然後直接就會發現他是卷積後的第(n-j)項,所以後面一半需要反轉a數組和反轉結果數組,這裏要註意下標,第二部分,第0項對應第n-1項。
代碼:
#include<bits/stdc++.h> #include<complex> #define db double #define ll long long #define mp make_pair #define fi first #define pb push_back #define se second #define rep(i,a,b)for(int i=a;i<=b;i++) using namespace std; const double pi=acos(-1); const int maxn=5e5+7; db spt(db x){return x*x;} int r[maxn]; complex<double>a[maxn],b[maxn],c[maxn],d[maxn]; db ans[maxn]; int M,N,n,l=0; db p[maxn]; void FFT(complex<double> f[],int op) { for(int i=0;i<N;i++)if(i<r[i])swap(f[i],f[r[i]]); for(int i=1;i<N;i<<=1) { complex<double >w(cos(pi/i),op*sin(pi/i)); for(int p=(i<<1),j=0;j<N;j+=p) { complex<double>W(1,0); for(int k=0;k<i;k++,W*=w) { complex<double>x=f[j+k],y=W*f[j+i+k]; f[j+k]=x+y;f[j+k+i]=x-y; } } } } int main() { scanf("%d",&n); for(int i=1;i<=n;i++) { scanf("%lf",&p[i]); } N=M=n-1; for(int i=0;i<n;i++){a[i]=p[i+1];if(i>=1)b[i]=1.0/spt(i);} M+=N; for(N=1;N<=M;N<<=1)l++; for(int i=0;i<N;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1)); FFT(a,1);FFT(b,1); for(int i=0;i<=N;i++)a[i]=a[i]*b[i]; FFT(a,-1); for(int i=0;i<n;i++)ans[i]+=a[i].real()*1.0/N; for(int i=0;i<n;i++){c[i]=p[n-i];if(i>=1)d[i]=1/spt(i);} FFT(c,1);FFT(d,1); for(int i=0;i<=N;i++)c[i]=c[i]*d[i]; FFT(c,-1); for(int i=0;i<n;i++)ans[i]-=c[n-i-1].real()*1.0/N; for(int i=0;i<n;i++)printf("%.3lf\n",ans[i]); return 0; }
bzoj3527: [Zjoi2014]力 卷積+FFT