Luogu P3338 [ZJOI2014]力

题目大意

给定 $n$ 个数 $q_i$，定义$$F_j=\sum _{i<j}\frac{q_i{q}_j}{(i-j{)}^2}-\sum _{i>j}\frac{q_i{q}_j}{(i-j{)}^2}$$令 $E_i=\frac{F_i}{q_i}$，求 $E_i$。

题解

首先需要牢记卷积的形式：$c_k=\sum _{i=1}^k{a}_i{b}_{k-i}$，即满足两个数组的下标和为 $\text{sigma}$ 的上标。
将 $F_i$ 展开得
$$E_j=\sum _{i<j}\frac{q_i}{(i-j{)}^2}-\sum _{i>j}\frac{q_i}{(i-j{)}^2}$$

更改枚举顺序
$$E_j=\sum _{i=1}^{j-1}\frac{q_i}{(i-j{)}^2}-\sum _{i=j+1}^n\frac{q_i}{(i-j{)}^2}$$

令 $f_i=q_i$，$g_i=\frac{1}{i^2}$，特别地，$g_0=0$，则有：
$$E_j=\sum _{i=1}^{j-1}{f}_i{g}_{i-j}-\sum _{i=j+1}^n{f}_i{g}_{i-j}$$

又 $g_{i-j}=g_{j-i}$ ，则有：
$$A_j=\sum _{i=1}^{j-1}{f}_i{g}_{i-j}=\sum _{i=1}^{j-1}{f}_i{g}_{j-i}=\sum _{i=1}^j{f}_i{g}_{j-i}$$

就是一个标准的卷积。然后考虑$$B_j=\sum _{i=j+1}^n{f}_i{g}_{i-j}=\sum _{i=1}^{n-j}{f}_{i+j}{g}_i$$

令 $p_i=f_{n-i}$，则 $p_{n-i-j}=f_{i+j}$，则有：
$$B_j=\sum _{i=1}^{n-j}{f}_{i+j}{g}_i=\sum _{i=1}^{n-j}{p}_{n-i-j}{g}_i=B_{n-j}^{\prime }$$

则可以用FFT求出 $B^{\prime }$，然后根据 $B_i=B_{n-i}^{\prime }$ 计算。

代码

注意精度问题（别问我怎么知道的）：

1. 不要写 $1.0/(i\ast i)$，写$1.0/i/i$
2. 不要ans[i]/=limit，直接输出就好。

#include<bits/stdc++.h>
using namespace std;
const int maxn=862154;
const double pi=acos((double)-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);}
}f[maxn],g[maxn],p[maxn],A[maxn],B[maxn];
int limit=1,n,l=0,r[maxn];
void fft(comp *t,int ty){
    for(int i=0;i<limit;i++)
        if(i<r[i])
            swap(t[i],t[r[i]]);
    for(int mid=1;mid<limit;mid<<=1){
        comp wn(cos(pi/mid),ty*sin(pi/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=t[j+k],y=w*t[j+k+mid];
                t[j+k]=x+y;
                t[j+k+mid]=x-y;
            }
        }
    }
}
void reset(int n,int m){
    while(limit<=n+m+2)
        limit<<=1,++l;
    for(int i=1;i<limit;i++)
        r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
}
int main(void)
{
    memset(f,0,sizeof f);
    memset(g,0,sizeof g);
    scanf("%d",&n);
    for(int i=1;i<=n;i++){
        scanf("%lf",&f[i].x);
        g[i].x=(double)1.0/i/i;
        p[n-i].x=f[i].x;
    }
    reset(n,n);
    fft(f,1);fft(g,1);fft(p,1);
    for(int i=0;i<limit;i++)
        A[i]=g[i]*f[i];
    fft(A,-1);
    for(int i=0;i<limit;i++)
        B[i]=g[i]*p[i];
    fft(B,-1);
    for(int i=1;i<=n;i++)
        printf("%.6lf\n",A[i].x/limit-B[n-i].x/limit);
    return 0;
}