首页 > 代码库 > [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]相关分析
声明：以上内容来自用户投稿及互联网公开渠道收集整理发布，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任，若内容有误或涉及侵权可进行投诉： 投诉/举报 工作人员会在5个工作日内联系你，一经查实，本站将立刻删除涉嫌侵权内容。