[题解] ZJOI2014 力

题目描述

给出 $n$ 个数 $q_i$ ,给出 $F_j$ 的定义如下:

$$F_j = \sum_{i<j} \frac{q_i q_j}{{(i-j)}^2} – \sum_{i>j} \frac{q_i q_j}{{(i-j)}^2} $$

令 $E_i = \frac{F_i}{q_i} $,求 $E_i$

输入格式

第一行一个整数 $n$
接下来 $n$ 行每行输入一个整数,第 $i$ 行表示 $q_i$

输出格式

$n$ 行,第 $i$ 行输出 $E_i$

样例输入

5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880

样例输出

-16838672.693
3439.793
7509018.566
4595686.886
10903040.872

数据范围及说明

时间限制 : $3s$
空间限制 : $128MB$

对于 $30 \%$ 的数据,$n \leq 1000$
对于 $50 \%$ 的数据,$n \leq 60000$
对于 $100 \%$ 的数据,$n \leq 100000, 0 < q_i < 1000000000 $

有 $spj$,误差不超过 $0.01$ 即可

解题思路

$$
\begin{aligned}
E_i &= \frac{F_i}{q_i} \\
&= \frac{\sum_{j < i} \frac{q_iq_j}{{(i-j)}^2} – \sum_{j>i} \frac{q_iq_j}{{(i-j)}^2}}{q_i} \\
&= \sum_{j < i} \frac{q_j}{{(i-j)}^2} – \sum_{j>i} \frac{q_j}{{(i-j)}^2} \\
&= \sum_{j=1}^{i-1} \frac{q_j}{{(i-j)}^2} – \sum_{j=i+1}^{n} \frac{q_j}{{(i-j)}^2} \\
\end{aligned}
$$

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

$$
\begin{aligned}
E_i &= \sum_{j=1}^{i-1} f[j]g[i-j] – \sum_{j=i+1}^{n} f[j]g[j-i] \\
&= \sum_{j=1}^{i-1} f[j]g[i-j] – \sum_{j=1}^{n-i} f[j+i]g[j] \\
&= \sum_{j=0}^{i-1} f[j]g[i-j] – \sum_{j=0}^{n-i} f[j+i]g[j] \\
\end{aligned}
$$

翻转 $f$ 数组得到 $f’$

$$
\begin{aligned}
E_i &= \sum_{j=0}^{i-1} f[j] g[i-j] – \sum_{j=0}^{n-i} f'[n-i-j+1]g[j] \\
\end{aligned}
$$

式子转化为两个卷积的形式,$FFT$ 即可。

代码

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>

using namespace std;
const int MAXN = 262144 + 10;
const long double PI = acos(-1.0);

struct complex{
    long double x, y;

    complex(long double _x = 0, long double _y = 0){
        x = _x, y = _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);
}

complex operator * (complex a, complex b){
    return complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
}

complex omega[MAXN], omegaInverse[MAXN];
int place[MAXN];

inline void init(int n){
    int k = 0;

    while((1 << k) < n)
        k++;

    for(register int i=0; i<n; i++){
        place[i] = (place[i >> 1] >> 1) | ((i & 1) << (k - 1));
    }

    long double single = 2 * PI / n;

    for(register int i=0; i<n; i++){
        omega[i] = complex(cos(single * i), sin(single * i));
        omegaInverse[i] = complex(omega[i].x, -omega[i].y);
    }
}

inline void transform(complex *a, int n, complex *omega){
    for(register int i=0; i<n; i++){
        if(i < place[i]){
            swap(a[i], a[place[i]]);
        }
    }

    for(register int range=2; range <= n; range <<= 1){
        int mid = range >> 1;
        int single = n / range;

        for(register complex *data=a; data != a + n; data += range){
            for(register int i=0; i<mid; i++){
                complex tmp = omega[single * i] * data[i + mid]; 

                data[i + mid] = data[i] - tmp;
                data[i] = data[i] + tmp;
            }
        }
    }
}

int n, m = 1;
long double f[MAXN], fInverse[MAXN], g[MAXN];
long double resultA[MAXN], resultB[MAXN];
complex tmp[MAXN], tmp2[MAXN];

inline void work(long double *a, long double *b, long double *result){
    for(register int i=0; i<m; i++){
        tmp[i].x = a[i];
        tmp2[i].x = b[i];

        tmp[i].y = 0;
        tmp2[i].y = 0;
    }

    transform(tmp, m, omega);
    transform(tmp2, m, omega);

    for(register int i=0; i<m; i++){
        tmp[i] = tmp[i] * tmp2[i];
    }

    transform(tmp, m, omegaInverse);

    for(register int i=0; i<n; i++){
        result[i] = tmp[i].x / m;
    }
}

int main(){

    scanf("%d", &n);

    for(register int i=1; i<=n; i++){
        scanf("%Lf", &f[i]);

        g[i] = (long double)1 / i / i;
        fInverse[i] = f[i];
    }

    n++;

    while(m < (n << 1)){
        m <<= 1;
    }

    init(m);
    reverse(fInverse + 1, fInverse + n);

    work(f, g, resultA);
    work(fInverse, g, resultB);

    n--;

    for(register int i=1; i<=n; i++){
        printf("%.10Lf\n", resultA[i] - resultB[n - i + 1]);
    }

    return 0;
}