题目大意

给定$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 \leq i \leq n$, 求出 $E_i$

Solution

推下公式: $$ E_j = \frac{F_j}{q_j} = \sum_{i = 1}^{j - 1}{\frac{q_i}{(i - j)^2}} - \sum_{i = j + 1}^{n}{\frac{q_i}{(i - j) ^2}} $$ 然后如果要用FFT加速的话, 要把这个推成卷积式

卷积式就是$C_k = \sum_{0}^{k}{A[i] * B[k - i]}$

继续推 $$ E_j = \sum_{i = 1}^{j}{\frac{q_i}{(i - j)^2}} - \sum_{i = j}^{n}{\frac{q_i}{(i - j) ^ 2}} \\ $$ 设$f[i] = q_i$, $g[i] = \frac{1}{i ^ 2}$

原式变成 $$ E_j = \sum_{i = 1}^{j}{f[i] * g[j - i]} - \sum_{i = j}^{n}{f[i] * g[i - j]} $$ 然后设$f[0] = g[0] = 0$ $$ E_j = \sum_{i = 0}^{j}f[i] * g[j - i] - \sum_{i = j}^{n}f[i] * g[i - j] $$ 发现左边已经是一个卷积式子了, 继续推右边的.

将他展开, 发现 $$ \sum_{i = j}^{n}f[i] * g[i - j] = \\ f[j] * g[0] + f[j + 1] * g[1] + \dots + f[j + (n - j)] * g[n - j] \\ = \sum_{i = 0}^{n - j}f[j + i] * g[i] $$ 然后发现$f$还不太行, 翻转一下.

设$f'[i] = f[n - i]$ $$ \sum_{i = 0}^{n - j}f[j + i] * g[i] = \sum_{i = 0}^{n - j}f'[n - (j + i)] * g[i] \\ = \sum_{i = 0}^{n - j}f'[n - j - i] * g[i] \\ $$ 然后设$t = n - j$

原式 $$ \sum_{i = 0}^{t}f'[t - i] * g[i] $$ 答案就是 $$ E_j = \sum_{i = 0}^{j}{f[i] * g[j - i]} - \sum_{i = 0}^{t}f[t - i] * g[i]\quad (t = n - j) $$ 然后设$A(x) = \sum_{i = 0}^{n}f[i], B(x) = \sum_{i = 0}^{n} g[i] , C(x) = \sum_{i = 0}^{n}f'[i]$

设$D(x) = A(x) * B(x), E(x) = B(x) * C(x)$

然后 $E_j = D(j) - E(n - j)$

就做完了

代码:

#include <bits/stdc++.h>

using namespace std;

#define int long long

const long double PI = 3.14159265358978323846264338;
const int N = 1e6 + 10;

complex<double> omega[N], omegaInv[N];
int rev[N];
double f[N], g[N], ff[N], res1[N], res2[N];

inline void init(const int n)
{
    int k = 0;
    while ((1 << k) < n) k++;
    for (int i = 0; i < n; i++)
    {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (k - 1));
        omega[i] = complex<double>(cos(2 * PI / n * i), sin(2 * PI / n * i));
        omegaInv[i] = conj(omega[i]);
    }
    rev[0] = 0;
}

void fft(complex<double> *a, const int n, const complex<double> *omega)
{
    for (int i = 0; i < n; i++)
    {
        if (i > rev[i]) continue;
        swap(a[i], a[rev[i]]);
    }
    for (int l = 2; l <= n; l *= 2)
    {
        int m = l / 2;
        for (complex<double> *p = a; p != a + n; p += l)
            for (int i = 0; i < m; i++)
            {
                complex<double> t = omega[n / l * i] * p[m + i];
                p[m + i] = p[i] - t;
                p[i] += t;
            }
    }
}

inline void dft(complex<double> *a, const int n)
{
    fft(a, n, omega);
}

inline void idft(complex<double> *a, const int n)
{
    fft(a, n, omegaInv);
    for (int i = 0; i < n; i++)
    {
        a[i] /= n;
    }
}

inline void squared(int *a1, const int n1)
{
    int n = 1;
    while (n < n1 + n1 + 2) n *= 2;
    static complex<double> c1[N];
    for (int i = 0; i <= n1; i++) c1[i].real(a1[i]);
    init(n);
    dft(c1, n);
    for (int i = 0; i < n; i++)
        c1[i] *= c1[i];
    idft(c1, n);
    for (int i = 0; i <= n1 * 2; i++)
        a1[i] = static_cast<int>(floor(c1[i].real() + 0.5));
}

inline void multiply(double *a1, const int n1, double *a2, const int n2, double *res)
{
    int n = 1;
    while (n < n1 + n2 + 2) n *= 2;
    static complex<double> c1[N], c2[N];
    memset(c1, 0, sizeof(c1));
    memset(c2, 0, sizeof(c2));
    for (int i = 0; i <= n1; i++) c1[i].real(a1[i]);
    for (int i = 0; i <= n2; i++) c2[i].real(a2[i]);
    init(n);
    dft(c1, n);
    dft(c2, n);
    for (int i = 0; i < n; i++)
        c1[i] *= c2[i];
    idft(c1, n);
    // idft(c2, n);
    for (int i = 0; i <= n1 + n2; i++)
        res[i] = c1[i].real();
}

signed main()
{
    int n;
    cin >> n;
    for (int i = 1; i <= n; i++)
    {
        cin >> f[i];
        ff[n - i] = f[i];
        g[i] = (double)(1.0 / i / i);
    }
    // for (int i = 1; i <= n; i++)
    //     cin >> g[i];
    multiply(f, n, g, n, res1);
    // for (int i = 0; i <= n * 2; i++) 
    //     cout << fixed << setprecision(10) << res1[i] << ' ';
    // cout << endl;
    // exit(0);
    multiply(ff, n, g, n, res2);
    for (int i = 1; i <= n; i++)
        cout << fixed << setprecision(3) << res1[i] - res2[n - i] << endl;
    return 0;
}