[BZOJ4821][SDOI2017]相关分析（线段树）

1、概述

一道比BZOJ1858还要恶心的线段树……
调了好几个小时……

2、维护内容

考虑拆开式子：
$\frac{\sum_{i=L}^R(x_i-\bar x)(y_i-\bar y)}{\sum_{i=L}^R(x_i-\bar x)^2}=\frac{\sum_{i=L}^Rx_iy_i-\frac{(\sum_{i=L}^Rx_i)*(\sum_{i=L}^Ry_i)}{R-L+1}}{\sum_{i=L}^Rx_i^2-\frac{(\sum_{i=L}^Rx_i)^2}{R-L+1}}$。
所以只要维护$\sum x,\sum y,\sum xy,\sum x^2$的值。
对于修改，维护$4$个标记：
$addX$：$x$区间加的标记。
$addY$：$y$区间加的标记。
$tagX$：$x$区间成段更新的标记，并且每一个$x$值加上$i$。
$tagY$：$y$区间成段更新的标记，并且每一个$y$值加上$i$。

3、标记更新

标记的更新一个比一个恶心……
以下$l$为左端点，$r$为右端点。
对于区间加操作，
$\sum xy+=T\sum x+S\sum y+(r-l+1)*S*T$
$\sum x^2+=2S\sum x+(r-l+1)*S^2$
$\sum x+=(r-l+1)*S$
$\sum y+=(r-l+1)*T$
对于区间更新操作，
$\sum x=(r-l+1)*S+\frac{(r-l+1)*(l+r)}{2}$
$\sum y=(r-l+1)*T+\frac{(r-l+1)*(l+r)}{2}$
$\sum x^2=(r-l+1)*S^2+S*(l+r)*(r-l+1)+\frac{r(r+1)(2r+1)}{6}-\frac{l(l-1)(2l-1)}{6}$
$\sum xy=(r-l+1)*S*T+\frac{(S+T)*(l+r)*(r-l+1)}{2}+\frac{r(r+1)(2r+1)}{6}-\frac{l(l-1)(2l-1)}{6}$

4、代码

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
inline int read() {
    int res = 0; bool bo = 0; char c;
    while (((c = getchar()) < '0' || c > '9') && c != '-');
    if (c == '-') bo = 1; else res = c - 48;
    while ((c = getchar()) >= '0' && c <= '9')
        res = (res << 3) + (res << 1) + (c - 48);
    return bo ? ~res + 1 : res;
}
const int N = 1e5 + 5, M = 4e5 + 5;
int n, m, X[N], Y[N], QAQ = 1;
struct cyx {
    double sumX, sumY, sumXY, sumX2, addX, addY,
        firX, firY; bool cn;
    int l, r, lc, rc, inc;
} T[M];
struct pyz {
    double sumX, sumY, sumXY, sumX2;
    friend inline pyz operator + (pyz a, pyz b) {
        pyz res; res.sumX = a.sumX + b.sumX;
        res.sumY = a.sumY + b.sumY;
        res.sumX2 = a.sumX2 + b.sumX2;
        res.sumXY = a.sumXY + b.sumXY;
        return res;
    }
};
inline double sqr(const double x) {return x * x;}
inline double sumx2(const double l, const double r) {
    return 1.0 * r * (r + 1) * (2.0 * r + 1) / 6.0
        - 1.0 * l * (l - 1) * (2.0 * l - 1) / 6.0;
}
inline void downdate(const int x) {
    int lc = T[x].lc, rc = T[x].rc;
    if (T[x].cn) {
        T[lc].firX = T[x].firX;
        T[rc].firX = T[x].firX;
        T[lc].cn = T[rc].cn = 1; T[x].cn = 0;
        T[lc].addX = T[rc].addX = T[lc].inc = 0;
        T[lc].firY = T[x].firY;
        T[rc].firY = T[x].firY;
        T[lc].cn = T[rc].cn = 1; T[x].cn = 0;
        T[lc].addY = T[rc].addY = T[rc].inc = 0;
    }
    T[lc].addX += T[x].addX; T[rc].addX += T[x].addX; T[x].addX = 0;
    T[lc].addY += T[x].addY; T[rc].addY += T[x].addY; T[x].addY = 0;
    T[lc].inc += T[x].inc; T[rc].inc += T[x].inc; T[x].inc = 0;
}
inline void upt(const int x) {
    int lc = T[x].lc, rc = T[x].rc, l = T[x].l, r = T[x].r;
    double x1, y1, xy1, x21, x2, y2, xy2, x22; int mid = l + r >> 1;
    x1 = T[lc].sumX; x2 = T[rc].sumX;
    y1 = T[lc].sumY; y2 = T[rc].sumY;
    x21 = T[lc].sumX2; x22 = T[rc].sumX2;
    xy1 = T[lc].sumXY; xy2 = T[rc].sumXY;
    x21 += 2.0 * x1 * T[lc].addX + T[lc].addX * T[lc].addX * (mid - l + 1);
    x22 += 2.0 * x2 * T[rc].addX + T[rc].addX * T[rc].addX * (r - mid);
    xy1 += x1 * T[lc].addY + T[lc].addX * y1 +
        T[lc].addX * T[lc].addY * (mid - l + 1);
    xy2 += x2 * T[rc].addY + T[rc].addX * y2 +
        T[rc].addX * T[rc].addY * (r - mid);
    x1 += T[lc].addX * (mid - l + 1); x2 += T[rc].addX * (r - mid);
    y1 += T[lc].addY * (mid - l + 1); y2 += T[rc].addY * (r - mid);
    if (T[lc].inc) {
        double lx = T[lc].firX + T[lc].addX, ly = T[lc].firY + T[lc].addY;
        x21 = 1.0 * (mid - l + 1) * sqr(lx) + lx * (l + mid) *
            (mid - l + 1) + sumx2(l, mid);
        xy1 = 1.0 * lx * ly * (mid - l + 1) + (lx + ly) * (l + mid)
            * (mid - l + 1) / 2.0 + sumx2(l, mid);
        x1 = lx * (mid - l + 1) + 1.0 * (l + mid) * (mid - l + 1) / 2.0;
        y1 = ly * (mid - l + 1) + 1.0 * (l + mid) * (mid - l + 1) / 2.0;
    }
    if (T[rc].inc) {
        double rx = T[rc].firX + T[rc].addX, ry = T[rc].firY + T[rc].addY;
        x22 = 1.0 * (r - mid) * sqr(rx) + rx * (mid + 1 + r) *
            (r - mid) + sumx2(mid + 1, r);
        xy2 = 1.0 * rx * ry * (r - mid) + (rx + ry) * (mid + 1 + r)
            * (r - mid) / 2.0 + sumx2(mid + 1, r);
        x2 = rx * (r - mid) + 1.0 * (mid + 1 + r) * (r - mid) / 2.0;
        y2 = ry * (r - mid) + 1.0 * (mid + 1 + r) * (r - mid) / 2.0;
    }
    T[x].sumX = x1 + x2; T[x].sumY = y1 + y2;
    T[x].sumXY = xy1 + xy2; T[x].sumX2 = x21 + x22;
}
inline void build(const int l, const int r, const int p) {
    T[p].l = l; T[p].r = r;
    if (l == r) {
        T[p].sumX = X[l]; T[p].sumY = Y[l];
        T[p].sumXY = T[p].sumX * T[p].sumY;
        T[p].sumX2 = T[p].sumX * T[p].sumX;
        return;
    }
    int p2 = (T[p].lc = ++QAQ), p3 = (T[p].rc = ++QAQ);
    int mid = l + r >> 1; build(l, mid, p2); build(mid + 1, r, p3);
    T[p].sumX = T[p2].sumX + T[p3].sumX;
    T[p].sumY = T[p2].sumY + T[p3].sumY;
    T[p].sumXY = T[p2].sumXY + T[p3].sumXY;
    T[p].sumX2 = T[p2].sumX2 + T[p3].sumX2;
}
inline void change(const int s, const int e,
const int Si, const int Ti, const int p) {
    int l = T[p].l, r = T[p].r;
    if (l == s && r == e) {
        T[p].addX += 1.0 * Si; T[p].addY += 1.0 * Ti;
        return;
    }
    int mid = l + r >> 1, p2 = T[p].lc, p3 = T[p].rc; downdate(p);
    if (e <= mid) change(s, e, Si, Ti, p2);
    else if (s >= mid + 1) change(s, e, Si, Ti, p3);
    else change(s, mid, Si, Ti, p2), change(mid + 1, e, Si, Ti, p3);
    upt(p);
}
inline void modify(const int s, const int e, const int Si,
const int Ti, const int p) {
    int l = T[p].l, r = T[p].r;
    if (T[p].l == s && T[p].r == e) {
        T[p].cn = 1; T[p].firX = 1.0 * Si; T[p].firY = 1.0 * Ti; T[p].inc = 1;
        T[p].addX = T[p].addY = 0;
        return;
    }
    int mid = l + r >> 1, p2 = T[p].lc, p3 = T[p].rc; downdate(p);
    if (e <= mid) modify(s, e, Si, Ti, p2);
    else if (s >= mid + 1) modify(s, e, Si, Ti, p3);
    else modify(s, mid, Si, Ti, p2), modify(mid + 1, e, Si, Ti, p3);
    upt(p);
}
inline pyz ask(const int s, const int e, const int p) {
    pyz res; int l = T[p].l, r = T[p].r;
    if (l == s && r == e) {
        double x, y;
        if (T[p].cn) {
            x = T[p].firX + T[p].addX; y = T[p].firY + T[p].addY;
            res.sumX = x * (r - l + 1) +
                1.0 * (l + r) * (r - l + 1)
                / 2.0;
            res.sumY = y * (r - l + 1) +
                1.0 * (l + r) * (r - l + 1)
                / 2.0;
            res.sumXY = x * y * (r - l + 1) +
                1.0 * (x + y) * (l + r) * (r - l + 1)
                / 2.0 + sumx2(l, r);
            res.sumX2 = sqr(x) * (r - l + 1) +
                x * (l + r) * (r - l + 1) + sumx2(l, r);
            return res;
        }
        res.sumX = T[p].sumX + T[p].addX * (r - l + 1);
        res.sumY = T[p].sumY + T[p].addY * (r - l + 1);
        res.sumXY = T[p].sumXY + T[p].sumX * T[p].addY
            + T[p].sumY * T[p].addX + T[p].addX * T[p].addY * (r - l + 1);
        res.sumX2 = T[p].sumX2 + 2.0 * T[p].addX * T[p].sumX
            + sqr(T[p].addX) * (r - l + 1);
        return res;
    }
    int mid = l + r >> 1, p2 = T[p].lc, p3 = T[p].rc; downdate(p);
    if (e <= mid) res = ask(s, e, p2);
    else if (s >= mid + 1) res = ask(s, e, p3);
    else res = ask(s, mid, p2) + ask(mid + 1, e, p3);
    return upt(p), res;
}
inline double trans(const int l, const int r, const pyz res) {
    double x = res.sumX, y = res.sumY, x2 = res.sumX2, xy = res.sumXY;
    return (xy - x * y / (r - l + 1)) / (x2 - x * x / (r - l + 1));
}
int main() {
    int i, op, x, y, u, v; n = read(); m = read();
    for (i = 1; i <= n; i++) X[i] = read();
    for (i = 1; i <= n; i++) Y[i] = read();
    build(1, n, 1);
    while (m--) {
        op = read(); if (op == 2) {
            x = read(); y = read(); u = read(); v = read();
            change(x, y, u, v, 1);
        }
        else if (op == 3) {
            x = read(); y = read(); u = read(); v = read();
            modify(x, y, u, v, 1);
        }
        else {
            x = read(); y = read();
            printf("%.10lf\n", trans(x, y, ask(x, y, 1)));
        }
    }
    return 0;
}