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_jx_{j-i} - sum_{j>i}q_jx_{j-i}$
即
$E_i = sum_{j=1}^{i-1}q_jx_{i-j}-sum_{j=i+1}^nq_jx_{i-j}$
这已经是卷积形式了
然而减数中的 $x_{i-j}$ 下标为负数, 所以转换成
$E_i = sum_{j=1}^{i-1}q_jx_{i-j}-sum_{j=i+1}^nq_jx_{j-i}$
为了保持卷积形式, 继续设 $b_i = q_{n-i+1}$
$E_i = sum_{j=1}^{i-1}q_jx_{i-j}-sum_{j=i+1}^nb_{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
", (B[i].x - A[N-i+1].x)/N1);
        return 0;
}