[BZOJ4821][SDOI2017]相关分析

[BZOJ4821][SDOI2017]相关分析

相关分析

Description

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)。

Input

第一行两个数n,m，表示观测数据组数和操作次数。
接下来一行n个数，第i个数是x_i。
接下来一行n个数，第i个数是y_i。
接下来m行，表示操作，格式见题目描述。
1<=n,m<=10^5,0<=|S|,|T|,|x_i|,|y_i|<=10^5
保证1操作不会出现分母为0的情况。

Output

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

Sample Input

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

Sample Output

1.0000000000
-1.5000000000
-0.6153846154

HINT

Source

鸣谢infinityedge上传

---

写完去洛谷自信一交结果80pts……
然后检查半天无果，看讨论才知道是爆longlong了……

(╯‵□′)╯︵┻━┻

---

思路:
纯粹的线段树。

首先拆式子。具体方式为，把平均数全部换成总和除以个数，然后整理一下。
可以发现会拆出来这样一个东西:

a=∑xi∗yi−∑x∗∑yr−l+1∑x2i−(∑x)2r−l+1 $$a=\frac{\sum x_i*y_i-\frac{\sum x*\sum y}{r-l+1}}{\sum x_i^2-\frac{(\sum x)^2}{r-l+1}}$$

那么维护 8 $8$个值即可:
∑x,∑y,∑x2,∑x∗y$\sum x,\sum y,\sum x^2,\sum x*y$，以及两种操作对应的 S $S$和T$T$的tag。

然后就是推 push_down() $push\_down()$函数了。
显然 ∑x $\sum x$和 ∑y $\sum y$是非常简单的。

考虑2操作，此时有:

∑(x+s)2=∑x2+2sx+s2 $\sum (x+s)^2=\sum x^2+2sx+s^2$
于是可以利用之前的 ∑x $\sum x$和 ∑x2 $\sum x^2$，加上 s2 $s^2$求出。

∑(x+s)∗(y+t)=∑xy+∑sy+∑tx+∑st $\sum (x+s)*(y+t)=\sum xy + \sum sy + \sum tx + \sum st$
对于二、三两项可以利用之前的 ∑x $\sum x$和 ∑y $\sum y$求出。

然后是魔幻的3操作:

∑i=lr(s+i)2=∑i=lrs2+2si+i2=∑s2+2s∗∑i+∑i2 $\sum\limits_{i=l}^r (s+i)^2=\sum\limits_{i=l}^r s^2+2si+i^2 =\sum s^2+2s*\sum i +\sum i^2$
于是套用前缀和公式和平方和公式即可。

∑i=lr(s+i)(t+i)=∑i=lrst+(s+t)i+i2=∑st+∑(s+t)i+∑i2 $\sum\limits_{i=l}^r (s+i)(t+i)=\sum\limits_{i=l}^r st+(s+t)i+i^2=\sum st + \sum (s+t)i +\sum i^2$
于是同样套用前缀和公式和平方和公式即可。

最后，把线段树的 longlong $long long$改成 double $double$。
然后就可以过了~

#include<bits/stdc++.h>
using namespace std;typedef long long ll;
typedef double db;const int N=1e5+9;
const int hath=19260817;
int n,m;
db x[N],y[N];inline db read()
{db x=0,f=1;char ch=getchar();while(ch<'0' || '9'<ch){if(ch=='-')f=-1;ch=getchar();}while('0'<=ch && ch<='9')x=x*10+(ch^48),ch=getchar();return x*f;
}inline db sum(db l,db r){return (l+r)*(r-l+1)/2;}
inline db sqr(db n){return n*(n+1)*(2*n+1)/6;}namespace segt
{db tsumx[N<<2],tsumy[N<<2];db tsqr[N<<2],tyx[N<<2];db tagax[N<<2],tagay[N<<2];db tagsx[N<<2],tagsy[N<<2];inline void marks(int x,db l,db r,db s,db t){tagsx[x]=s;tagsy[x]=t;tagax[x]=tagay[x]=hath;tsumx[x]=(s+s+r-l)*(r-l+1)/2;tsumy[x]=(t+t+r-l)*(r-l+1)/2;tsqr[x]=sqr(s+r-l)-sqr(s-1);tyx[x]=(r-l+1)*s*t+sum(1,r-l)*(s+t)+sqr(r-l);}inline void marka(int x,db l,db r,db s,db t){if(tagax[x]==hath)tagax[x]=s,tagay[x]=t;else tagax[x]+=s,tagay[x]+=t;tsqr[x]+=tsumx[x]*2*s+s*s*(r-l+1);tyx[x]+=s*tsumy[x]+t*tsumx[x]+s*t*(r-l+1);tsumx[x]+=s*(r-l+1);tsumy[x]+=t*(r-l+1);}inline void push(int x,int l,int r){int mid=l+r>>1;if(tagsx[x]!=hath){marks(x<<1,l,mid,tagsx[x],tagsy[x]);marks(x<<1|1,mid+1,r,tagsx[x]+(mid-l+1),tagsy[x]+(mid-l+1));tagsx[x]=tagsy[x]=hath;}if(tagax[x]!=hath){marka(x<<1,l,mid,tagax[x],tagay[x]);marka(x<<1|1,mid+1,r,tagax[x],tagay[x]);tagax[x]=tagay[x]=hath;}}inline void update(int x){tsumx[x]=tsumx[x<<1]+tsumx[x<<1|1];tsumy[x]=tsumy[x<<1]+tsumy[x<<1|1];tsqr[x]=tsqr[x<<1]+tsqr[x<<1|1];tyx[x]=tyx[x<<1]+tyx[x<<1|1];}inline void build(int rt,int l,int r){tagax[rt]=tagay[rt]=hath;tagsx[rt]=tagsy[rt]=hath;if(l==r){tsumx[rt]=x[l];tyx[rt]=x[l]*y[l];tsumy[rt]=y[l];tsqr[rt]=x[l]*x[l];return;}int mid=l+r>>1;build(rt<<1,l,mid);build(rt<<1|1,mid+1,r);update(rt);}inline void modify(int x,int l,int r,int dl,int dr,int s,int t){if(l!=r)push(x,l,r);if(dl==l && r==dr){marks(x,l,r,s+l,t+l);return;}int mid=l+r>>1;if(dr<=mid)modify(x<<1,l,mid,dl,dr,s,t);else if(mid<dl)modify(x<<1|1,mid+1,r,dl,dr,s,t);else{modify(x<<1,l,mid,dl,mid,s,t);modify(x<<1|1,mid+1,r,mid+1,dr,s,t);}update(x);}inline void add(int x,int l,int r,int dl,int dr,int s,int t){if(l!=r)push(x,l,r);int mid=l+r>>1;if(dl==l && r==dr){marka(x,l,r,s,t);return;}if(dr<=mid)add(x<<1,l,mid,dl,dr,s,t);else if(mid<dl)add(x<<1|1,mid+1,r,dl,dr,s,t);else add(x<<1,l,mid,dl,mid,s,t),add(x<<1|1,mid+1,r,mid+1,dr,s,t);update(x);}inline db query(int x,int l,int r,int dl,int dr,db *t){if(l!=r)push(x,l,r);int mid=l+r>>1;if(dl==l && r==dr)return t[x];if(dr<=mid)return query(x<<1,l,mid,dl,dr,t);if(mid<dl)return query(x<<1|1,mid+1,r,dl,dr,t);return query(x<<1,l,mid,dl,mid,t)+query(x<<1|1,mid+1,r,mid+1,dr,t);}inline db check(int x,int l,int r,int p,db *t){if(l!=r)push(x,l,r);if(l==r)return t[x];int mid=l+r>>1;if(p<=mid)return check(x<<1,l,mid,p,t);else return check(x<<1|1,mid+1,r,p,t);}
}int main()
{n=read();m=read();for(int i=1;i<=n;i++)x[i]=read();for(int i=1;i<=n;i++)y[i]=read();segt::build(1,1,n);while(m--){int ty=read();ll l=read(),r=read();if(ty==1){db sumx=segt::query(1,1,n,l,r,segt::tsumx);db sumy=segt::query(1,1,n,l,r,segt::tsumy);db up=0,down=0;up+=segt::query(1,1,n,l,r,segt::tyx);up-=(db)sumx*sumy/((db)r-l+1.0);down+=segt::query(1,1,n,l,r,segt::tsqr);down-=(db)sumx*sumx/((db)r-l+1.0);printf("%.8f\n",up/down);}else{ll s=read(),t=read();if(ty==2)segt::add(1,1,n,l,r,s,t);elsesegt::modify(1,1,n,l,r,s,t);}}return 0;
}