BZOJ3527:「ZJOI2014」力

主要思路

首先知道:

\[E_i = \sum_{j = 1}^{i – 1}\frac{q_j}{(i – j)^2} – \sum_{j = i + 1}^{n}\frac{q_j}{(i – j)^2} \]

现在要求所有的。我们可以考虑分成前后两个部分分别计算。先考虑前面的部分:

\[a_i = \sum_{j = 1}^{i – 1}\frac{q_j}{(i – j)^2}\]

我们发现这一部分还是很好计算的:当多项式指数为\(i\)时,对应\(a_i\),发现\(q_j\)和\((i – j)^2\)的下标相加恰好为\(i\)。正好是卷积的形式。后面那一个部分也可以直接这样算。

代码

// BZ3527.cpp
#include <bits/stdc++.h>

using namespace std;

const int MAX_N = 400100;
const double pi = acos(-1.0);
typedef complex<double> cd;

int n, rev[MAX_N], poly_bit, poly_siz;
double qi[MAX_N], ans[MAX_N];
cd A[MAX_N], B[MAX_N], C[MAX_N], D[MAX_N];

void fft_initialize()
{
    for (int i = 0; i < poly_siz; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (poly_bit - 1));
}

void fft(cd *arr, int dft)
{
    for (int i = 0; i < poly_siz; i++)
        if (i < rev[i])
            swap(arr[i], arr[rev[i]]);
    for (int step = 1; step < poly_siz; step <<= 1)
    {
        cd omega_n = exp(cd(0, dft * pi / step));
        for (int j = 0; j < poly_siz; j += (step << 1))
        {
            cd omega_nk(1, 0);
            for (int k = j; k < j + step; k++, omega_nk *= omega_n)
            {
                cd t = arr[k + step] * omega_nk;
                arr[k + step] = arr[k] - t;
                arr[k] += t;
            }
        }
    }
    if (dft == -1)
        for (int i = 0; i < poly_siz; i++)
            arr[i] /= poly_siz;
}

double pow2(double x) { return x * x; }

int main()
{
    scanf("%d", &n);
    for (int i = 1; i <= n; i++)
        scanf("%lf", &qi[i]);

    int len = 2 * n - 1;
    while ((1 << poly_bit) < len)
        poly_bit++;
    poly_siz = (1 << poly_bit);
    fft_initialize();
    for (int i = 0; i < n; i++)
        A[i] = qi[i + 1], B[i] = 1.0 / pow2(i + 1);
    fft(A, 1), fft(B, 1);
    for (int i = 0; i < poly_siz; i++)
        A[i] = A[i] * B[i];
    fft(A, -1);
    for (int i = 2; i <= n; i++)
        ans[i] += A[i - 2].real();
    for (int i = 0; i <= poly_siz; i++)
        A[i] = B[i] = cd(0, 0);
    for (int i = 0; i < n; i++)
        A[i] = qi[n - i], B[i] = 1.0 / pow2(i + 1);
    fft(A, 1), fft(B, 1);
    for (int i = 0; i < poly_siz; i++)
        A[i] = A[i] * B[i];
    fft(A, -1);
    for (int i = n - 1; i >= 1; i--)
        ans[i] -= A[n - i - 1].real();
    for (int i = 1; i <= n; i++)
        printf("%.3lf\n", ans[i]);
    return 0;
}

 

Leave a Reply

Your email address will not be published. Required fields are marked *