[bzoj4821][Sdoi2017]相关分析

来自FallDream的博客，未经允许，请勿转载，谢谢。

---

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

 

展开式子，答案是

$frac{sum(xiyi -xyi - yxi + xy)}{sum (xi^{2} -2xxi + x^{2})}$

所以维护区间的x,y,xy的和就行了。

x^2维不维护比较随意。

计算记得用longdouble  不然会爆掉

#include<iostream>
#include<cstdio>
#define MN 100000
#define ld long double
using namespace std;
inline int read()
{
    int x = 0 , f = 1; char ch = getchar();
    while(ch < '0' || ch > '9'){ if(ch == '-') f = -1;  ch = getchar();}
    while(ch >= '0' && ch <= '9'){x = x * 10 + ch - '0';ch = getchar();}
    return x * f;
}

struct data{ld x,y,sqx,xy;
    friend data operator + (data a,data b)
    {
        return (data){a.x+b.x,a.y+b.y,a.sqx+b.sqx,a.xy+b.xy};    
    }
}res;
struct Mark{int op,s,t;}M[MN+5];
struct Tree{int l,r,tag;ld s,t;data x;}T[MN*4+5]; 
int n,m,X[MN+5],Y[MN+5],tms=0;
void build(int x,int l,int r)
{
    if((T[x].l=l)==(T[x].r=r))
    {
        T[x].x=(data){X[l],Y[l],(ld)X[l]*X[l],(ld)X[l]*Y[l]};
        return;
    }
    int mid=l+r>>1;
    build(x<<1,l,mid);build(x<<1|1,mid+1,r);
    T[x].x=T[x<<1].x+T[x<<1|1].x;
}
inline ld Sum(int x){return (ld)x*(x+1)/2;}
inline ld Sqr(int x){return (ld)x*(x+1)*(2*x+1)/6;}
void _Mark(int x,int op,int s,int t)
{
    if(op==3){T[x].tag=3;T[x].s=s;T[x].t=t;}
    else if(T[x].tag) T[x].s+=s,T[x].t+=t;
    else T[x].tag=2,T[x].s=s,T[x].t=t;
    if(op==2)
    {
        T[x].x.sqx+=(ld)2*s*T[x].x.x+(ld)(T[x].r-T[x].l+1)*s*s;
        T[x].x.xy+=(ld)s*T[x].x.y+(ld)t*T[x].x.x+(ld)(T[x].r-T[x].l+1)*s*t;
        T[x].x.x+=(ld)(T[x].r-T[x].l+1)*s;
        T[x].x.y+=(ld)(T[x].r-T[x].l+1)*t;
    }
    else
    {
        T[x].x.x=Sum(T[x].r)-Sum(T[x].l-1)+(ld)(T[x].r-T[x].l+1)*s;
        T[x].x.y=Sum(T[x].r)-Sum(T[x].l-1)+(ld)(T[x].r-T[x].l+1)*t;
        T[x].x.sqx=(ld)(T[x].r-T[x].l+1)*s*s+Sqr(T[x].r)-Sqr(T[x].l-1)+(ld)2*s*(Sum(T[x].r)-Sum(T[x].l-1));
        T[x].x.xy=(ld)(T[x].r-T[x].l+1)*s*t+(ld)(s+t)*(Sum(T[x].r)-Sum(T[x].l-1))+Sqr(T[x].r)-Sqr(T[x].l-1); 
    }
}

void pushdown(int x)
{
    if(T[x].tag)
    {
        int l=x<<1,r=l|1; 
        _Mark(l,T[x].tag,T[x].s,T[x].t);
        _Mark(r,T[x].tag,T[x].s,T[x].t);    
        T[x].tag=0;
    }
}

void Modify(int x,int l,int r,int op)
{
    if(T[x].l==l&&T[x].r==r)
    {
        _Mark(x,M[op].op,M[op].s,M[op].t);
        return;    
    }
    pushdown(x);
    int mid=T[x].l+T[x].r>>1;
    if(r<=mid) Modify(x<<1,l,r,op);
    else if(l>mid) Modify(x<<1|1,l,r,op);
    else Modify(x<<1,l,mid,op),Modify(x<<1|1,mid+1,r,op);
    T[x].x=T[x<<1].x+T[x<<1|1].x;
}    

void Query(int x,int l,int r)
{
    if(T[x].l==l&&T[x].r==r){res=res+T[x].x;return;}    
    pushdown(x);
    int mid=T[x].l+T[x].r>>1;
    if(r<=mid) Query(x<<1,l,r);
    else if(l>mid) Query(x<<1|1,l,r);
    else Query(x<<1,l,mid),Query(x<<1|1,mid+1,r);
}

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();
    build(1,1,n);
    for(int i=1;i<=m;++i)
    {
        int op=read(),l=read(),r=read();        
        if(op==1)
        {
            res=(data){0,0,0,0};Query(1,l,r); 
            ld _x=(ld)res.x/(r-l+1),_y=(ld)res.y/(r-l+1);
            ld u=res.xy-_x*res.y-_y*res.x+_x*_y*(r-l+1);
            ld d=res.sqx-2*_x*res.x+_x*_x*(r-l+1);
            printf("%.8lf
",(double)u/(double)d);
        }
        else
        {
            int s=read(),t=read();
            M[i]=(Mark){op,s,t};
            Modify(1,l,r,i);
        }
    }
    return 0;
}