题意
E j = ∑ i = 1 j − 1 q i ( i − j ) 2 − ∑ i = j + 1 n 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} Ej=∑i=1j−1(i−j)2qi−∑i=j+1n(i−j)2qi
给你 q i q_i qi,求 E i E_i Ei。
分析
E j = ∑ i = 1 j − 1 q i ( i − j ) 2 − ∑ i = j + 1 n 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} Ej=∑i=1j−1(i−j)2qi−∑i=j+1n(i−j)2qi
令 f [ i ] = q i , g [ i ] = 1 i 2 f[i]=q_i,g[i]=\frac{1}{i^2} f[i]=qi,g[i]=i21
E j = ∑ i = 1 j f [ i ] g [ j − i ] − ∑ i = j n f [ i ] g [ i − j ] E_j=\sum_{i=1}^{j}f[i]g[j-i]-\sum_{i=j}^{n}f[i]g[i-j] Ej=∑i=1jf[i]g[j−i]−∑i=jnf[i]g[i−j]
前一项已经为卷积的形式,化简第二项
∑ 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=j+1}^{n}f[i]g[i-j]=f[j]g[0]+f[j+1]g[1]+...+f[j+(n-j)]g[n-j] ∑i=j+1nf[i]g[i−j]=f[j]g[0]+f[j+1]g[1]+...+f[j+(n−j)]g[n−j]
即 ∑ i = 0 n − j f [ j + i ] g [ i ] \sum_{i=0}^{n-j}f[j+i]g[i] ∑i=0n−jf[j+i]g[i]
令 f ′ [ i ] = f [ n − i ] f'[i]=f[n-i] f′[i]=f[n−i]
∑ i = 0 n − j f [ j + i ] g [ i ] = ∑ i = 0 n − j f [ n − j − i ] g [ i ] = ∑ i = 0 t f [ t − i ] g [ 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] ∑i=0n−jf[j+i]g[i]=∑i=0n−jf[n−j−i]g[i]=∑i=0tf[t−i]g[i]
故答案就行两项卷积相减,使用 f f t fft 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;
}