1. 程式人生 > >【Bzoj3527】【Luogu3338】[Zjoi2014]力(FFT)

【Bzoj3527】【Luogu3338】[Zjoi2014]力(FFT)

題面

Bzoj

Luogu

題解

先來頹柿子

$$ F_i=\sum_{j<i}\frac{q_iq_j}{(i-j)^2}-\sum_{j>i}\frac{q_iq_j}{(i-j)^2} \\=q_i(\sum_{j<i}\frac{q_j}{(i-j)^2}-\sum_{j>i}\frac{q_j}{(i-j)^2}) $$

所以

$$ E_i=\sum_{j<i}\frac{q_j}{(i-j)^2}-\sum_{j>i}\frac{q_j}{(j-i)^2} $$

令$x_i=\frac{1}{i^2}$

則原式可化為:

$$ E_i=\sum_{j=0}^{i}q_jx_{i-j}-\sum_{j=i+1}^{n-1}q_jx_{j-i} $$

令$p_{n-j-1}=q_j$所以原式又等於

$$ E_i=\sum_{j=0}^{i}q_jx_{i-j}-\sum_{j=i}^{n-1}p_{n-j-1}x_{j-i} \\=\sum_{j=0}^{i}q_jx_{i-j}-\sum_{j=0}^{n-i-1}p_{n-j-i-1}x_{j} $$

然後就變成兩個卷積相減了吧。(為了方便,先將$n$減去$1$):

#include <cmath>
#include <cstdio>
#include <algorithm>
using std::swap;

const int N = 3e5 + 10;
const double Pi = acos(-1);
int n, m, r[N], P;
struct C { double x, y; } q[N], p[N], x[N];
C operator + (C a, C b) { return (C){ a.x + b.x, a.y + b.y }; }
C operator - (C a, C b) { return (C){ a.x - b.x, a.y - b.y }; }
C operator * (C a, C b) { return (C){ a.x * b.x - a.y * b.y, a.x * b.y + b.x * a.y }; }

void FFT(C f[], int opt) {
    for(int i = 0; i < n; ++i) if(i < r[i]) swap(f[i], f[r[i]]);
    for(int len = 1, nl = 2; len < n; len = nl, nl <<= 1) {
        C rot = (C){cos(Pi / len), opt * sin(Pi / len)};
        for(int l = 0; l < n; l += nl) {
            C w = (C){1, 0}; int r = l + len;
            for(int k = l; k < r; ++k, w = w * rot) {
                C x = f[k], y = w * f[k + len];
                f[k] = x + y, f[k + len] = x - y;
            }
        }
    }
} 

int main() {
    scanf("%d", &n); int tmp = (--n); m = n << 1;
    for(int i = 0; i <= n; ++i) scanf("%lf", &q[i].x), p[n - i].x = q[i].x;
    for(int i = 1; i <= n; ++i) x[i].x = (double)1.0 / i / i;
    for(n = 1; n <= m; n <<= 1, ++P);
    for(int i = 0; i < n; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (P - 1));
    FFT(q, 1), FFT(p, 1), FFT(x, 1);
    for(int i = 0; i < n; ++i) q[i] = q[i] * x[i], p[i] = p[i] * x[i];
    FFT(q, -1), FFT(p, -1);
    for(int i = 0; i <= tmp; ++i)
        printf("%.3lf\n", q[i].x / n - p[tmp - i].x / n);
    return 0;
}