时间戳

首页 > 解决方案 > 稳健准确地计算两个浮点数的商的自然对数

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

问题描述

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

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

可以决定使用对数的商规则来计算log (a/b)as log(a) - log(b)，但这会使计算在和彼此接近时面临减法抵消的风险，从而导致非常大的误差。ab

如何既准确地计算给定精度的两个浮点数的商的对数，例如误差小于 2 ulps，又如何稳健地计算，即中间计算中没有下溢和上溢，而不诉诸高于原生精度计算？

标签： calgorithmfloating-pointlogarithm

解决方案

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

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

#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);
    }
}

推荐阅读