c++ - 实现 Bin Fu 的近似求和算法

Question

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

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

但是，我一定做错了什么，因为运行算法会导致错误(我也不是很精通伪算法语言，并且数组绑定(bind)检查等一些东西似乎隐含在这段代码中)。

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

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 的建议，我用适当的中断替换了这种情况，产生了以下代码:

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] 返回 numeric(0))。同样，循环中的测试必须防范 m1**2 > n + 1 , 否则你将下标为负数。

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