LuoguP3338 [ZJOI2014]力

题意

已知\(n\),\(q_1\),\(q_2\),...,\(q_n\)，定义

\(F_j= \sum_{i=1}^{j-1}\frac{q_i\times q_j}{(i-j)^2}- \sum_{i=j+1}^n\frac{q_i\times q_j}{(i-j)^2}\)，
\(E_i=\frac{F_i}{q_i}\)，求\(1\)~\(n\)所有\(E\)的值。

---

把\(F\)和\(E\)的式子联立起来可以得到：

\[ E_j=\sum_{i=1}^{j-1}\frac{q_i}{(i-j)^2}- \sum_{i=j+1}^{n}\frac{q_i}{(i-j)^2} \]

把式子分成两部分来看。

先看\(\sum_{i=1}^{j-1}\frac{q_i}{(i-j)^2}\)。设\(a_i=q_i\)，\(b_i=\frac{1}{i^2}\)，那么式子可以写成：\(\sum_{i=1}^{j-1}a_ib_{j-i}\)。

如果我们令\(a_0=b_0=0\)，那么式子可以进一步写成：\(\sum_{i=0}^{j}a_ib_{j-i}\)，可以发现是卷积的形式，用FFT做。

再看\(\sum_{i=j+1}^n\frac{q_i}{(i-j)^2}\)。设\(a_i=q_i\)，\(b_i=\frac{1}{i^2}\)，那么式子可以写成：

\[ \sum_{i=j+1}^na_ib_{j-i}=\sum_{i=1}^{n-j}a_{i+j}b_i \]

这里是一个很套路的翻转操作：引入\(a'_i=a_{n-i}\)，那么可以把式子写成\(\sum_{i=1}^{n-j}a'_{n-j-i}b_i\)。如果我们令\(b_0=0\)，那么式子可以进一步写成：\(\sum_{i=0}^{n-j}a'_{n-j-i}b_i\)，发现也是卷积的形式，用FFT做。

#include<bits/stdc++.h>
#define rg register
#define il inline
#define cn const
#define fp(i,a,b) for(rg int i=(a),ed=(b);i<=ed;++i)
using namespace std;
typedef cn int cint;
typedef double db;
cint maxn=100010; cn db pi=acos(-1.0);

int n,lim=1,l,r[maxn<<2];
db Sin[maxn<<2],Cos[maxn<<2];

struct cpx{db x,y;}a[maxn<<2],b[maxn<<2],c[maxn<<2];
il cpx operator+(cn cpx &x,cn cpx &y){return (cpx){x.x+y.x,x.y+y.y};}
il cpx operator-(cn cpx &x,cn cpx &y){return (cpx){x.x-y.x,x.y-y.y};}
il cpx operator*(cn cpx &x,cn cpx &y){return (cpx){x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x};}

il void fft(cpx *a,cint &f){
    fp(i,0,lim-1)if(i<r[i])swap(a[i],a[r[i]]);
    for(rg int md=1;md<lim;md<<=1){
        rg int len=md<<1; rg cpx wn=(cpx){Cos[len],f*Sin[len]};
        for(rg int l=0;l<lim;l+=len){
            rg cpx w=(cpx){1,0};
            for(rg int nw=0;nw<md;++nw,w=w*wn){
                rg cpx x=a[l+nw],y=a[l+nw+md]*w;
                a[l+nw]=x+y,a[l+nw+md]=x-y;
            }
        }
    }
}

int main(){
    scanf("%d",&n); fp(i,1,n)scanf("%lf",&a[i].x),b[i].x=1.0/i/i,c[n-i].x=a[i].x;
    while(lim<=n<<1)lim<<=1,++l;
    fp(i,0,lim-1)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    fp(i,0,lim)Cos[i]=cos(2*pi/i),Sin[i]=sin(2*pi/i);
    fft(a,1),fft(b,1),fft(c,1);
    fp(i,0,lim)a[i]=a[i]*b[i],c[i]=b[i]*c[i];
    fft(a,-1),fft(c,-1);
    fp(i,1,n)printf("%.3lf\n",(a[i].x-c[n-i].x)/lim);
    return 0;
}