[BZOJ3527][ZJOI2014]力 FFT+数学

Posted 2020-10-10 halfrot

题目链接：http://www.lydsy.com/JudgeOnline/problem.php?id=3527

首先卷积的形式是$h(i)=\sum_{i=0}^jf(i)g(i-j)$，如果我们可以把式子整理成这个样子再套上FFT就成功了。

$$E_i=\sum_{j<i}\frac{q_j}{(j-i)^2}-\sum_{j>i}\frac{q_j}{(i-j)^2}$$

$$E_i=\sum_{j=0}^{i-1}\frac{q_j}{(j-i)^2}^2-\sum_{j=0}^{n-i-1}\frac{q_{n-j}}{(n-i-j)^2}$$

令$f(i)=q_i$，$g(i)=\frac{1}{i^2}$，$t(i)=q_{n-i}$，则有

$$E_i=\sum_{j=0}^{i-1}f(i)g(i-j)-\sum_{j=0}^{n-i-1}t(j)g(n-i-j)$$

因为$j$无法取到$i$或者$n-i$，所以令$g(0)=0$来消除最后一项的影响。

 

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 #include<cmath>
 5 using namespace std;
 6 const double pi=acos(-1.0);
 7 struct complex{
 8     double r,i;
 9     complex (double _r=0,double _i=0){
10         r=_r,i=_i;
11     }
12     complex operator + (const complex &_){
13         return complex(r+_.r,i+_.i);
14     }
15     complex operator - (const complex &_){
16         return complex(r-_.r,i-_.i);
17     }
18     complex operator * (const complex &_){
19         return complex(r*_.r-i*_.i,r*_.i+_.r*i);
20     }
21 };
22 void reverse(complex y[],int len){
23     for(int i=1,j=len>>1,k;i+1<len;i++){
24         if(i<j) swap(y[i],y[j]);
25         k=len>>1;
26         while(j>=k){
27             j-=k;
28             k>>=1;
29         }
30         if(j<k) j+=k;
31     }
32 }
33 void FFT(complex y[],int len,int on){
34     reverse(y,len);
35     for(int h=2;h<=len;h<<=1){
36         complex wn(cos(-on*2*pi/h),sin(-on*2*pi/h));
37         for(int j=0;j<len;j+=h){
38             complex w(1,0);
39             for(int k=j;k<j+(h>>1);k++){
40                 complex u=y[k],t=w*y[k+(h>>1)];
41                 y[k]=u+t;
42                 y[k+(h>>1)]=u-t;
43                 w=w*wn;
44             }
45         }
46     }
47     if(on==-1) for(int i=0;i<len;i++) y[i].r/=len;
48 }
49 double q[100010],ans[100010];
50 complex a[300010],b[300010],c[300010];
51 int main(){
52     int n,len=1;
53     scanf("%d",&n);
54     while(len<n) len<<=1;len<<=1;
55     for(int i=0;i<n;i++) scanf("%lf",&q[i]);
56     for(int i=0;i<n;i++) a[i]=complex(q[i],0);
57     for(int i=1;i<n;i++) b[i]=complex(1.0/i/i,0);
58     for(int i=n;i<len;i++) a[i]=b[i]=complex();
59     b[0]=complex();
60     FFT(a,len,1);
61     FFT(b,len,1);
62     for(int i=0;i<len;i++) c[i]=a[i]*b[i];
63     FFT(c,len,-1);
64     for(int i=0;i<n;i++) ans[i]=c[i].r;
65     for(int i=0;i<n;i++) a[i]=complex(q[n-i-1],0);
66     for(int i=1;i<n;i++) b[i]=complex(1.0/i/i,0);
67     for(int i=n;i<len;i++) a[i]=b[i]=complex();
68     b[0]=complex();
69     FFT(a,len,1);
70     FFT(b,len,1);
71     for(int i=0;i<len;i++) c[i]=a[i]*b[i];
72     FFT(c,len,-1);
73     for(int i=0;i<n;i++) ans[i]-=c[n-i-1].r;
74     for(int i=0;i<n;i++) printf("%.3lf\n",ans[i]);
75     return 0;
76 }