Tsr and Variance

题目描述

Tsr is a cute boy with handsome moustache.

You are given a sequence with length $n$. Tsr wants you to calculate the sum of variance of each successive subsequence. Note: The variance in this problem should’t be divided by length.

Recall $\overline{a}{l, r}=\frac{1}{r-l+1} \sum{i=l}^r a_i$. Then you are supposed to calculate $\sum_{i=1}^n \sum_{j=i}^n \sum_{k=i}^j (a_k-\overline{a}_{i, j})^2$.

题意概述

给定一个长度为$n$的序列,求$\sum_{i=1}^n \sum_{j=i}^n \sum_{k=i}^j (a_k-\overline{a}_{i, j})^2$。有$T$组数据。

数据范围:$1 \le T \le 20, ; 1 \le n \le 10000, ; 1 \le a_i \le 10$。

算法分析

$$\begin{align} \sum_{i=1}^n \sum_{j=i}^n \sum_{k=i}^j (a_k-\overline{a}{i, j})^2 &= \sum{i=1}^n \sum_{j=i}^n \sum_{k=i}^j (a_k^2-2a_k \overline{a}{i, j}+\overline{a}{i, j}^2) \\ &= \sum_{k=1}^n \sum_{i=1}^k \sum_{j=k}^n a_k^2-2\sum_{i=1}^n \sum_{j=i}^n \overline{a}{i, j} \sum{k=i}^j a_k+ \sum_{i=1}^n \sum_{j=i}^n (j-i+1) \overline{a}{i, j}^2 \\ &= \sum{k=1}^n k(n-k+1)a_k^2-\sum_{i=1}^n \sum_{j=i}^n (j-i+1) \overline{a}_{i, j}^2 \end{align}$$

前一项可以在$O(n)$时间内求出,只需考虑如何求后一项。令$s_i=\sum_{j=1}^i a_j$。

$$\begin{align} \sum_{i=1}^n \sum_{j=i}^n (j-i+1) \overline{a}{i, j}^2 &= \sum{j=1}^n \sum_{i=1}^j \frac{(s_j-s_{i-1})^2}{j-i+1} \\ &= \sum_{j=1}^n \sum_{i=1}^j \frac{s_j^2-2s_js_{i-1}+s_{i-1}^2}{j-i+1} \\ &= \sum_{j=1}^n s_j^2 \sum_{i=1}^j \frac{1}{i}-2\sum_{j=1}^n s_j \sum_{i=1}^j \frac{s_{i-1}}{j-(i-1)}+\sum_{j=1}^n \sum_{i=1}^j \frac{s_{i-1}^2}{j-(i-1)} \end{align}$$

第一项也可以在$O(n)$时间内求出,而后面两项都是卷积的形式,可以用 FFT 在$O(n\log n)$时间内求出。

代码

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
#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>

int const N = 10005, T = 400005;
double const PI = acos(-1);

int a[N], rev[T];

struct Point {
double x, y;
Point(double _x = 0, double _y = 0) : x(_x), y(_y) {}
friend Point const operator + (Point const &a, Point const &b) {
return Point(a.x + b.x, a.y + b.y);
}
friend Point const operator - (Point const &a, Point const &b) {
return Point(a.x - b.x, a.y - b.y);
}
friend Point const operator * (Point const &a, Point const &b) {
return Point(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
}
friend Point const operator / (Point const &a, double const &n) {
return Point(a.x / n, a.y / n);
}
} wn[T], A[T], B[T], C[T];

int init(int n) {
int m = n, l = 0;
for (n = 1; n <= m; n <<= 1, ++l) ;
for (int i = 1; i < n; ++i) {
rev[i] = rev[i >> 1] >> 1 | (i & 1) << l - 1;
}
for (int i = 0; i < n >> 1; ++i) {
wn[i] = Point(cos(2 * PI / n * i), sin(2 * PI / n * i));
}
return n;
}

void fft(Point *a, int n, int inv = 0) {
for (int i = 0; i < n; ++i) {
if (i < rev[i]) {
std::swap(a[i], a[rev[i]]);
}
}
for (int i = 1; i < n; i <<= 1) {
for (int j = 0; j < n; j += i << 1) {
for (int k = 0; k < i; ++k) {
Point x = a[j + k], y = wn[n / (i << 1) * k] * a[j + k + i];
a[j + k] = x + y;
a[j + k + i] = x - y;
}
}
}
if (inv) {
std::reverse(a + 1, a + n);
for (int i = 0; i < n; ++i) {
a[i] = a[i] / n;
}
}
}

int main() {
int T;
scanf("%d", &T);
for (; T--;) {
int n;
scanf("%d", &n);
double ans = 0;
for (int i = 1; i <= n; ++i) {
scanf("%d", &a[i]);
ans += 1. * a[i] * a[i] * i * (n - i + 1);
a[i] += a[i - 1];
}
double sum = 0;
for (int i = 1; i <= n; ++i) {
sum += 1. / i;
ans -= sum * a[i] * a[i];
}
int len = init(n << 1);
for (int i = 0; i < n; ++i) {
A[i] = Point(a[i]);
B[i] = Point(1. * a[i] * a[i]);
C[i] = Point(1. / (i + 1));
}
for (int i = n; i < len; ++i) {
A[i] = B[i] = C[i] = Point(0);
}
fft(A, len);
fft(B, len);
fft(C, len);
for (int i = 0; i < len; ++i) {
A[i] = A[i] * C[i];
B[i] = B[i] * C[i];
}
fft(A, len, 1);
fft(B, len, 1);
for (int i = 0; i < n; ++i) {
ans += 2 * a[i + 1] * A[i].x - B[i].x;
}
printf("%.10f\n", ans);
}
return 0;
}

本博客所有文章除特别声明外,均采用 CC BY-SA 4.0 协议 ,转载请注明出处!