P3338 [ZJOI2014]力 [FFT]

$$P3338[ZJOI2014]力$$

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

$$F_j=\sum _{i<j}\frac{q_i{q}_j}{(i-j{)}^2}-\sum _{i>j}\frac{q_i{q}_j}{(i-j{)}^2}$$

令$Ei=Fi/qi$，求$Ei$.

---

Solution

$$Ei=F_i/q_i=\sum _{j<i}\frac{q_j}{(j-i{)}^2}-\sum _{j>i}\frac{q_j}{(j-i{)}^2}$$

设$x_i=i^2$
$$Ei=\sum _{j<i}{q}_j{x}_{j-i}-\sum _{j>i}{q}_j{x}_{j-i}$$
即
$$E_i=\sum _{j=1}^{i-1}{q}_j{x}_{i-j}-\sum _{j=i+1}^n{q}_j{x}_{i-j}$$
这已经是卷积形式了
然而 减数 中的$x_{i-j}$下标为负数, 所以转换成
$$E_i=\sum _{j=1}^{i-1}{q}_j{x}_{i-j}-\sum _{j=i+1}^n{q}_j{x}_{j-i}$$
为了保持卷积形式, 继续设$b_i=q_{n-i+1}$
$$E_i=\sum _{j=1}^{i-1}{q}_j{x}_{i-j}-\sum _{j=i+1}^n{b}_{n-j+1}{x}_{j-i}$$
这又是一个卷积形式了!!!
于是$FFT$走起~

---

Attention

$FFT$要开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;
}