首页 > 代码库 > [bzoj4821][Sdoi2017]相关分析

[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 doubleusing 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\n",(double)u/(double)d);        }        else        {            int s=read(),t=read();            M[i]=(Mark){op,s,t};            Modify(1,l,r,i);        }    }    return 0;}

[bzoj4821][Sdoi2017]相关分析