Luogu P3338 [ZJOI2014]力

题目描述

https://www.luogu.com.cn/problem/P3338

Solution

首先我们令 $f_k=\frac{1}{k^2},g_k=q_k$，那么我们得到 $E_k=\sum_{i=0}^kf_{k-i}g_i-\sum_{i=k}^nf_{i-k}g_i$

左边显然是一个卷积，我们只考虑右边

$$\begin{aligned} \sum_{i=k}^nf_{i-k}g_i&=\sum_{i=0}^{n-k}f_{i}g_{i+k}\\ &=\sum_{i=0}^{n-k}f_{n-k-i}g_{n-i}\\ &=\sum_{i=0}^{n-k}f_{n-k-i}g'_{i} \end{aligned}$$

其中 $g’_k=g_{n-k}$，能够发现右边也是一个卷积的式子，所以我们做两遍 $FFT$ 即可

#include <iostream>
#include <cstdio>
#include <cmath>
#include <iomanip>
#define maxn 2097152
using namespace std;

const double pi = 3.141592653589793;

int n;

double q[maxn];

struct Complex {
    double x, y;

    friend Complex operator + (const Complex &u, const Complex &v) { return { u.x + v.x, u.y + v.y }; }
    friend Complex operator - (const Complex &u, const Complex &v) { return { u.x - v.x, u.y - v.y }; }
    friend Complex operator * (const Complex &u, const Complex &v) {
        return { u.x * v.x - u.y * v.y, u.x * v.y + u.y * v.x };
    }
} a[maxn], b[maxn], A[maxn], B[maxn];

int P[maxn];
void init_P(int n) {
    int l = 0; while ((1 << l) < n) ++l;
    for (int i = 0; i < n; ++i) P[i] = (P[i >> 1] >> 1) | ((i & 1) << l - 1);
}

void FFT(Complex *a, int n, int type) {
    for (int i = 0; i < n; ++i) if (i < P[i]) swap(a[i], a[P[i]]);
    for (int i = 2, m = 1; i <= n; m = i, i *= 2) 
        for (int j = 0; j < n; j += i) {
            Complex wn = { cos(2 * pi / i), type * sin(2 * pi / i) }, w = { 1, 0 }; 
            for (int k = 0; k < m; ++k, w = w * wn) {
                Complex t1 = a[j + k], t2 = a[j + k + m] * w;
                a[j + k] = t1 + t2; a[j + k + m] = t1 - t2; 
            } 
        }
    if (type < 0) for (int i = 0; i < n; ++i) a[i].x /= n; 
}

void Mul(Complex *a, Complex *B, int n1, int n2) {
    int n = 1; while (n < n1 + n2 - 1) n <<= 1; init_P(n);
    for (int i = 0; i < n2; ++i) b[i] = B[i];
    fill(a + n1, a + n, Complex { 0, 0 }); fill(b + n2, b + n, Complex { 0, 0 });
    FFT(a, n, 1); FFT(b, n, 1);
    for (int i = 0; i < n; ++i) a[i] = a[i] * b[i];
    FFT(a, n, -1); 
}

double ans[maxn];
int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr); cout.tie(nullptr);

    cin >> n; ++n; 
    q[0] = 0; for (int i = 1; i < n; ++i) cin >> q[i], A[i].x = q[i];
    B[0] = { 0, 0 }; for (int i = 1; i < n; ++i) B[i].x = 1.0 / i / i;
    Mul(A, B, n, n);
    for (int i = 1; i < n; ++i) ans[i] += A[i].x;
    for (int i = 0; i < n; ++i) A[i].x = q[n - 1 - i], A[i].y = 0; 
    Mul(A, B, n, n);
    for (int i = 1; i < n; ++i) ans[i] -= A[n - 1 - i].x;
    for (int i = 1; i < n; ++i) cout << fixed << setprecision(3) << ans[i] << "\n";
    return 0;
}