题解方差

Posted 2020-11-25 ctjcalc

原文链接：https://cnblogs.com/ctjcalc/p/post6.html

这道题其实除了一些细节，就是个线段树板子题。虽然我WA了不知道多少次，只是因为不等号方向写反。

Copyright ? 2019 ctjcalc，转载请注明URL，并给出原文链接，谢谢。

算法分析

基本思路

以下是本蒟蒻的推导过程。

首先对方差公式进行变形：
[ egin{aligned} s^2&=frac{sum_{i=1}^{n}(x_i-ar{x})^2}{n} \\ &=frac{sum_{i=1}^{n}(x_i^2+ar{x}^2-2x_iar{x})}{n} \\ &=frac{sum_{i=1}^{n}x_i^2+sum_{i=1}^{n}ar{x}^2-2ar{x}sum_{i=1}^{n}x_i}{n} end{aligned} ]
这是对于整个序列来说的方差，很显然，在对一个区间上求方差也是可以的。

可以发现，我们只要维护区间和，区间平方和即可，果断上线段树。

Copyright ? 2019 ctjcalc，转载请注明URL，并给出原文链接，谢谢。

但是对于区间加的操作，平方和又要怎么维护呢？请看下面：
[ egin{aligned} sum_{i=l}^{r}(x_i+k)^2&=sum_{i=l}^{r}(x_i^2+k^2+2x_ik) &=sum_{i=l}^{r}x_i^2+sum_{i=l}^{r}k^2+2ksum_{i=l}^{r}x_i &=sum_{i=l}^{r}x_i^2+2ksum_{i=l}^{r}x_i+k^2(r-l+1) end{aligned} ]

所以就可以写出只要的pushdown函数：

inline void update(int x,int len,double v){
    tree[x].add+=v;
    tree[x].sqr+=len*square(v)+2*v*tree[x].sum; // Warning 1
    tree[x].sum+=len*v; 
}

inline void pushdown(int x,int l,int r) {
    if(abs(tree[x].add)>eps){ // Warning 2
        int mid=(l+r)>>1;
        update(ls(x),mid-l+1,tree[x].add);
        update(rs(x),r-mid,tree[x].add);
        tree[x].add=0.0;
    }
}

浮点数比较

Copyright ? 2019 ctjcalc，转载请注明URL，并给出原文链接，谢谢。

浮点数是个很坑的东西，由于计算机架构的特殊原因，浮点数很不精确，可以自己运行一下下面的代码：

#include <bits/stdc++.h>
using namespace std;

int main(){
    double x=3.14;
    printf("%.20lf",x);
    return 0;
}

我的输出为(3.14000000000000012434)，直接比较两个浮点数是很容易出锅的。

我们可以引入一个极小值(epsilon)，比如(10^{-8})，本题是可以接受的，然后判断两数差值的绝对值是否小于(epsilon)，如果小于，则两数相等，否则不等。

const double eps=1e-8;
bool equal(double x,double y){
    return abs(x-y)<eps;
}

Copyright ? 2019 ctjcalc，转载请注明URL，并给出原文链接，谢谢。

代码实现

Talk is cheap, show you the code.

#include <bits/stdc++.h>
#define FILE_IO true
#define TIME_LIMIT 1000
using namespace std;

struct node
{
    double sum,sqr,add;
};

const int maxn=100000+10;
const double eps=1e-8;
node tree[maxn<<2];
int n,m;
double a[maxn];

inline int ls(int x) { return x<<1; }

inline int rs(int x) { return x<<1|1; }

inline void pushup(int x){
    tree[x].sum=tree[ls(x)].sum+tree[rs(x)].sum;
    tree[x].sqr=tree[ls(x)].sqr+tree[rs(x)].sqr;
}

inline double square(double x) { return x*x; }

inline void update(int x,int len,double v){
    tree[x].add+=v;
    tree[x].sqr+=len*square(v)+2*v*tree[x].sum;
    tree[x].sum+=len*v;
}

inline void pushdown(int x,int l,int r) {
    if(abs(tree[x].add)>eps){
        int mid=(l+r)>>1;
        update(ls(x),mid-l+1,tree[x].add);
        update(rs(x),r-mid,tree[x].add);
        tree[x].add=0.0;
    }
}

void build(int x,int l,int r){
    if(l==r){
        tree[x].sum=a[l];
        tree[x].sqr=square(a[l]);
        return;
    }
    int mid=(l+r)>>1;
    build(ls(x),l,mid);
    build(rs(x),mid+1,r);
    pushup(x);
}

void modify(int x,int l,int r,int ml,int mr,double v){
    if(ml<=l&&r<=mr){
        update(x,r-l+1,v);
        return;
    }
    int mid=(l+r)>>1;
    pushdown(x,l,r);
    if(ml<=mid)
        modify(ls(x),l,mid,ml,mr,v);
    if(mid<mr)
        modify(rs(x),mid+1,r,ml,mr,v);
    pushup(x);
}

double querySum(int x,int l,int r,int ql,int qr){
    if(ql<=l&&r<=qr)
        return tree[x].sum;
    int mid=(l+r)>>1;
    double res=0.0;
    pushdown(x,l,r);
    if(ql<=mid)
        res+=querySum(ls(x),l,mid,ql,qr);
    if(mid<qr)
        res+=querySum(rs(x),mid+1,r,ql,qr);
    return res;
}

double querySquare(int x,int l,int r,int ql,int qr){
    if(ql<=l&&r<=qr)
        return tree[x].sqr;
    int mid=(l+r)>>1;
    double res=0.0;
    pushdown(x,l,r);
    if(ql<=mid)
        res+=querySquare(ls(x),l,mid,ql,qr);
    if(mid<qr)
        res+=querySquare(rs(x),mid+1,r,ql,qr);
    return res;
}

inline double queryAvarage(int l,int r){
    return querySum(1,1,n,l,r)/(r-l+1);
}

inline double queryVariance(int l,int r){
    double sum=querySum(1,1,n,l,r);
    double sqr=querySquare(1,1,n,l,r);
    double avg=sum/(r-l+1);
    double res=(sqr+(r-l+1)*square(avg)-2*avg*sum)/(r-l+1);
    return res;
}

int main(int argc, char *argv[]) {
#if !defined(ONLINE_JUDGE) && FILE_IO
    freopen("project.in", "r", stdin);
    freopen("project.out", "w", stdout);
    clock_t stime = clock();
#endif
    // [ Codes ] ////////////////////////
    
    scanf("%d%d",&n,&m);
    for(int i=1;i<=n;++i) scanf("%lf",&a[i]);
    build(1,1,n);
    while(m--){
        int op,x,y; double k;
        scanf("%d%d%d",&op,&x,&y);
        if(op==1){
            scanf("%lf",&k);
            modify(1,1,n,x,y,k);
        }else if(op==2){
            printf("%.4lf
",queryAvarage(x,y));
        }else{
            printf("%.4lf
",queryVariance(x,y));
        }
    }

    // [ Codes ] ////////////////////////
#if !defined(ONLINE_JUDGE) && FILE_IO
    clock_t etime = clock();
    printf("
----------------------------
");
    printf("Time : %dms
", etime - stime);
    if (etime - stime >= TIME_LIMIT) {
        printf("<<< Warning >>> Time Limited Exceeded
");
    }
#endif
    return 0;
}

//
//        ^ y
//       1|                            y = sin x (0 < x < 2π)
//  - - - | - - - + - - - - - - - - - - - - - -
//        |  +         +
// -------+---------------+---------------+----------> x
//       O|               π  +         +  2π
//  - - - | - - - - - - - - - - - + - - - - - -
//      -1|
//        |
//

Copyright ? 2019 ctjcalc，转载请注明URL，并给出原文链接，谢谢。