世界中心量子力学及红石自动化研究中心

4821: [Sdoi2017]相关分析

Time Limit: 10 Sec Memory Limit: 128 MB Special Judge

Description

Frank对天文学非常感兴趣，他经常用望远镜看星星，同时记录下它们的信息，比如亮度、颜色等等，进而估算出

星星的距离，半径等等。Frank不仅喜欢观测，还喜欢分析观测到的数据。他经常分析两个参数之间（比如亮度和

半径）是否存在某种关系。现在Frank要分析参数$X$与$Y$之间的关系。他有$n$组观测数据，第$i$组观测数据记录了$x_i$和

$y_i$。他需要一下几种操作：

$1$ $L$,$R$：用直线拟合第$L$组到底$R$组观测数据。用$\overline x$表示这些观测数据中$x$的平均数，用$\overline y$

表示这些观测数据中$y$的平均数，即：

$$\overline x=\frac1{R-L+1}\sum_{i=L}^Rx_i$$

$$\overline y=\frac1{R-L+1}\sum_{i=L}^Ry_i$$

如果直线方程是$y=ax+b$，那么$a$应当这样计算：

$$\left\{\begin{array}{l}a=\frac{{\displaystyle\sum_{i=L}^R}\left(x_i-\overline x\right)\left(y_i-\overline y\right)}{{\displaystyle\sum_{i=L}^R}\left(x_i-\overline x\right)^2}\\b=\overline y-a\overline x\end{array}\right.$$

你需要帮助Frank计算$a$。

$2$ $L$,$R$,$S$,$T$：

Frank发现测量数据第$L$组到第$R$组数据有误差，对每个$i$满足$L \leqslant i\leqslant R$，$x_i$需要加上$S$，$y_i$需要加上$T$。

$3$ $L$,$R$,$S$,$T$：

Frank发现第$L$组到第$R$组数据需要修改，对于每个$i$满足$L \leqslant i\leqslant R$，$x_i$需要修改为$S+i$，$y_i$需要修改为$T+i$。

Input

第一行两个数$n$,$m$，表示观测数据组数和操作次数。

接下来一行$n$个数，第$i$个数是$x_i$。

接下来一行$n$个数，第$i$个数是$y_i$。

接下来$m$行，表示操作，格式见题目描述。

保证$1$操作不会出现分母为$0$的情况。

Output

对于每个$1$操作，输出一行，表示直线斜率$a$。

选手输出与标准输出的绝对误差不超过$10^{-5}$即为正确。

Sample Input

3 5

1 2 3

1 1 3

2 2 3 -3 2

1 1 2

3 1 2 2 1

1 1 3

Sample Output

1.0000000000

-1.5000000000

-0.6153846154

Solution

思考过如何快速求区间标准差的读者应该可以很容易地想到此题的思路：用线段树维护一些值，根据这些值计算答案。

观察a的表达式，容易得到：

$$a=\frac{\displaystyle\sum_{i=L}^R\left(x_iy_i-\overline xy_i-\overline yx_i-\overline x\cdot\overline y\right)}{\displaystyle\sum_{i=L}^R\left(x_i^2-2x_i\overline x+\overline x^2\right)}$$

平均值可以通过维护区间和得到，为了计算答案，我们还需要知道$\sum_{i=L}^Rx_i^2$和$\sum_{i=L}^Rx_iy_i$。于是我们可以用线段树维护4个值：$sumx$、$sumy$、$sumxx$、$sumxy$分别表示$x$的区间和、$y$的区间和、$x^2$的区间和、$xy$的区间和。

代码如下：

struct query_t

{

query_t operator+(const query_t b)const

{

return (query_t){

sumx+b.sumx,sumy+b.sumy,

sumxx+b.sumxx,sumxy+b.sumxy

};

}

long double sumx,sumy,sumxx,sumxy;

};

struct

{

int l,r;

query_t d;

long double flag_addx,flag_addy;

long double flag_updatax,flag_updatay;

}tree[400001];

int l,r;

read(l);read(r);

query_t result=query(l,r);

long double n=r-l+1;

long double ave_x=result.sumx/n;

long double ave_y=result.sumy/n;

long double a=

(result.sumxy-result.sumy*ave_x-result.sumx*ave_y+n*ave_x*ave_y)/

(result.sumxx-2*ave_x*result.sumx+n*ave_x*ave_x);

cout<<a<<"\n";

我们来考虑如何实现区间修改。

首先考虑操作2：设$x$的增加量为$k$，$y$的增加量为$g$，修改后$sumxx'=\sum(x_i+k)^2=\sum x_i^2+2k\sum x_i+\sum k^2=sumxx+2k\cdot sumx+nk^2$，$n$为区间长度，类似地，$sumxy'=\sum(x_i+{k)(y_i+g)}=\sum(x_i+y_i)+g\sum x_i+k\sum y_i+\sum k\cdot g=sumxy+g\cdot sumx+k\cdot sumy+n\cdot k\cdot g$。

操作3的实现比较简单，$sumxx'=\sum_{i=l}^r(k+i)^2$，$sumxx=\sum_{i=l}^r\left(k+i\right)\left(g+i\right)=n\cdot k\cdot g+\left(k+g\right)\frac{n\left(l+r\right)}2+\sum_{i=l}^ri^2$。pushdown函数如下：

inline void pushdown(int now)

{

if(tree[now].flag_updatax!=INF)

{

for(int child=now<<1;child<=(now<<1|1);child++)

{

long double l=tree[child].l,r=tree[child].r;

long double k=tree[now].flag_updatax,g=tree[now].flag_updatay;

long double n=r-l+1;

tree[child].flag_updatax=tree[now].flag_updatax;

tree[child].flag_updatay=tree[now].flag_updatay;

////////////////////////////////////

tree[child].d.sumxx=sumsqr(k+l,k+r);

tree[child].d.sumxy=n*k*g+(k+g)*(n*(l+r)/2)+sumsqr(l,r);

tree[child].d.sumx=(l+k+r+k)*n/2;

tree[child].d.sumy=(l+g+r+g)*n/2;

}

else if(tree[now].flag_addx!=0||tree[now].flag_addy!=0)

{

for(int child=now<<1;child<=(now<<1|1);child++)

{

long double l=tree[child].l,r=tree[child].r;

long double k=tree[now].flag_addx,g=tree[now].flag_addy;

long double n=r-l+1;

if(tree[child].flag_updatax!=INF)tree[child].flag_updatax+=k;

else tree[child].flag_addx+=k;

if(tree[child].flag_updatay!=INF)tree[child].flag_updatay+=g;

else tree[child].flag_addy+=g;

////////////////////////

tree[child].d.sumxx+=2ll*k*tree[child].d.sumx+n*k*k;

tree[child].d.sumxy+=k*tree[child].d.sumy+g*tree[child].d.sumx+n*k*g;

tree[child].d.sumx+=n*k;

tree[child].d.sumy+=n*g;

}

tree[now].flag_addx=0;

tree[now].flag_addy=0;

tree[now].flag_updatax=INF;

tree[now].flag_updatay=INF;

}

这里我们看到，flag_add被优先下传至flag_updata，这是因为操作3造成的影响可“抹去”之前的操作2的影响。

其他部分都是线段树的基本操作，不再详细讲解。

Code

#include <iostream>

#include <stdio.h>

#include <string.h>

#include <algorithm>

#include <vector>

#include <iomanip>

using namespace std;

template <typename tp>

void read(tp &x)

{

x=0;char ch=getchar();bool f=0;

while(ch>'9'||ch<'0'){ if(ch=='-')f=1;ch=getchar(); }

while(ch>='0'&&ch<='9'){ x=x*10+ch-48;ch=getchar(); }

if(f)x=-x;

}

long double sumsqr(long double l,long double r)

{

if(l==0)return r*(r+1)*(2*r+1)/6;

if(l>0)

{

long double sum2=r*(r+1)/2;

sum2*=(2*r+1);

sum2/=3;

long double sum1=(l-1)*l/2;

sum1*=(2*l-1);

sum1/=3;

return sum2-sum1;

}

if(r<=0)return sumsqr(-r,-l);

return sumsqr(0,-l)+sumsqr(0,r);

}

//我并不知道这个值对于long double究竟等于多少，所以不建议读者这样写

const long double INF=0x3f3f3f3f3f3f3f3f;

struct query_t

{

query_t operator+(const query_t b)const

{

return (query_t){

sumx+b.sumx,sumy+b.sumy,

sumxx+b.sumxx,sumxy+b.sumxy

};

}

long double sumx,sumy,sumxx,sumxy;

};

struct

{

int l,r;

query_t d;

long double flag_addx,flag_addy;

long double flag_updatax,flag_updatay;

}tree[400001];

long double x_tmp[100001],y_tmp[100001];

void build(int l,int r,int now=1)

{

tree[now].l=l,tree[now].r=r;

tree[now].flag_addx=0,tree[now].flag_updatax=INF;

tree[now].flag_addy=0,tree[now].flag_updatay=INF;

if(l==r)

{

tree[now].d.sumx=x_tmp[l];

tree[now].d.sumy=y_tmp[l];

tree[now].d.sumxx=tree[now].d.sumx*tree[now].d.sumx;

tree[now].d.sumxy=tree[now].d.sumx*tree[now].d.sumy;

return;

}

int mid=l+r>>1;

build(l,mid,now<<1);

build(mid+1,r,now<<1|1);

tree[now].d=tree[now<<1].d+tree[now<<1|1].d;

}

inline void pushdown(int now)

{

if(tree[now].flag_updatax!=INF)

{

for(int child=now<<1;child<=(now<<1|1);child++)

{

long double l=tree[child].l,r=tree[child].r;

long double k=tree[now].flag_updatax,g=tree[now].flag_updatay;

long double n=r-l+1;

tree[child].flag_updatax=tree[now].flag_updatax;

tree[child].flag_updatay=tree[now].flag_updatay;

////////////////////////////////////

tree[child].d.sumxx=sumsqr(k+l,k+r);

tree[child].d.sumxy=n*k*g+(k+g)*(n*(l+r)/2)+sumsqr(l,r);

tree[child].d.sumx=(l+k+r+k)*n/2;

tree[child].d.sumy=(l+g+r+g)*n/2;

}

else if(tree[now].flag_addx!=0||tree[now].flag_addy!=0)

{

for(int child=now<<1;child<=(now<<1|1);child++)

{

long double l=tree[child].l,r=tree[child].r;

long double k=tree[now].flag_addx,g=tree[now].flag_addy;

long double n=r-l+1;

if(tree[child].flag_updatax!=INF)tree[child].flag_updatax+=k;

else tree[child].flag_addx+=k;

if(tree[child].flag_updatay!=INF)tree[child].flag_updatay+=g;

else tree[child].flag_addy+=g;

////////////////////////

tree[child].d.sumxx+=2ll*k*tree[child].d.sumx+n*k*k;

tree[child].d.sumxy+=k*tree[child].d.sumy+g*tree[child].d.sumx+n*k*g;

tree[child].d.sumx+=n*k;

tree[child].d.sumy+=n*g;

}

tree[now].flag_addx=0;

tree[now].flag_addy=0;

tree[now].flag_updatax=INF;

tree[now].flag_updatay=INF;

}

query_t query(int l,int r,int now=1)

{

if(tree[now].l==l&&tree[now].r==r)

{

return tree[now].d;

}

int mid=tree[now].l+tree[now].r>>1;

pushdown(now);

if(r<=mid)return query(l,r,now<<1);

else if(mid<l)return query(l,r,now<<1|1);

else return query(l,mid,now<<1)+query(mid+1,r,now<<1|1);

}

void add(int l,int r,long double dx,long double dy,int now=1)

{

if(tree[now].l==l&&tree[now].r==r)

{

if(tree[now].flag_updatax!=INF)tree[now].flag_updatax+=dx;

else tree[now].flag_addx+=dx;

if(tree[now].flag_updatay!=INF)tree[now].flag_updatay+=dy;

else tree[now].flag_addy+=dy;

///////////////

long double n=r-l+1;

tree[now].d.sumxx+=2ll*dx*tree[now].d.sumx+n*dx*dx;

tree[now].d.sumxy+=dx*tree[now].d.sumy+dy*tree[now].d.sumx+n*dx*dy;

tree[now].d.sumx+=n*dx;

tree[now].d.sumy+=n*dy;

return;

}

int mid=tree[now].l+tree[now].r>>1;

pushdown(now);

if(r<=mid)add(l,r,dx,dy,now<<1);

else if(mid<l)add(l,r,dx,dy,now<<1|1);

else

{

add(l,mid,dx,dy,now<<1);

add(mid+1,r,dx,dy,now<<1|1);

}

tree[now].d=tree[now<<1].d+tree[now<<1|1].d;

}

void updata(int l,int r,long double dx,long double dy,int now=1)

{

if(tree[now].l==l&&tree[now].r==r)

{

tree[now].flag_updatax=dx;

tree[now].flag_updatay=dy;

///////////////

long double n=r-l+1;

tree[now].d.sumxx=sumsqr(dx+l,dx+r);

tree[now].d.sumxy=n*dx*dy+(dx+dy)*(n*(l+r)/2)+sumsqr(l,r);

tree[now].d.sumx=(l+dx+r+dx)*n/2;

tree[now].d.sumy=(l+dy+r+dy)*n/2;

return;

}

int mid=tree[now].l+tree[now].r>>1;

pushdown(now);

if(r<=mid)updata(l,r,dx,dy,now<<1);

else if(mid<l)updata(l,r,dx,dy,now<<1|1);

else

{

updata(l,mid,dx,dy,now<<1);

updata(mid+1,r,dx,dy,now<<1|1);

}

tree[now].d=tree[now<<1].d+tree[now<<1|1].d;

}

int main()

{

int n,Q;

read(n);read(Q);

for(int i=1;i<=n;i++)read(x_tmp[i]);

for(int i=1;i<=n;i++)read(y_tmp[i]);

build(1,n);

cout<<setprecision(15)<<fixed;

while(Q--)

{

int type;

read(type);

if(type==1)

{

int l,r;

read(l);read(r);

query_t result=query(l,r);

long double n=r-l+1;

long double ave_x=result.sumx/n;

long double ave_y=result.sumy/n;

long double a=

(result.sumxy-result.sumy*ave_x-result.sumx*ave_y+n*ave_x*ave_y)/

(result.sumxx-2*ave_x*result.sumx+n*ave_x*ave_x);

cout<<a<<"\n";

}

else if(type==2)

{

int l,r;

long long s,t;

read(l);read(r);read(s);read(t);

add(l,r,s,t);

}

else

{

int l,r;

long long s,t;

read(l);read(r);read(s);read(t);

updata(l,r,s,t);

}

return 0;

}

BZOJ 4821 相关分析 解题报告

BZOJ 4821 相关分析解题报告