计算没有上溢或下溢的 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);
    }
}

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

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

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

Answer 1

对于 scalar 的正值，可以将其解释为 "extra dimension" 并按以下方式调用 hypot 函数

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

编辑：

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

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

Calculating 3D Euclidean Distance without overflows or underflows

c

3d

floating-accuracy

euclidean-distance