「SDOI2017」相关分析

#题面

#题目描述

Frank 对天文学非常感兴趣，他经常用望远镜看星星， 同时记录下它们的信息，比如亮度、颜色等等，进而估算出星星的距离、半径等等。

Frank 不仅喜欢观测，还喜欢分析观测到的数据。他经常分析两个参数之间 （比如亮度和半径） 是否存在某种关系。

现在 Frank 要分析参数 $X$ 与 $Y$ 之间的关系。他有 $n$ 组观测数据，第 $i$ 组观测数据记录了 $x_i$ 和 $y_i$。他需耍进行以下几种操作：

• $\texttt{1 L R}$

  用直线拟合第 $L$ 组到第 $R$ 组观测数据。用 $\bar x$ 表示这些观测数据中 $x$ 的平均数，用 $\bar y$ 表示这些观测数据中 $y$ 的平均数，即

  $$\begin{aligned} \bar{x} &= \frac{1}{R - L + 1} \sum_{i = L} ^ R x_i \\ \bar{y} &= \frac{1}{R - L + 1} \sum_{i = L} ^ R y_i \end{aligned}$$

  如果直线方程是 $y = ax + b$，那么 $a$、$b$ 应该这样计算：

  $$\begin{aligned} a &= \frac{\sum_{i = L} ^ R (x_i - \bar{x})(y_i - \bar{y})}{\sum_{i = L} ^ R (x_i - \bar{x}) ^ 2} \\ b &= \bar{y} - a \bar{x} \end{aligned}$$

  你需要帮助 Frank 计算 $a$。

• $\texttt{2 L R S T}$

  Frank 发现测量第 $L$ 组到第 $R$ 组数据时有误差，对于每个 $i$ 满足 $L \leq i \leq R$，$x_i$ 需要加上 $S$，$y_i$ 需要加上 $T$。

• $\texttt{3 L R S T}$

  Frank 发现第 $L$ 组到第 $R$ 组数据需要修改，对于每个 $i$ 满足 $L \leq i \leq R$，$x_i$ 需要修改为 $(S + i)$，$y_i$ 需要修改为 $(T + i)$。

#输入格式

第一行两个数 $n$、$m$，表示观测数据组数和操作次数。

接下来一行 $n$ 个数，第 $i$ 个数是 $x_i$。

接下来一行 $n$ 个数，第 $i$ 个数是 $y_i$。

接下来 $m$ 行，表示操作，格式见题目描述。

#输出格式

对于每个 1 操作，输出一行，表示直线斜率 $a$。选手输出与标准输出的绝对误差或相对误差不超过 $10^{-5}$ 即为正确。

#样例输入输出

#样例输入 #1

3 5
1 2 3
1 2 3
1 1 3
2 2 3 -3 2
1 1 2
3 1 2 2 1
1 1 3

#样例输出 #1

1.0000000000
-1.5000000000
-0.6153846154

#数据范围与约定

• 对于 $20\%$ 的数据，$1 \leq n, m \leq 1000$；
• 对于另外 $20\%$ 的数据没有 3 操作，且 2 操作中 $S = 0$；
• 对于另外 $30\%$ 的数据没有 3 操作；
• 对于 $100\%$ 的数据，$1 \leq n, m \leq 10^5$，$1 \leq L \leq R \leq n$，$0 \leq |S|, |T| \leq 10^5$，$0 \leq |x_i|, |y_i| \leq 10^5$，1 操作中不会出现分母为 $0$ 这类特殊情况。

#思路

首先需要转化 $a$ 式。

$$\begin{aligned} a &= \frac{\sum_{i = L}^{R} (x_i - \bar{x}) \times (y_i - \bar{y})}{\sum_{i = L}^{R} (x_i - \bar{x})^2} \\ &= \frac{\sum_{i = L}^{R} (x_i \times y_i - x_i \times \bar{y} - y_i \times \bar{x} + \bar{x} \times \bar{y})}{\sum_{i = L}^{R} x_i^2 - 2 x_i \bar{x} + \bar{x}^2} \\ &= \frac{(\sum_{i = L}^{R} x_i y_i) - (\bar{y} \times \sum_{i = L}^{R} x_i) - (\bar{x} \times \sum_{i = L}^{R} y_i) + \Big((R - L + 1) \times \bar{x} \times \bar{y} \Big)}{(\sum_{i = L}^{R} x_i^2) - (2 \bar{x} \sum_{i = L}^{R} x_i) + \Big((R - L + 1) \bar{x}^2 \Big)} \\ &= \frac{\sum_{i = L}^{R} x_i y_i - \frac{\sum_{i = L}^{R} x_i \sum_{i = L}^{R} y_i}{R - L + 1} - \frac{\sum_{i = L}^{R} x_i \sum_{i = L}^{R} y_i}{R - L + 1} + \frac{\sum_{i = L}^{R} x_i \sum_{i = L}^{R} y_i}{R - L + 1}}{\sum_{i = L}^{R} x_i^2 - 2 \times \frac{(\sum_{i = L}^{R} x_i) ^ 2}{R - L + 1} + \frac{(\sum_{i = L}^{R} x_i) ^ 2}{R - L + 1}} \\ &= \frac{\sum_{i = L}^{R} x_i y_i - \frac{\sum_{i = L}^{R} x_i \sum_{i = L}^{R} y_i}{R - L + 1}}{\sum_{i = L}^{R} x_i^2 - \frac{(\sum_{i = L}^{R} x_i) ^ 2}{R - L + 1}} \end{aligned}$$

这样只需要维护 $\sum x_i$、$\sum y_i$、$\sum x_i y_i$、$\sum x_i^2$ 即可，根据上式即可计算出斜率 $a$，免去了计算平均数的麻烦。

显然可以用线段树维护这些值。

1. 对于操作 1，求出区间和后代入上式即可计算出答案。

2. 对于操作 2，与线段树区间加类似，只不过处理懒标记时比较麻烦：

   • $\displaystyle \sum_{i = L}^{R} (x_i + S) = \sum_{i = L}^{R} x_i + S \times (R - L + 1)$；
   • $\displaystyle \sum_{i = L}^{R} (y_i + T) = \sum_{i = L}^{R} y_i + T \times (R - L + 1)$；
   • $\displaystyle \sum_{i = L}^{R} (x_i + S) \times (y_i + T) = \sum_{i = L}^{R} (x_i y_i + T x_i + S y_i + ST) = \sum_{i = L}^{R} x_i y_i + T \sum_{i = L}^{R} x_i + S \sum_{i = L}^{R} y_i + S \times T \times (R - L + 1)$；
   • $\displaystyle \sum_{i = L}^{R} (x_i + S)^2 = \sum_{i = L}^{R} (x_i^2 + 2 x_i S + S^2) = \sum_{i = L}^{R} x_i^2 + 2 \times S \times \sum_{i = L}^{R} x_i + S^2 \times (R - L + 1)$。

3. 对于操作 3，可以先进行区间覆盖，再进行一次操作 2 将 $S, T$ 添加到区间上。

   在区间覆盖时，利用公式 $\displaystyle \sum_{i = 1}^{n} i^2 = \frac{n \cdot (n + 1) \cdot (2n + 1)}{6}$ 可以快速计算出区间平方和。

本题数据会爆 long long，请开 __int128。

#代码

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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202

#include <iostream>
#include <iomanip>

using std::cin;
using std::cout;
const char endl = '\n';

const int N = 1e5 + 5;

struct node {
    int l, r;
    __int128 sum_x, sum_y, sum_xx, sum_xy;
    __int128 lazy_s, lazy_t;
    bool covered;

    node(const int &_l = 0, const int &_r = 0)
        : l(_l),
          r(_r),
          sum_x(0),
          sum_y(0),
          sum_xx(0),
          sum_xy(0),
          lazy_s(0),
          lazy_t(0),
          covered(false) {}

    node operator+(const node &rhs) const {
        node res;

        res.l = l;
        res.r = rhs.r;
        res.sum_x = sum_x + rhs.sum_x;
        res.sum_y = sum_y + rhs.sum_y;
        res.sum_xx = sum_xx + rhs.sum_xx;
        res.sum_xy = sum_xy + rhs.sum_xy;

        return res;
    }
} tr[N << 2];

int n, m, x[N], y[N];

__int128 square_sum(__int128 x) {
    return x * (x + 1) * (2 * x + 1) / 6;
}

__int128 square_sum(__int128 l, __int128 r) {
    return square_sum(r) - square_sum(l - 1);
}

/// Add (s, t) to each element under tr[u]
void modify_node(int u, __int128 s, __int128 t) {
    tr[u].lazy_s += s;
    tr[u].lazy_t += t;

    tr[u].sum_xx += 2 * s * tr[u].sum_x + s * s * (tr[u].r - tr[u].l + 1);
    tr[u].sum_xy += t * tr[u].sum_x + s * tr[u].sum_y + s * t * (tr[u].r - tr[u].l + 1);

    tr[u].sum_x += s * (tr[u].r - tr[u].l + 1);
    tr[u].sum_y += t * (tr[u].r - tr[u].l + 1);
}

/// x[i] = y[i] = i
void change_node(int u) {
    tr[u].lazy_s = tr[u].lazy_t = 0;
    tr[u].sum_x = tr[u].sum_y = static_cast<__int128>(tr[u].l + tr[u].r) * (tr[u].r - tr[u].l + 1) / 2;
    tr[u].sum_xx = tr[u].sum_xy = square_sum(tr[u].l, tr[u].r);
    tr[u].covered = true;
}

void pushup(int u) {
    tr[u] = tr[u << 1] + tr[u << 1 | 1];
}

void pushdown(int u) {
    if (tr[u].covered) {
        change_node(u << 1);
        change_node(u << 1 | 1);

        tr[u].covered = false;
    }

    modify_node(u << 1, tr[u].lazy_s, tr[u].lazy_t);
    modify_node(u << 1 | 1, tr[u].lazy_s, tr[u].lazy_t);

    tr[u].lazy_s = tr[u].lazy_t = 0;
}

void build(int u, int l, int r) {
    tr[u] = node(l, r);

    if (l == r) {
        tr[u].sum_x = x[l];
        tr[u].sum_y = y[l];
        tr[u].sum_xx = static_cast<__int128>(x[l]) * x[l];
        tr[u].sum_xy = static_cast<__int128>(x[l]) * y[l];

        return;
    }

    int mid = l + r >> 1;

    build(u << 1, l, mid);
    build(u << 1 | 1, mid + 1, r);

    pushup(u);
}

void modify(int u, int l, int r, __int128 s, __int128 t) {
    if (l <= tr[u].l && tr[u].r <= r) {
        modify_node(u, s, t);

        return;
    }

    int mid = tr[u].l + tr[u].r >> 1;

    pushdown(u);

    if (l <= mid) modify(u << 1, l, r, s, t);
    if (r > mid) modify(u << 1 | 1, l, r, s, t);

    pushup(u);
}

void change(int u, int l, int r, __int128 s, __int128 t) {
    if (l <= tr[u].l && tr[u].r <= r) {
        change_node(u);
        modify_node(u, s, t);

        return;
    }

    int mid = tr[u].l + tr[u].r >> 1;

    pushdown(u);

    if (l <= mid) change(u << 1, l, r, s, t);
    if (r > mid) change(u << 1 | 1, l, r, s, t);

    pushup(u);
}

node query(int u, int l, int r) {
    if (l <= tr[u].l && tr[u].r <= r) return tr[u];

    int mid = tr[u].l + tr[u].r >> 1;
    node res;

    pushdown(u);

    if (l <= mid) res = res + query(u << 1, l, r);
    if (r > mid) res = res + query(u << 1 | 1, l, r);

    return res;
}

int main() {
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr);

    cin >> n >> m;

    for (int i = 1; i <= n; i++) {
        cin >> x[i];
    }

    for (int i = 1; i <= n; i++) {
        cin >> y[i];
    }

    build(1, 1, n);

    while (m--) {
        int op, l, r;

        cin >> op >> l >> r;

        if (op == 1) {
            auto res = query(1, l, r);

            double x = static_cast<double>(res.sum_xy) - static_cast<double>(res.sum_x) * res.sum_y / (r - l + 1),
                   y = static_cast<double>(res.sum_xx) - static_cast<double>(res.sum_x) * res.sum_x / (r - l + 1);

            cout << std::fixed << std::setprecision(10) << x / y << endl;
        } else if (op == 2) {
            int s, t;

            cin >> s >> t;

            modify(1, l, r, s, t);
        } else {  // op == 3
            int s, t;

            cin >> s >> t;

            change(1, l, r, s, t);
        }
    }

    return 0;
}