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即为正确。

式子可以化为a=((Σx_iy_i)-(R-L+1)*xx*yy)/((Σx_i²)-(R-L+1)*xx²) (L<=i<=R)

线段树维护一下区间内 x的和，y的和，xy的和，x²的和，然后打一下加法、覆盖标记就行了

#include<cstdio>
typedef long double i64;
const int N=100007;
int _l,_r,_a,_b;
i64 as[4];
i64 s2(i64 x){
    return x*(x+1)*(x*2+1)/6;
}
struct node{
    node*lc,*rc;
    int L,R,M,D;
    bool f;
    i64 x,y,x2,xy,xa,ya,xf,yf;
    void _add(i64 ax,i64 ay){
        xa+=ax;ya+=ay;
        x2+=2*ax*x+D*ax*ax;
        xy+=ay*x+ax*y+D*ax*ay;
        x+=D*ax;y+=D*ay;
    }
    void _fil(i64 fx,i64 fy){
        f=1;
        xf=fx;yf=fy;
        xa=ya=0;
        x2=s2(R+fx)-s2(L+fx-1);
        xy=D*fx*fy+(fx+fy)*(L+R)*D/2+s2(R)-s2(L-1);
        x=(L+R+fx*2)*D/2;
        y=(L+R+fy*2)*D/2;
    }
    void dn(){
        if(f){
            lc->_fil(xf,yf);
            rc->_fil(xf,yf);
            f=0;
        }
        if(xa||ya){
            lc->_add(xa,ya);
            rc->_add(xa,ya);
            xa=ya=0;
        }
    }
    void up(){
        x=lc->x+rc->x;
        y=lc->y+rc->y;
        x2=lc->x2+rc->x2;
        xy=lc->xy+rc->xy;
    }
    void add(){
        if(_l<=L&&R<=_r){
            _add(_a,_b);
            return;
        }
        dn();
        if(_l<=M)lc->add();
        if(_r>M)rc->add();
        up();
    }
    void fil(){
        if(_l<=L&&R<=_r){
            _fil(_a,_b);
            return;
        }
        dn();
        if(_l<=M)lc->fil();
        if(_r>M)rc->fil();
        up();
    }
    void get(){
        if(_l<=L&&R<=_r){
            as[0]+=x;
            as[1]+=y;
            as[2]+=x2;
            as[3]+=xy;
            return;
        }
        dn();
        if(_l<=M)lc->get();
        if(_r>M)rc->get();
    }
}ns[N*2],*np=ns,*rt;
int xs[N],ys[N];
node*build(int L,int R){
    node*w=np++;
    w->L=L;w->R=R;w->D=R-L+1;
    if(L<R){
        int M=w->M=L+R>>1;
        w->lc=build(L,M);
        w->rc=build(M+1,R);
        w->up();
    }else w->_add(xs[L],ys[L]);
    return w;
}
char buf[N*100],*ptr=buf-1;
int _(){
    int x=0,f=1,c=*++ptr;
    while(c<48)c=='-'&&(f=-1),c=*++ptr;
    while(c>47)x=x*10+c-48,c=*++ptr;
    return x*f;
}
int n,m;
int main(){
    fread(buf,1,sizeof(buf),stdin)[buf]=0;
    n=_();m=_();
    for(int i=1;i<=n;++i)xs[i]=_();
    for(int i=1;i<=n;++i)ys[i]=_();
    rt=build(1,n);
    for(int i=0,o;i<m;++i){
        o=_();_l=_();_r=_();
        if(o==1){
            as[0]=as[1]=as[2]=as[3]=0;
            rt->get();
            i64 d=_r-_l+1;
            i64 r=(d*as[3]-as[0]*as[1])/(d*as[2]-as[0]*as[0]);
            printf("%.10Lf\n",r);
        }else if(o==2){
            _a=_();_b=_();
            rt->add();
        }else{
            _a=_();_b=_();
            rt->fil();
        }
    }
    return 0;
}