Luogu3338 [ZJOI2014]力

Posted 2021-02-07 xunzhen

题目蓝链

Description

给定一个长度为(n)的数列(q_i)，定义
[ 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)

Solution

直接上推导过程
[ E_j = frac{sum_{i < j} frac{q_i q_j}{(i - j)^2} - sum_{i > j} frac{q_i q_j}{(i - j)^2}}{q_j} = sum_{i < j} frac{q_i}{(i - j)^2} - sum_{i > j} frac{q_i}{(i - j)^2} ]
于是我们可以令
[ F(x) = q_x G(x) = egin{cases} -frac{1}{x^2}~~~(x < 0) ~~~~~~~~(x = 0) \frac{1}{x^2}~~~~~~(x > 0) end{cases} ]
则(下标默认从0开始)
[ E_j = sum_{i = 0}^{n - 1} F(j) cdot G(j - i) ]
所以这显然是一个卷积形式，所以我们就直接FFT求出卷积的系数即(E_i)即可

Code

#include <bits/stdc++.h>

using namespace std;

#define fst first
#define snd second
#define mp make_pair
#define squ(x) ((LL)(x) * (x))
#define debug(...) fprintf(stderr, __VA_ARGS__)

typedef long long LL;
typedef pair<int, int> pii;

template<typename T> inline bool chkmax(T &a, const T &b) { return a < b ? a = b, 1 : 0; }
template<typename T> inline bool chkmin(T &a, const T &b) { return a > b ? a = b, 1 : 0; }

inline int read() {
    int sum = 0, fg = 1; char c = getchar();
    for (; !isdigit(c); c = getchar()) if (c == ‘-‘) fg = -1;
    for (; isdigit(c); c = getchar()) sum = (sum << 3) + (sum << 1) + (c ^ 0x30);
    return fg * sum;
}

namespace FFT {

    const int MAX_LEN = 1 << 19;
    const double PI = acos(-1.0);

    struct com {
        double a, b;
        com (double _a = 0.0, double _b = 0.0): a(_a), b(_b) { }
        com operator + (const com &t) const { return com(a + t.a, b + t.b); }
        com operator - (const com &t) const { return com(a - t.a, b - t.b); }
        com operator * (const com &t) const { return com(a * t.a - b * t.b, a * t.b + b * t.a); }
    };

    int cnt, len, rev[MAX_LEN];
    com g[MAX_LEN];

    void init(int N) {
        for (cnt = -1, len = 1; len <= N; len <<= 1) ++cnt;
        for (int i = 0; i < len; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << cnt);
        g[0] = com(1.0, 0.0);
        com G(cos(PI * 2 / len), sin(PI * 2 / len));
        for (int i = 1; i < len; i++) g[i] = g[i - 1] * G;
    }

    void DFT(com *x, int op) {
        for (int i = 0; i < len; i++) if (i < rev[i]) swap(x[i], x[rev[i]]);
        for (int k = 2; k <= len; k <<= 1)
            for (int j = 0; j < len; j += k)
                for (int i = 0; i < k / 2; i++) {
                    com X = x[j + i], Y = x[j + i + k / 2] * g[len / k * (~op ? i : (i ? k - i : i))];
                    x[j + i] = X + Y, x[j + i + k / 2] = X - Y;
                }
        if (op < 0) for (int i = 0; i < len; i++) x[i].a /= len;
    }

    void mul(double *a, int n, double *b, int m, double *c) {
        init(n + m);
        static com F[MAX_LEN], G[MAX_LEN], S[MAX_LEN];
        for (int i = 0; i < len; i++) F[i] = com(i <= n ? a[i] : 0.0, 0.0);
        for (int i = 0; i < len; i++) G[i] = com(i <= m ? b[i] : 0.0, 0.0);
        DFT(F, 1), DFT(G, 1);
        for (int i = 0; i < len; i++) S[i] = F[i] * G[i];
        DFT(S, -1);
        for (int i = 0; i <= n + m; i++) c[i] = S[i].a;
    }

}

const int maxn = 3e5 + 10;

int n;
double f[maxn], g[maxn], e[maxn];

int main() {
#ifdef xunzhen
    freopen("force.in", "r", stdin);
    freopen("force.out", "w", stdout);
#endif

    n = read();
    for (int i = 0; i < n; i++) scanf("%lf", &f[i]);
    for (int i = 0; i < (n << 1); i++) if (i != n) g[i] = (i < n ? -1.0 : 1.0) / squ(i - n);

    FFT::mul(f, n, g, (n << 1) - 1, e);

    for (int i = n; i < (n << 1); i++) printf("%.3lf
", e[i]);

    return 0;
}