BZOJ3527[Zjoi2014]力——FFT
题目描述
给出n个数qi,给出Fj的定义如下:
令Ei=Fi/qi,求Ei.
输入
第一行一个整数n。
接下来n行每行输入一个数,第i行表示qi。
n≤100000,0<qi<1000000000
输出
n行,第i行输出Ei。与标准答案误差不超过1e-2即可。
样例输入
5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880
样例输出
-16838672.693
3439.793
7509018.566
4595686.886
10903040.872
3439.793
7509018.566
4595686.886
10903040.872
将$q_{j}$移到等式左边,可以得到:$E_{k}=\sum\limits_{i<k}^{ }\frac{q_{i}}{(i-k)^2}-\sum\limits_{i>k}^{ }\frac{q_{i}}{(i-k)^2}$
将等式右边分成两部分看。
对于第一部分设$j=k-i$,可以得到:$E_{k}'=\sum\limits_{i+j==k}^{ }q_{i}*\frac{1}{j^2}$,设$F(i)=q_{i},G(i)=\frac{1}{i^2},H(i)=E_{i}'$,将$F$和$G$卷积即可得到$H$。
对于第二部分设$j=i-k$,可以得到:$E_{k}''=\sum\limits_{i-j==k}^{ }q_{i}*\frac{1}{j^2}$,即$E_{k}''=\sum\limits_{j+n-i+1==n-k+1}^{ }q_{i}*\frac{1}{j^2}$。
将$n-k+1$和$n-i+1$看作新的$k$和$i$,也就是将第一部分中的$H$和$F$翻转,然后再卷积即可。即$F(n-i+1)=q_{i},H(n-i+1)=E_{i}''$。
将两部分对应做减法即可。
#include<set> #include<map> #include<queue> #include<stack> #include<cmath> #include<vector> #include<cstdio> #include<bitset> #include<cstring> #include<iostream> #include<algorithm> #define ll long long using namespace std; const double pi=acos(-1.0); struct lty { double x,y; lty(double X=0,double Y=0){x=X,y=Y;} }f[400000],g[400000],h[400000]; int n,mask=1; double ans[400000]; double q[400000]; lty operator +(lty a,lty b){return lty(a.x+b.x,a.y+b.y);} lty operator -(lty a,lty b){return lty(a.x-b.x,a.y-b.y);} lty operator *(lty a,lty b){return lty(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);} void DFT(lty *a,int len) { for(int i=0,k=0;i<len;i++) { if(i>k) { swap(a[i],a[k]); } for(int j=len>>1;(k^=j)<j;j>>=1); } for(int k=2;k<=len;k<<=1) { int t=k>>1; lty x(cos(pi/t),sin(pi/t)); for(int i=0;i<len;i+=k) { lty w(1,0); for(int j=i;j<i+t;j++) { lty res=a[j+t]*w; a[j+t]=a[j]-res; a[j]=a[j]+res; w=w*x; } } } } void IDFT(lty *a,int len) { for(int i=0,k=0;i<len;i++) { if(i>k) { swap(a[i],a[k]); } for(int j=len>>1;(k^=j)<j;j>>=1); } for(int k=2;k<=len;k<<=1) { int t=k>>1; lty x(cos(pi/t),-1.0*sin(pi/t)); for(int i=0;i<len;i+=k) { lty w(1,0); for(int j=i;j<i+t;j++) { lty res=a[j+t]*w; a[j+t]=a[j]-res; a[j]=a[j]+res; w=w*x; } } } for(int i=0;i<len;i++) { a[i].x=a[i].x/len; } } int main() { scanf("%d",&n); for(int i=1;i<=n;i++) { scanf("%lf",&q[i]); } while(mask<=(n<<1)) { mask<<=1; } for(int i=1;i<=n;i++) { f[i].x=q[i]; } for(int i=1;i<=n;i++) { g[i].x=(double)1/((double)i*i); } DFT(f,mask); DFT(g,mask); for(int i=0;i<mask;i++) { h[i]=g[i]*f[i]; } IDFT(h,mask); for(int i=1;i<=n;i++) { ans[i]=h[i].x; } for(int i=0;i<mask;i++) { f[i].x=f[i].y=0; g[i].x=g[i].y=0; } for(int i=1;i<=n;i++) { f[i].x=q[n-i+1]; } for(int i=1;i<=n;i++) { g[i].x=(double)1/((double)i*i); } DFT(f,mask); DFT(g,mask); for(int i=0;i<mask;i++) { h[i]=g[i]*f[i]; } IDFT(h,mask); for(int i=1;i<=n;i++) { ans[i]-=h[n-i+1].x; } for(int i=1;i<=n;i++) { printf("%.3f\n",ans[i]); } }