[ZJOI2014] 力

Posted 2020-11-26 hock

题目大意

给定长度为(n)的序列(q_i)
[F_j = sum_{i = 1}^{j - 1} frac{q_i * q_j}{(i - j) ^ 2} - sum_{i = j+1}^{n} frac{q_i * q_j}{(i - j) ^ 2} ]
[E_i = frac {F_i}{q_i}]
对于 (1 leq i leq n) 求(E_i)

精度为(10^{-2})

正解

很明显的变形 (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}))
相应的 (E_i = sum_{i = 1}^{j - 1} frac{q_i}{(i - j) ^ 2} - sum_{i = j+1}^{n} frac{q_i}{(i - j) ^ 2})
我们尝试用卷积的形式表示(E_i)
令(g(x) = sum_{i = 1}^{i=n} frac {1}{i*i} * x ^ i) (f(x) = sum_{i=1}^{n} q_i*x^i)
那么(E_i = sum_{j=1}^{i-1}f(j) * g(i-j) - sum_{j=i+1}^{n}f(j)*g(j-i))
前面的式子是个很裸的卷积 后面式子 我们也尝试将它变为卷积
令(h(j) = f(n-j+1))
那么(E_i = sum_{j=1}^{i-1}f(j) * g(i-j) - sum_{j=i+1}^{n}f(n-j+1)*g(j-i))
两边都是卷积 (FFT)就行了 (4)次(FFT)被卡精度了 只好写(5)次了 哭唧唧

Code

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <iostream>
#include <set>
#include <map>
#include <queue>

using namespace std;

template <typename T> void chkmax(T &x, T y) {x = x > y ? x : y;}
template <typename T> void chkmin(T &x, T y) {x = x > y ? y : x;}

typedef long long ll;

const int INF = 2139062143;

#define DEBUG(x) std::cerr << #x << " = " << x << std::endl

template <typename T> void read (T &x) {
    x = 0; bool f = 1; char ch;
    do {ch = getchar(); if (ch == '-') f = 0;} while (ch > '9' || ch < '0');
    do {x = x * 10 + ch - '0'; ch = getchar();} while (ch >= '0' && ch <= '9');
    x = f ? x : -x;
}   

template <typename T> void write (T x) {
    if (x < 0) x = ~x + 1, putchar ('-');
    if (x > 9) write (x / 10);
    putchar (x % 10 + '0');
}

const long double PI = acos (-1.0);
const int N = 262144 + 5;


struct complex {
    long double x, y;
    complex (long double xx = 0.0, long double yy = 0.0) {x = xx, y = yy;}
} a[N], b[N], c[N];

complex operator * (complex a, complex b) {return complex (a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);}
complex operator + (complex a, complex b) {return complex (a.x + b.x, a.y + b.y);}
complex operator - (complex a, complex b) {return complex (a.x - b.x, a.y - b.y);}
complex operator / (complex a, long double b) {return complex ((long double)a.x / b, (long double)a.y / b);}

int n, r[N];

inline void FFT (complex *a, int limit, int opt) {
    for (int i = 0; i < limit; i ++ ) if (i < r[i]) swap (a[i], a[r[i]]);
    for (int i = 1; i < limit; i <<= 1) {
        complex Wn = complex (cos (PI / i), sin (PI / i) * opt);
        for (int j = 0; j < limit; j += (i << 1)) {
            complex w = complex (1, 0);
            for (int k = 0; k < i; k ++ , w = w * Wn) {
                complex p = a[j + k], q = w * a[j + k + i];
                a[j + k] = p + q; a[j + k + i] = p - q;
            }
        }
    }
    if (opt == -1) for (int i = 0; i < limit; i ++ ) a[i] = a[i] / (long double)limit;
}

inline void File () {
    read (n);
    for (int i = 0; i < n; i ++ ) scanf ("%Lf", &a[i].x), b[n - i - 1].x = a[i].x;
    for (int i = 1; i < n; i ++ ) c[i].x = (long double) 1 / i / i;

}

inline void Poly () {
    int limit = 1, l = 0; while (limit <= (n << 1)) limit <<= 1, l ++ ;
    for (int i = 0; i < limit; i ++ ) r[i] = ((r[i >> 1] >> 1) | ((i & 1) << (l - 1)));
    FFT (a, limit, 1); FFT (b, limit, 1); FFT (c, limit, 1);
    for (int i = 0; i < limit; i ++ ) a[i] = a[i] * c[i], b[i] = b[i] * c[i];
    FFT (a, limit, -1); FFT (b, limit, -1);
}

inline void Print () {
    for (int i = 0; i < n; i ++ ) printf ("%.10Lf
", (long double) a[i].x - b[n - i - 1].x);
}

int main() {
    File ();
    Poly ();
    Print ();
    return 0;
}