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

Posted 2023-02-16

【中文标题】计算没有上溢或下溢的 3D 欧几里得距离

【英文标题】：Calculating 3D Euclidean Distance without overflows or underflows

【发布时间】：2015-11-13 22:46:43

【问题描述】：

你好，

我编写了代码来计算两个向量之间的 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()函数通常用于解决这个问题，但是由于在计算平方根之前添加了标量，因此不能在这种情况下使用它。

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

【问题讨论】：

“很小，也很大”。实际限制是多少？这些限制会导致问题吗？似乎您是在尝试解决问题而不首先确定问题是否存在。

没有限制。可以可能使用全范围的 32 位浮点值。确实，在一般情况下，输入的值不会引起问题，但它们可能会引起问题。为了编写健壮的代码，我试图保护自己免受那些不太可能但可能发生的情况。

将 dx,dy,dz,dx2 等声明为 double 并使用 sqrt 代替 sqrtf。

【参考方案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);

【讨论】：

是的，这会起作用，但是它会非常慢，因为我必须计算 4 个平方根。我认为，为了保持以相当快的速度运行代码，应该避免使用 hypot() 解决方案。好主意……

在这种情况下，也许可以从 BLAS 函数snrm2 中获得灵感，它似乎执行某种缩放以减轻浮点错误（而且，它只计算一个平方根）-@987654322 @

是的，这就是我要找的。谢谢。