LuoguP3338「ZJOI2014」题解

发表于

基础FFT应用题

Description

给出$n$个数$q_1, q_2, \cdots, 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}$$

求$E_1, E_2, \cdots E_n$, $n \leq 10^5, q_i \leq 10^9$.

Samples

Sample Input #1

1
2
3
4
5
6
5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880

Sample Output #1

1
2
3
4
5
-16838672.693
3439.793
7509018.566
4595686.886
10903040.872

Solution

显然要推柿子，发现把$j$都加入两边的Sigma不影响答案，那么

$$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}$, 那么

$$\sum_{i=1}^j f_ig_{j-i} - \sum_{i=j}^n f_ig_{i-j}$$

定一下边界$f_0 = g_0 = 0$, 可以发现

$$\sum_{i=0}^j f_ig_{j-i} - \sum_{i=j}^n f_ig_{i - j}$$

哎， 左半边柿子已经变成了线性卷积的形式，这样就可以用$\texttt{FFT}$在$\mathcal O(n\log n)$的时间内计算了。

那么我们尝试把右半边的柿子也化成卷积形式，考虑翻转$f_i$, 设$h_i = f_{n - i}$，那么右半边

$$\sum_{i=j}^n f_ig_{i - j} = \sum_{i=0}^{n-j} f_{i+j}g_i$$

$$\sum_{i=0}^{n-j} h_{n - i - j}g_i$$

令$m = n - j$, 那么整个柿子就化成了两个都是卷积形式的部分

$$\sum_{i=0}^j f_ig_{j-i} - \sum_{i=0}^m g_ih_{m-i}$$

那么分别用$\texttt{FFT}$计算即可啦。

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
/* Headers */
#include <algorithm>
#include <iostream>
#include <cstring>
#include <climits>
#include <cstdio>
#include <cctype>
#include <vector>
#include <cmath>
#include <queue>
#include <stack>
#include <map>
#include <set>
#define FOR(i,a,b,c) for(int i=(a);i<=(b);i+=(c))
#define ROF(i,a,b,c) for(int i=(a);i>=(b);i-=(c))
#define FORL(i,a,b,c) for(long long i=(a);i<=(b);i+=(c))
#define ROFL(i,a,b,c) for(long long i=(a);i>=(b);i-=(c))
#define FORR(i,a,b,c) for(register int i=(a);i<=(b);i+=(c))
#define ROFR(i,a,b,c) for(register int i=(a);i>=(b);i-=(c))
#define RevEdge(x) x^1
#define lowbit(x) x&(-x)
#define LeftChild(x) x<<1
#define RightChild(x) (x<<1)+1
#define CLOSE_IN() fclose(stdin);
#define CLOSE_OUT() fclose(stdout);
#define FILE_IN(x) freopen(x,"r",stdin);
#define FILE_OUT(x) freopen(x,"w",stdout);
#define IOS(x) std::ios::sync_with_stdio(x)
#define Dividing() printf("-----------------------------------\n");
namespace FastIO {
#define gc() (iS == iT ? (iT = (iS = ibuff) + fread(ibuff, 1, SIZ, stdin), (iS == iT ? EOF : *iS++)) : *iS++)
    const int SIZ = 1 << 21 | 1;
    char* iS, * iT, ibuff[SIZ], obuff[SIZ], * oS = obuff, * oT = oS + SIZ - 1, fu[110], cc;
    int fr;
    inline void out() {
        fwrite(obuff, 1, oS - obuff, stdout);
        oS = obuff;
    }
    template<class Type>
    inline void read(Type& x) {
        x = 0;
        Type y = 1;
        for (cc = gc(); (cc > '9' || cc < '0') && cc != '-'; cc = gc());
        cc == '-' ? y = -1 : x = (cc & 15);
        for (cc = gc(); cc >= '0' && cc <= '9'; cc = gc())
            x = x * 10 + (cc & 15);
        x *= y;
    }
    template<class Type>
    inline void print(Type x, char text = '\n') {
        if (x < 0)
            * oS++ = '-', x *= -1;
        if (x == 0)
            * oS++ = '0';
        while (x)
            fu[++fr] = x % 10 + '0', x /= 10;
        while (fr)
            * oS++ = fu[fr--];
        * oS++ = text;
        out();
    }
    inline void prints(char x[], char text = '\n') {
        for (register int i = 0; x[i]; ++i)
            * oS++ = x[i];
        * oS++ = text;
        out();
    }
} using namespace FastIO;
template<typename T>
inline T max(T a, T b) {return (a > b) ? a : b;}
template<typename T>
inline T min(T a, T b) {return (a > b) ? b : a;}
/* definitions */
const int MAXN = 1e6 + 10;
const double pi = std::acos(-1.0);
struct Complex {
    double x, y;
    friend Complex operator + (const Complex &a, const Complex &b) {
        return ((Complex) {a.x + b.x, a.y + b.y});
    } friend Complex operator - (const Complex &a, const Complex &b) {
        return ((Complex) {a.x - b.x, a.y - b.y});
    } friend Complex operator * (const Complex &a, const Complex &b) {
        return ((Complex) {a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x});
    }
} f[MAXN], g[MAXN], h[MAXN];
int n, rev[MAXN], len = 1, lgm;
/* functions */
using std::cos;
using std::sin;
inline void FFT(Complex *f, int len, int type) {
    for(int i = 0; i < len; ++i) if(i < rev[i]) std::swap(f[i], f[rev[i]]);
    for(int i = 1; i < len; i <<= 1) {
        Complex Wn = (Complex) {cos(pi / (double)i), (double)type * sin(pi / (double)i)};
        for(int j = 0; j < len; j += (i << 1)) {
            Complex e = (Complex) {1, 0};
            for(int k = 0; k < i; ++k, e = e * Wn) {
                Complex x = f[j + k], y = e * f[i + j + k];
                f[j + k] = x + y; f[i + j + k] = x - y;
            }
        }
    } if(type == -1) for(int i = 0; i < len; ++i) f[i].x /= len;
}
int main(int argc, char *argv[]) {
    /*FILE_IN("data.in");*/ scanf("%d", &n);
    FOR(i, 1, n, 1) scanf("%lf", &f[i].x), h[n - i].x = f[i].x, g[i].x = (double)(1.0 / i / i);
    while(len <= (n << 1)) {len <<= 1; ++lgm;}
    FOR(i, 0, len - 1, 1) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (lgm - 1));
    FFT(f, len, 1); FFT(g, len, 1); FFT(h, len, 1);
    FOR(i, 0, len - 1, 1) f[i] = f[i] * g[i], h[i] = h[i] * g[i];
    FFT(f, len, -1); FFT(h, len, -1);
    FOR(i, 1, n, 1) printf("%.3lf\n", f[i].x - h[n - i].x);
    return 0;
}

THE END