实现傅斌的近似和算法

Question

我正在尝试实现 Bin Fu's approximate sum algorithm 用真实的语言更好地了解它的工作原理。

In a nutshell，这是一种有效计算 $(1+\epsilon)$-bounds 值的算法 $s(x)=\sum_{i=1}^n x_i$ 其中 $x$ 是已排序浮点数的向量。

但是，我一定做错了，因为运行算法会导致错误 （我对伪算法语言也不是很精通，而且像数组边界检查这样的东西似乎隐含在这段代码中）。

这是我到目前为止的非工作代码，欢迎任何有关该问题的提示/帮助-我不知道语言，我只是使用 R，因为它是 1-index（算法是 1- index) 开源解释语言：

ApproxRegion<-function(x,n,b,delta){
    if(x[n]<b)  return(NULL)
    if(x[n-1]<b)    return(c(n,n))
    if(x[1]>=b) reurn(c(1,n))
    m1<-2
    while(x[n-m1**2+1]>=b)  m1<-m1**2
    i<-1
    m1<-m1
    r1<-m1
    while(m1>(1+delta)){
        m1<-sqrt(m1)
        if(x[n-floor(m1*r1)+1]>=b){
            r1<-m1*r1
        } else {
            r1=r1
        }
        i=i+1
    }
    return(c(n-floor(r1*m1)+1,n))
}       
ApproxSum<-function(x,n,epsilon){
    if(x[n]==0) return(0)
    delta<-3*epsilon/4
    r1p<-n
    s<-0
    i<-1
    b1<-x[n]/(1+delta)
    while(b1>=((delta*x[n])/(3*n))){
        Ri<-ApproxRegion(x=x,n=r1p,b=b1,delta=delta)
        r1p<-Ri[1]-1
        b1<-x[r1p]/(1+delta)
        s1<-(Ri[2]-Ri[1]+1)*b1
        s<-s+s1
        i<-i+1
    }
    return(s)
}
n<-100;
x<-sort(runif(n));
ApproxSum(x=x,n=length(x),epsilon=1/10);
sum(x)

作者提到了一个c++版本，但我在网上找不到（前面的任何帮助也很好）。

Modo：我把问题放在这里（而不是理论上的 CS stackexchange 站点），因为它是关于实现问题。随意移动。

编辑

原始代码有一个“毛茸茸”的退出条件（x[i]=$-\infty$ for $i\leq 0$）。 按照 Martin Morgan 的建议，我用适当的 break 替换了 this 的出现，产生了以下代码：

ApproxRegion<-function(x,b,delta,n){
    if(n<=1)            return(NULL)
    if(x[n]<b)          return(NULL)
    if(x[n-1]<b)            return(c(n,n))
    if(x[1]>=b)         return(c(1,n))
    m<-2
    xit<-0
    while(!xit){
        if(n-m**2+1<1)      break
        if(x[n-m**2+1]<b)   break
        m<-m**2
    }
    i<-1
    r<-m
    while(m>=(1+delta)){
        m<-sqrt(m)  
        if(n-floor(m*r)+1>0){
            if(x[n-floor(m*r)+1]>=b)    r=m*r   
        }   
        i<-i+1      
    }
    return(c(n-floor(m*r)+1,n))
}       
ApproxSum<-function(x,n,epsilon){
    if(x[n]==0) return(0)
    delta=3*epsilon/4
    rp<-n
    s<-0
    i<-1
    b<-x[n]/(1+delta)
    while(b>=delta*x[n]/(3*n)){
        R<-ApproxRegion(x,b,delta,rp)
            if(is.null(R))  break   
        if(R[1]<=1) break
        rp<-R[1]-1
        b<-x[rp]/(1+delta)
        si<-(R[2]-R[1]+1)*b
        s<-s+si
        i<-i+1
    }
    return(s)
}

现在，它可以工作了：

n<-100;
set.seed(123)
x<-sort(runif(n));
ApproxSum(x=x,n=length(x),epsilon=1/10);
sum(x)

Answer 1

作为部分答案...存在算法未明确处理的边缘条件。例如在ApproxRegion 中，需要注意n = 0（返回值应该为NULL？）或1（c（n，n）？）否则第一个或第二个条件x[n] < b，x[n - 1] < b 不会评估为预期（例如，x[0] 返回数字（0））。同样，循环中的测试必须防范m1**2 > n + 1，否则下标为负数。

我认为ApproxSum 中也存在类似问题，尤其是当ApproxRegion 返回时，例如c(1, 1)（因此 r1p == 0，b1 = integer()）。看到更新的实现会很有趣。