P3338 [ZJOI2014]力 快速傅里叶变换fft，推式子

题意

$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$,求$E_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}$

令$f[i]=q_i,g[i]=\frac{1}{i^2}$

$E_j={\sum }_{i=1}^{j}f[i]g[j-i]-{\sum }_{i=j}^{n}f[i]g[i-j]$

前一项已经为卷积的形式，化简第二项

${\sum }_{i=j+1}^{n}f[i]g[i-j]=f[j]g[0]+f[j+1]g[1]+...+f[j+(n-j)]g[n-j]$

即${\sum }_{i=0}^{n-j}f[j+i]g[i]$

令$f^{\prime }[i]=f[n-i]$

${\sum }_{i=0}^{n-j}f[j+i]g[i]={\sum }_{i=0}^{n-j}f[n-j-i]g[i]={\sum }_{i=0}^{t}f[t-i]g[i]$

故答案就行两项卷积相减，使用$fft$加速可得到答案。

代码

#include <bits/stdc++.h>
using namespace std;

const int N = 1e6+10, M = N * 2, inf = 1e8;
const double PI = acos(-1);
int n, m;
int rev[N], bit, tot;
struct Complex{
    double x, y;
    Complex operator+ (const Complex& t) const { return {x + t.x, y + t.y}; }
    Complex operator- (const Complex& t) const{ return {x - t.x, y - t.y}; }
    Complex operator* (const Complex& t) const{ return {x * t.x - y * t.y, x * t.y + y * t.x}; }
}a[N], b[N], c[N];
// a[i] = j 相当于ja^i，a的i次幂前系数为j
void fft(Complex a[], int inv){
    for(int i = 0; i < tot; i++)
        if(i < rev[i]) swap(a[i], a[rev[i]]);
    for(int mid = 1; mid < tot; mid <<= 1){
        Complex w1 = Complex({cos(PI/mid), inv * sin(PI/mid)});
        for(int i = 0; i < tot; i += mid * 2){
            Complex wk = Complex({1, 0});
            for(int j = 0; j < mid; j ++, wk = wk * w1){
                Complex x = a[i + j], y = wk * a[i + j + mid];
                a[i + j] = x + y, a[i + j + mid] = x - y;
            }
        }
    }
}
int main(){
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    cin>>n;
    while((1 << bit) < n + n + 1) bit++; tot = 1 << bit;
    for(int i = 0; i < tot; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
    
    for(int i = 1; i <= n; i++) cin>>a[i].x, b[i].x = a[i].x, c[i].x = 1.0 / (double) i / (double) i;
    
    reverse(b + 1, b + n + 1);
    fft(a, 1), fft(b, 1); fft(c, 1);
    for(int i = 0; i < tot; i++) a[i] = a[i] * c[i];
    for(int i = 0; i < tot; i++) b[i] = b[i] * c[i];
    fft(a, -1); fft(b, -1);
    
    for(int i = 0; i <= n + m; i++) a[i].x = (a[i].x/tot),b[i].x = (b[i].x/tot);
    
    for(int i = 1; i <= n; i++) printf("%.6lf\n", a[i].x - b[n-i+1].x);
    return 0;
}