bzoj 4821: [Sdoi2017]相关分析 线段树

题目:

Frank对天文学非常感兴趣，他经常用望远镜看星星，同时记录下它们的信息，比如亮度、颜色等等，进而估算出星星的距离，半径等等。Frank不仅喜欢观测，还喜欢分析观测到的数据。他经常分析两个参数之间（比如亮度和半径）是否存在某种关系。现在Frank要分析参数X与Y之间的关系。他有n组观测数据，第i组观测数据记录了x_i和y_i。他需要一下几种操作
1 L,R：用直线拟合第L组到底R组观测数据。
用\(xx\)表示这些观测数据中\(x\)的平均数，用\(yy\)表示这些观测数据中\(y\)的平均数，即
xx = Σx_i/(R-L+1)(L<=i<=R)
yy = Σy_i/(R-L+1)(L<=i<=R)
如果直线方程是y=ax+b，那么a应当这样计算：
a = (Σ(x_i-xx)(y_i-yy))/(Σ(x_i-xx)(x_i-xx)) (L<=i<=R)
你需要帮助Frank计算a。
2 L,R,S,T：
Frank发现测量数据第L组到底R组数据有误差，对每个i满足L <= i <= R，x_i需要加上S，y_i需要加上T。
3 L,R,S,T：
Frank发现第L组到第R组数据需要修改，对于每个i满足L <= i <= R，x_i需要修改为(S+i)，y_i需要修改为(T+i)。
1<=n,m<=10^{5,0<=|S|,|T|,|x_i|,|y_i|<=10}5

题解:

\[\begin{align} ans & = \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}^r(x_iy_i-x_i\bar{y}-y_i\bar{x}+\bar{x}\bar{y})}{\sum_{i=l}^r(x_i^2-2x_i\bar{x}+\bar{x}^2)} \\ & = \frac{\sum_{i=1}^rx_iy_i-\bar{y}\sum_{i=l}^rx_i - \bar{x}\sum_{i=l}^ry_i+(r-l+1)\bar{x}\bar{y}}{\sum_{i=l}^rx_i^2-2\bar{x}\sum_{i=l}^rx_i + (r-l+1)\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}} \end{align}\]

所以我们需要维护\(\sum x_i,\sum y_i,\sum x_iy_i,\sum x_i^2\)

很容易处理操作二:

• 对于\(\sum x_iy_i\)加上\(val\sum y_i\)
• 对于\(\sum x_i^2\),有\((x + val)^2 = x^2 + 2x*val + vval^2\)所以加上增量\(2*val\sum x,val^2\)即可

对于操作三,转化成一个赋值操作和一个操作二.
对于赋值操作..可以发现赋值后有\(\sum x_iy_i = \sum x_i^2\)
其实就是一个k次前缀和.k = 1 or 2
我们有:
\(\sum_{i=1}^ni = \frac{n(n+1)}{2}\)
\(\sum_{i=1}^ni^2 = \frac{n(n+1)(2n+1)}{6}\)
所以直接公式计算即可.

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
inline void read(int &x){
	x = 0;static char ch,f;f = 1;
	while(ch=getchar(),ch<'!');if(ch == '-') ch=getchar(),f = -1;
	while(x=(x<<1)+(x<<3)+ch-'0',ch=getchar(),ch>'!');x *= f;
}
#define rg register int 
#define rep(i,a,b) for(rg i=(a);i<=(b);++i)
#define per(i,a,b) for(rg i=(a);i>=(b);--i)
const int maxn = 100010;
struct Node{
	double sumx,sumy,misum,xysum;
	double lazy_x,lazy_y;bool tag;
	Node(){
		sumx = sumy = misum = xysum = lazy_x = lazy_y = 0;
		tag = 0;
	}
	friend Node operator + (const Node &a,const Node &b){
		Node c;
		c.sumx = a.sumx + b.sumx;
		c.sumy = a.sumy + b.sumy;
		c.misum = a.misum + b.misum;
		c.xysum = a.xysum + b.xysum;
		return c;
	}
}T[maxn<<2];
#define sum2(x) ( 1LL*(x)*((x)+1)*(2*(x)+1)/6 )
#define sum(x) (1LL*(x)*((x)+1)/2)
inline void add_tag(int rt,int l,int r){
	T[rt].misum = T[rt].xysum = sum2(r) - sum2(l-1);
	T[rt].sumx = T[rt].sumy = sum(r) - sum(l-1);
	T[rt].lazy_x = T[rt].lazy_y = 0;T[rt].tag = true;
}
inline void add_lazy_x(int rt,int l,int r,double val){
	T[rt].lazy_x += val;
	T[rt].misum += 2LL*val*T[rt].sumx + 1LL*val*val*(r-l+1);
	T[rt].xysum += 1LL*val*T[rt].sumy;
	T[rt].sumx += 1LL*val*(r-l+1); 
}
inline void add_lazy_y(int rt,int l,int r,double val){
	T[rt].lazy_y += val;
	T[rt].xysum += 1LL*val*T[rt].sumx;
	T[rt].sumy += 1LL*val*(r-l+1);
}
inline void pushdown(int rt,int l,int r){
	if(rt == 0 || l == r) return ;
	int mid = l+r >> 1;
	if(T[rt].tag){
		add_tag(rt<<1,l,mid);
		add_tag(rt<<1|1,mid+1,r);
		T[rt].tag = 0;
	}
	if(T[rt].lazy_x){
		add_lazy_x(rt<<1,l,mid,T[rt].lazy_x);
		add_lazy_x(rt<<1|1,mid+1,r,T[rt].lazy_x);
		T[rt].lazy_x = 0;
	}
	if(T[rt].lazy_y){
		add_lazy_y(rt<<1,l,mid,T[rt].lazy_y);
		add_lazy_y(rt<<1|1,mid+1,r,T[rt].lazy_y);
		T[rt].lazy_y = 0;
	}
	return ;
}
int L,R,vx,vy;
inline Node query(int rt,int l,int r){
	if(L <= l && r <= R) return T[rt];
	int mid = l+r >> 1;pushdown(rt,l,r);
	if(R <= mid) return query(rt<<1,l,mid);
	if(L >  mid) return query(rt<<1|1,mid+1,r);
	return query(rt<<1,l,mid) + query(rt<<1|1,mid+1,r);
}
inline void modify(int rt,int l,int r){
	if(L <= l && r <= R){
		add_lazy_x(rt,l,r,vx);
		add_lazy_y(rt,l,r,vy);
		return ;
	}
	int mid = l+r >> 1;pushdown(rt,l,r);
	if(L <= mid) modify(rt<<1,l,mid);
	if(R >  mid) modify(rt<<1|1,mid+1,r);
	T[rt] = T[rt<<1] + T[rt<<1|1];
}
inline void cover(int rt,int l,int r){
	if(L <= l && r <= R){
		add_tag(rt,l,r);
		return ;
	}
	int mid = l+r >> 1;pushdown(rt,l,r);
	if(L <= mid) cover(rt<<1,l,mid);
	if(R >  mid) cover(rt<<1|1,mid+1,r);
	T[rt] = T[rt<<1] + T[rt<<1|1];
}
int x[maxn],y[maxn];
inline void build(int rt,int l,int r){
	if(l == r){
		T[rt].sumx = x[l];
		T[rt].sumy = y[l];
		T[rt].misum = 1LL*x[l]*x[l];
		T[rt].xysum = 1LL*x[l]*y[l];
		return ;
	}
	int mid = l+r >> 1;
	build(rt<<1,l,mid);build(rt<<1|1,mid+1,r);
	T[rt] = T[rt<<1] + T[rt<<1|1];
}
int main(){
	int n,m;read(n);read(m);
	rep(i,1,n) read(x[i]);
	rep(i,1,n) read(y[i]);
	build(1,1,n);
	int op;double X,Y,len,a,b;
	Node tmp;
	while(m--){
		read(op);
		read(L);read(R);
		if(op == 1){
			tmp = query(1,1,n);
			X = tmp.sumx;
			Y = tmp.sumy;
			len = R - L + 1;
			a = tmp.xysum;
			b = tmp.misum;
			printf("%.10lf\n",(len*a - X*Y)/(len*b - X*X));
			continue;
		}
		read(vx);read(vy);
		if(op == 2){
			modify(1,1,n);
		}else{
			cover(1,1,n);
			modify(1,1,n);
		}
	}
	return 0;
}