计算没有上溢或下溢的 3D 欧几里德距离

Question

计算没有上溢或下溢的 3D 欧氏距离

你好，

我编写了代码来计算两个 vector 之间的 3D 欧几里德距离。
我知道这是一个非常常见的操作，但是我的情况有点特殊，因为我们需要在计算平方根之前添加一个标量。

以下是 C 代码中的例程:

void calculateAdjustedDistances3D(float *Ax, float *Ay, float *Az,
                                  float *Bx, float *By, float *Bz,
                                  float scalar, float *distances, int N)
{
    int i;

    for (i = 0; i < N; i++) {

        float dx = Ax[i] - Bx[i];
        float dy = Ay[i] - By[i];
        float dz = Az[i] - Bz[i];

        float dx2 = dx * dx;
        float dy2 = dy * dy;
        float dz2 = dz * dz;

        // potential for overflow/underflow
        float adjustedSquaredDistance = dx2 + dy2 + dz2 + scalar;

        distances[i] = sqrtf(adjustedSquaredDistance);
    }
}

对于我的应用程序，输入值的范围可以非常小，也可以非常大。 因此我现在正在考虑在计算这些距离时是否有必要防止上溢和下溢。 我知道有一些技术可以消除上溢/下溢的危险，但代价是使代码稍微慢一些。

例如，hypot()函数通常用于解决这个问题，但在这种情况下不能使用它，因为在计算平方根之前要添加标量。

如何重写代码以减轻或理想地消除计算中溢出和下溢的可能性？

Answer 1

对于标量的正值，人们也许可以将其解释为“额外维度”并通过以下方式调用hypot函数

distances[i] = hypot(hypot(hypot(dx, dy), dz), sqrt(scalar))

编辑:

为了避免调用函数 hypot，可以遵循 BLAS function snrm2 的实现它应该以“稳健”的方式计算 vector 2-范数:

float scale, ssq, ax;
const float x[4] = {dx, dy, dz, sqrt(scalar)};

scale = 0;
ssq = 1;

for(int j=0; j<4; j++){
    if(x[j] != 0){
        ax = fabs(x[j]);
        if(scale < ax){
            ssq = 1 + ssq*(scale/ax)*(scale/ax);
            scale = ax;
        }else{
            ssq += (ax/scale)*(ax/scale);
        }
    }
}
distances[i] = scale*sqrt(ssq);