题目传送门:
https://www.luogu.org/problemnew/show/P3338
题意:
给定一个公式,求值。
思路:
考虑化简公式。
F j = ∑ i < j q i ∗ q j ( i − j ) 2 − ∑ i > j q i ∗ q j ( i − j ) 2 F_j=\sum\limits_{i<j}\dfrac{q_i*q_j}{(i-j)^2}-\sum\limits_{i>j} \dfrac{q_i*q_j}{(i-j)^2} Fj=i<j∑(i−j)2qi∗qj−i>j∑(i−j)2qi∗qj
F
j
=
∑
i
=
1
j
−
1
q
i
∗
q
j
(
i
−
j
)
2
−
∑
i
=
j
+
1
n
q
i
∗
q
j
(
i
−
j
)
2
F_j=\sum\limits_{i=1}^{j-1}\dfrac{q_i*q_j}{(i-j)^2}-\sum\limits_{i=j+1}^{n}\dfrac{q_i*q_j}{(i-j)^2}
Fj=i=1∑j−1(i−j)2qi∗qj−i=j+1∑n(i−j)2qi∗qj
带入
E
j
=
F
j
q
j
E_j=\dfrac{F_j}{q_j}
Ej=qjFj,从而消除
q
j
q_j
qj,得:
E
j
=
∑
i
=
1
j
−
1
q
i
(
i
−
j
)
2
−
∑
i
=
j
+
1
n
q
i
(
i
−
j
)
2
E_j=\sum\limits_{i=1}^{j-1}\dfrac{q_i}{(i-j)^2}-\sum\limits_{i=j+1}^{n}\dfrac{q_i}{(i-j)^2}
Ej=i=1∑j−1(i−j)2qi−i=j+1∑n(i−j)2qi
令
A
i
=
1
i
2
A_i=\dfrac{1}{i^2}
Ai=i21,
q
i
′
q'_i
qi′为
q
i
q_i
qi反向取的序列,则有:
E j = ∑ i = 1 j − 1 q i ( i − j ) 2 − ∑ i = 1 j − 1 q i ′ ( i − j ) 2 E_j=\sum\limits_{i=1}^{j-1}\dfrac{q_i}{(i-j)^2}-\sum\limits_{i=1}^{j-1}\dfrac{q'_i}{(i-j)^2} Ej=i=1∑j−1(i−j)2qi−i=1∑j−1(i−j)2qi′
E j = ∑ i = 1 j − 1 q i ( j − i ) 2 − ∑ i = 1 j − 1 q i ′ ( j − i ) 2 E_j=\sum\limits_{i=1}^{j-1}\dfrac{q_i}{(j-i)^2}-\sum\limits_{i=1}^{j-1}\dfrac{q'_i}{(j-i)^2} Ej=i=1∑j−1(j−i)2qi−i=1∑j−1(j−i)2qi′
E
j
=
∑
i
=
1
j
−
1
q
i
∗
A
j
−
i
−
∑
i
=
1
j
−
1
∗
q
i
′
∗
A
j
−
i
E_j=\sum\limits_{i=1}^{j-1}q_i*A_{j-i}- \sum\limits_{i=1}^{j-1}*q'_i*A_{j-i}
Ej=i=1∑j−1qi∗Aj−i−i=1∑j−1∗qi′∗Aj−i
容易发现这是一个类似于卷积的形式,便可以用
F
F
T
FFT
FFT来加速。
代码:
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define LL long long
using namespace std;
const double pai=acos(-1.0);
struct comp
{
double x,y;
comp(double xx=0,double yy=0):x(xx),y(yy) {}
friend comp operator+(const comp &x,const comp &y) {return comp(x.x+y.x,x.y+y.y);}
friend comp operator-(const comp &x,const comp &y) {return comp(x.x-y.x,x.y-y.y);}
friend comp operator*(const comp &a,const comp &b) {return comp(a.x*b.x-a.y*b.y,a.x*b.y+b.x*a.y);}
}a[500000],b[500000];
int limit=1,n,l=0,ma=0;
int r[500000];
double d[500000][2],p[500000],f[500000];
void init()
{
while(limit<=ma)
limit<<=1,l++;
for(int i=1;i<=limit;i++)
r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
}
void FFT(comp *now,int ty)
{
for(int i=0;i<limit;i++)
if(i<r[i]) swap(now[i],now[r[i]]);
for(int mid=1;mid<limit;mid<<=1)
{
comp wn(cos(pai/mid),ty*sin(pai/mid));
for(int j=0,R=(mid<<1);j<limit;j+=R)
{
comp w(1,0);
for(int k=0;k<mid;k++,w=w*wn)
{
comp x=now[j+k],y=w*now[j+k+mid];
now[j+k]=x+y;
now[j+k+mid]=x-y;
}
}
}
}
void work(int id)
{
for(int i=0;i<=limit;i++)
a[i]=(comp){0,0},b[i]=(comp){0,0};
for(int i=1;i<=n;i++)
a[i].x=d[i][id],b[i].x=p[i];
FFT(a,1);
FFT(b,1);
for(int i=0;i<=limit;i++)
a[i]=a[i]*b[i];
FFT(a,-1);
if(id==1)
for(int i=0;i<=ma;i++)
f[i]=a[i].x/limit;
else
for(int i=0;i<=ma;i++)
f[i]-=a[n-i+1].x/limit;
}
int main()
{
scanf("%d",&n);
for(int i=1;i<=n;i++)
{
scanf("%lf",&d[i][1]);
p[i]=1.0/i/i;
}
for(int i=1;i<=n;i++)
d[i][2]=d[n-i+1][1];
ma=(n<<1);
init();
work(1);
work(2);
for(int i=1;i<=n;i++)
printf("%.3lf\n",f[i]);
}