P3338 [ZJOI2014]力 [FFT]

P3338[ZJOI2014]力P3338 [ZJOI2014]力P3338[ZJOI2014]力

给出n个数qi，给出Fj的定义如下：

Fj=∑i<jqiqj(i−j)2−∑i>jqiqj(i−j)2F_j = \sum_{i<j}\frac{q_i q_j}{(i-j)^2 }-\sum_{i>j}\frac{q_i q_j}{(i-j)^2 }Fj=i<j∑(i−j)2qiqj−i>j∑(i−j)2qiqj

令Ei=Fi/qiEi=Fi/qiEi=Fi/qi，求EiEiEi.

Solution

Ei=Fi/qi=∑j<iqj(j−i)2−∑j>iqj(j−i)2Ei = F_i/q_i = \sum_{j<i}\frac{q_j}{(j-i)^2 }-\sum_{j>i}\frac{q_j}{(j-i)^2 }Ei=Fi/qi=j<i∑(j−i)2qj−j>i∑(j−i)2qj

设xi=i2x_i = i^2xi=i2
Ei=∑j<iqjxj−i−∑j>iqjxj−iEi = \sum_{j<i}q_jx_{j-i} - \sum_{j>i}q_jx_{j-i}Ei=j<i∑qjxj−i−j>i∑qjxj−i
即
Ei=∑j=1i−1qjxi−j−∑j=i+1nqjxi−jE_i = \sum_{j=1}^{i-1}q_jx_{i-j}-\sum_{j=i+1}^nq_jx_{i-j}Ei=j=1∑i−1qjxi−j−j=i+1∑nqjxi−j
这已经是卷积形式了
然而减数中的xi−jx_{i-j}xi−j下标为负数, 所以转换成
Ei=∑j=1i−1qjxi−j−∑j=i+1nqjxj−iE_i = \sum_{j=1}^{i-1}q_jx_{i-j}-\sum_{j=i+1}^nq_jx_{j-i}Ei=j=1∑i−1qjxi−j−j=i+1∑nqjxj−i
为了保持卷积形式, 继续设bi=qn−i+1b_i = q_{n-i+1}bi=qn−i+1
Ei=∑j=1i−1qjxi−j−∑j=i+1nbn−j+1xj−iE_i = \sum_{j=1}^{i-1}q_jx_{i-j}-\sum_{j=i+1}^nb_{n-j+1}x_{j-i}Ei=j=1∑i−1qjxi−j−j=i+1∑nbn−j+1xj−i
这又是一个卷积形式了!!!
于是FFTFFTFFT走起~

Attention

FFTFFTFFT要开2倍数组
不要漏掉蝴蝶变化

Code

#include<cstdio>
#include<cmath>
#include<algorithm>const int maxn = 1e6 + 5;
const double Pi = acos(-1);int N, N1, rev[maxn];struct complex{ double x, y; complex(double x = 0, double y = 0):x(x), y(y) {} } q[maxn], A[maxn], C[maxn], B[maxn];complex operator + (complex a, complex b){ return complex(a.x+b.x, a.y+b.y); }
complex operator - (complex a, complex b){ return complex(a.x-b.x, a.y-b.y); }
complex operator * (complex a, complex b){ return complex(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x); }void FFT(complex *F, short Opt){for(int i = 0; i < N1; i ++)if(i < rev[i]) std::swap(F[i], F[rev[i]]);for(int p = 2; p <= N1; p <<= 1){int len = p >> 1;complex tmp(cos(Pi/len), Opt*sin(Pi/len));for(int i = 0; i < N1; i += p){complex Buf(1, 0);for(int k = i; k < i+len; k ++){complex Temp = Buf * F[k+len];F[k+len] = F[k] - Temp;F[k] = F[k] + Temp;Buf = Buf * tmp;}}}
}int main(){scanf("%d", &N);for(int i = 1; i <= N; i ++) scanf("%lf", &q[i].x), A[N-i+1] = q[i];for(int i = 1; i <= N; i ++) C[i].x = (1.0/(double)(i))/(double)(i);N1 = 1;int bit_num = 0;while(N1 <= (N<<1)) N1 <<= 1, bit_num ++;for(int i = 0; i < N1; i ++) rev[i] = (rev[i>>1]>>1)|((i&1) << bit_num-1);FFT(A, 1), FFT(q, 1), FFT(C, 1);for(int i = 0; i <= N1; i ++) B[i] = q[i] * C[i], A[i] = A[i] * C[i];FFT(A, -1), FFT(B, -1);for(int i = 1; i <= N; i ++) printf("%.3lf\n", (B[i].x - A[N-i+1].x)/N1);return 0;
}