稳健准确地计算两个浮点数商的自然对数

Question

计算log (a/b)时的一个明显问题，其中a和b是给定精度（这里称为本机精度）的两个非零正有限浮点操作数, 是商 a/b 可能无法表示为该精度的浮点数。此外，当源操作数的比例接近于一时，精度会下降。

这可以通过暂时切换到更高精度的计算来解决。但是这种更高的精度可能并不容易获得，例如当本机精度为 double 并且 long double 只是映射到 double 时。使用更高精度的计算也可能对性能产生非常显着的负面影响，例如在 GPU 上，float 计算的吞吐量可能比 double 计算的吞吐量高 32 倍。

当 a 和 b 彼此接近时，可以决定使用 quotient rule of logarithms to compute log (a/b) as log(a) - log(b), but this exposes the computation to the risk of subtractive cancellation，从而导致非常大的误差。

如何准确计算两个给定精度的浮点数的商的对数，例如误差小于 2 ulps 且稳健，即在中间计算中没有下溢和上溢，无需诉诸高于本机精度的计算？

Answer 1

到目前为止，我确定的最佳方法区分了三种情况，这三种情况基于较大的源操作数除以较小的源操作数的商。这个比率告诉我们操作数相距多远。如果它太大以至于超过了本机精度的最大可表示数，则必须使用商规则，结果计算为 log(a) - log(b)。如果比例接近于一，计算应该利用log1p()函数来提高精度，计算结果为log1p ((a - b) / b)。 Sterbenz Lemma 表明 2.0 是一个很好的切换点，因为如果比率 ≤ 2 将准确计算 a-b。对于所有其他情况，直接计算 log (a/b) 可以使用。

下面，我展示了接受 float 个参数的函数的此设计的实现。使用 float 可以更轻松地评估准确性，因为这允许对可能的测试用例进行更密集的采样。显然，整体准确性将取决于数学库中 logf() 和 logpf() 的执行质量。使用具有几乎正确舍入的函数的数学库（logf() 中的最大误差 < 0.524 ulp，log1pf() 中的最大误差 < 0.506 ulp），在 log_quotient() 中观察到的最大误差是< 1.5 ULPS。使用具有忠实舍入函数实现的不同库（logf() 中的最大误差 < 0.851 ulp，log1pf() 中的最大误差 < 0.874 ulp），在 log_quotient() 中观察到的最大误差< 1.7 ulps。

#include <float.h>
#include <math.h>

/* Compute log (a/b) for a, b ∈ (0, ∞) accurately and robustly, i.e. avoiding
   underflow and overflow in intermediate computations. Using a math library 
   that provides log1pf() and logf() with a maximum error close to 0.5 ulps,
   the maximum observed error was 1.49351 ulp.
*/
float log_quotient (float a, float b)
{
    float ratio = fmaxf (a, b) / fminf (a, b);
    if (ratio > FLT_MAX) {
        return logf (a) - logf (b);
    } else if (ratio > 2.0f) {
        return logf (a / b);
    } else {
        return log1pf ((a - b) / b);
    }
}

稳健准确地计算两个浮点数商的自然对数

Robustly and accurately computing natural logarithm of quotient of two floating-point numbers

c

algorithm

floating-point

logarithm