P3338 [ZJOI2014]力

P3338 [ZJOI2014]力

题解

题意

给出$n$个数$q_1,q_2, \dots, q_n$定义

$$F_j= \sum_{i = 1} ^ {j - 1} \frac{q_i \times q_j}{(i - j)^2} - \sum_{i = j + 1} ^ n \frac{q_i \times q_j}{(i - j) ^2} \\ E_i = \frac{F_i} {q_i}$$

对于$1 \le i \le n$, 求$E_i$的值。

$1 \le n \le 10 ^ 5$

解法

首先这个题的数据范围有很大的问题， 原题目中给出的数据范围为$0 < q_i <10^ 9$， 但是数据中存在$q_i = 0$的情况， 所以需要在计算时先将$q_i$消去，（很让人不爽）。

现在要求计算的就是

$$F_i = \sum_{j = 1} ^ {i - 1} \frac{q_j}{(i - j)^2} - \sum_{i = j + 1} ^ n \frac{q_j}{(i - j)^2}$$

首先考虑拆开计算，设为$F_i = A_i - B_i$， $A_i = \sum_{j = 1} ^ {j - 1} \frac{q_j}{(i - j)^2}$，观察可知，是一个卷积的形式，可以变为

$$A_i = \sum_{j + k = i}q_j \times \frac{1}{k ^2}$$

然后用FFT就可以计算出$A$了，然后考虑$B$的形式和$A$的一样， 可以直接将$q_i$,reverse一下，然后计算出$B$即可。

#include <bits/stdc++.h>

using namespace std;

const int N =4e5 + 10;
const double pi = acos(-1.0);

struct Complex
{
    double x, y;
    Complex(double a = 0, double b = 0) : x(a), y(b) {}
    friend Complex operator + (Complex a, Complex b) {return Complex(a.x + b.x, a.y + b.y);}
    friend Complex operator - (Complex a, Complex b) {return Complex(a.x - b.x, a.y - b.y);}
    friend Complex operator * (Complex a, Complex b) {return Complex(a.x * b.x - a.y * b.y, a.y * b.x + b.y * a.x);}
};

int n;

double q[N], a[N], b[N];

Complex F[N], G[N];

int rev[N], len = 1;

void FFT(Complex *a, int len, int type)
{
    for(int i = 0; i < len; i++)
        if(i < rev[i])swap(a[i], a[rev[i]]);
    for(int k = 1; k < len; k <<= 1)
    {
        Complex x(cos(pi / k), type * sin(pi / k));
        for(int i = 0; i < len; i += k << 1)
        {
            Complex w(1, 0);
            for(int j = 0; j < k; j++)
            {
                Complex y = a[i + j];
                Complex z = w * a[i + j + k];
                a[i + j] = y + z;
                a[i + j + k] = y - z;
                w = w * x;
            }
        }
    }
    if(type == -1)
        for(int i = 0; i < len; i++)
            a[i].x /= len;
}

void init()
{
    memset(F, 0, sizeof(F));
    memset(G, 0, sizeof(G));
}

void Get()
{
    for(int i = 1; i <= n; i++)
        F[i].x = q[i];
    for(int i = 1; i <= n; i++)
        G[i].x = 1.0 / i / i;
    FFT(F, len, 1), FFT(G, len, 1);
    for(int i = 0; i <= len; i++)
        F[i] = F[i] * G[i];
    FFT(F, len, -1);
}

int main()
{
    scanf("%d", &n);
    for(int i = 1; i <= n; i++)
        scanf("%lf", &q[i]);
    int cnt = 0;
    while(len <= (n << 1))len <<= 1, cnt++;
    for(int i = 0; i <= len; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (cnt - 1));
    init(); Get();
    for(int i = 1; i <= n; i++)
        a[i] = F[i].x;
    reverse(q + 1, q + n + 1);
    init(); Get();
    for(int i = 1; i <= n; i++)
        b[i] = F[i].x;
    reverse(q + 1, q + n + 1);
    for(int i = 1; i <= n; i++)
        printf("%.8lf\n", a[i] - b[n - i + 1]);
    return 0;
}