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

Question

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

你好，

我编写了代码来计算两个向量之间的 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

对于scalar 的正值，也许可以将其解释为“额外维度”并以下列方式调用hypot 函数

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

编辑：

为了避免调用函数hypot，可以遵循BLAS function snrm2 的实现，它应该以“稳健”的方式计算向量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);