题解 P3338 【[ZJOI2014]力】

题目链接

Solution [ZJOI2014]力

题目大意：给定长为(n)的数列(p)，(F_j=sum_{i=1}^{j-1}frac{q_iq_j}{(i-j)^2}-sum_{i=j+1}^{n}frac{q_iq_j}{(i-j)^2})，对于(iin[1,n])，求(E_i=frac{F_i}{q_i})

FFT、卷积

---

分析：

[F_j=sum_{i=1}^{j-1}frac{q_iq_j}{(i-j)^2}-sum_{i=j+1}^{n}frac{q_iq_j}{(i-j)^2} ]

首先最后要除(q_i)，那么我们直接约去

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

正着刚这个式子十分困难，我们可以从侧面，考虑每一个(q_i)会对哪些位置的答案产生怎样的贡献

不难发现，对于式子的前半部分，(q_i)会对(q_{i+k})产生(frac{q_i}{k^2})的贡献

把负号移进来，(q_i)会对(q_{i-k})产生(frac{-q_i}{k^2})的贡献

观察发现，这个东西与卷积十分相似

假设我们有

(A[i]=q_i)

(B[i]=egin{cases}0 quad i = 0 \ frac{1}{i^2} quad i > 0,iin[1,n] \ frac{-1}{i^2} quad i < 0,iin[-n,-1]end{cases})

那么(A,B)做个卷积就可以得到答案

下标没法为负数我们直接平移一下(B)就可以了

卷积直接(FFT)

#include <iostream>
#include <iomanip>
#include <cmath>
#include <algorithm>
using namespace std;
const int maxn = 4e5 + 100;
typedef double type;
const type pi = acos(-1);
struct complex{
    double x,y;
    complex operator + (const complex &rhs)const{return complex{x + rhs.x,y + rhs.y};}
    complex operator - (const complex &rhs)const{return complex{x - rhs.x,y - rhs.y};}
    complex operator * (const complex &rhs)const{return complex{x * rhs.x - y * rhs.y,x * rhs.y + y * rhs.x};}
}A[maxn << 1],B[maxn << 1];
type q[maxn];
int n,N,tr[maxn << 1];
inline void fft(complex *f,int flag){
    for(int i = 0;i < N;i++)
        if(i < tr[i])swap(f[i],f[tr[i]]);
    for(int p = 2;p <= N;p <<= 1){
        int len = p >> 1;
        const complex unit{cos(2 * pi / p),sin(2 * pi / p) * flag};
        for(int k = 0;k < N;k += p){
            complex g{1,0};
            for(int l = k;l < k + len;l++){
                complex t = f[l + len] * g;
                f[l + len] = f[l] - t;
                f[l] = f[l] + t;
                g = g * unit;
            }
        }
    }
}
int main(){
    ios::sync_with_stdio(false);
    cin >> n;
    for(int i = 1;i <= n;i++)cin >> A[i].x,q[i] = A[i].x;
    for(int i = 1;i <= n;i++){
        B[n + 1 - i].x = (type(-1) / i) / i;
        B[n + 1 + i].x = (type(1) / i) / i;
    }
    for(N = 1;N <= 3 * n + 1;N <<= 1);
    for(int i = 0;i < N;i++)
        tr[i] = (tr[i >> 1] >> 1) | ((i & 1) ? (N >> 1) : 0);
    fft(A,1),fft(B,1);
    for(int i = 0;i < N;i++)A[i] = A[i] * B[i];
    fft(A,-1);
    for(int i = 0;i < N;i++)A[i].x /= type(N);
    for(int i = n + 2;i <= 2 * n + 1;i++)cout << fixed << A[i].x << '
';
    return 0;
}