[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 =  , f = ; char ch = getchar();
    while(ch < '' || ch > ''){ if(ch == '-') f = -;  ch = getchar();}
    while(ch >= '' && ch <= ''){x = x *  + ch - '';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+];
struct Tree{int l,r,tag;ld s,t;data x;}T[MN*+];
int n,m,X[MN+],Y[MN+],tms=;
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>>;
    build(x<<,l,mid);build(x<<|,mid+,r);
    T[x].x=T[x<<].x+T[x<<|].x;
}
inline ld Sum(int x){return (ld)x*(x+)/;}
inline ld Sqr(int x){return (ld)x*(x+)*(*x+)/;}
void _Mark(int x,int op,int s,int t)
{
    if(op==){T[x].tag=;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=,T[x].s=s,T[x].t=t;
    if(op==)
    {
        T[x].x.sqx+=(ld)*s*T[x].x.x+(ld)(T[x].r-T[x].l+)*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+)*s*t;
        T[x].x.x+=(ld)(T[x].r-T[x].l+)*s;
        T[x].x.y+=(ld)(T[x].r-T[x].l+)*t;
    }
    else
    {
        T[x].x.x=Sum(T[x].r)-Sum(T[x].l-)+(ld)(T[x].r-T[x].l+)*s;
        T[x].x.y=Sum(T[x].r)-Sum(T[x].l-)+(ld)(T[x].r-T[x].l+)*t;
        T[x].x.sqx=(ld)(T[x].r-T[x].l+)*s*s+Sqr(T[x].r)-Sqr(T[x].l-)+(ld)*s*(Sum(T[x].r)-Sum(T[x].l-));
        T[x].x.xy=(ld)(T[x].r-T[x].l+)*s*t+(ld)(s+t)*(Sum(T[x].r)-Sum(T[x].l-))+Sqr(T[x].r)-Sqr(T[x].l-);
    }
}

void pushdown(int x)
{
    if(T[x].tag)
    {
        int l=x<<,r=l|;
        _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=;
    }
}

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>>;
    if(r<=mid) Modify(x<<,l,r,op);
    else if(l>mid) Modify(x<<|,l,r,op);
    else Modify(x<<,l,mid,op),Modify(x<<|,mid+,r,op);
    T[x].x=T[x<<].x+T[x<<|].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>>;
    if(r<=mid) Query(x<<,l,r);
    else if(l>mid) Query(x<<|,l,r);
    else Query(x<<,l,mid),Query(x<<|,mid+,r);
}

main()
{
    n=read();m=read();
    for(int i=;i<=n;++i)X[i]=read();
    for(int i=;i<=n;++i)Y[i]=read();
    build(,,n);
    for(int i=;i<=m;++i)
    {
        int op=read(),l=read(),r=read();
        if(op==)
        {
            res=(data){,,,};Query(,l,r);
            ld _x=(ld)res.x/(r-l+),_y=(ld)res.y/(r-l+);
            ld u=res.xy-_x*res.y-_y*res.x+_x*_y*(r-l+);
            ld d=res.sqx-*_x*res.x+_x*_x*(r-l+);
            printf("%.8lf\n",(double)u/(double)d);
        }
        else
        {
            int s=read(),t=read();
            M[i]=(Mark){op,s,t};
            Modify(,l,r,i);
        }
    }
    return ;
}