題目連結

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}\)，對於\(i\in[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]=\begin{cases}0 \quad i = 0 \\ \frac{1}{i^2} \quad i > 0,i\in[1,n] \\ \frac{-1}{i^2} \quad i < 0,i\in[-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 << '\n';
    return 0;
}