代码之家 › 专栏 › 技术社区 › njuffa

两个浮点数商的自然对数的鲁棒精确计算

logarithm floating-point algorithm c

njuffa · 技术社区 · 7 年前

计算时的一个明显问题 log (a/b) ,在哪里 a 和 b 是两个给定精度的非零正有限浮点操作数(这里称为本机精度),是指商 a/b 在该精度中不能表示为浮点数。此外,当源操作数的比值接近于一时,精度会降低。

这可能通过暂时切换到更高精度的计算来解决。但这种更高的精度可能不容易获得,例如当本机精度为 double 和 long double 只需映射到 双重的 . 使用更高精度的计算也可能对性能产生非常显著的负面影响,例如对GPU的吞吐量 float 计算可能比 双重的 计算。

我们可以决定使用 quotient rule 要计算的对数 日志(A/B) 作为 log(a) - log(b) 但这就暴露了计算的风险 subtractive cancellation 什么时候? 一 和 乙 彼此接近,造成很大的误差。

给定精度的两个浮点数的商的对数如何既能精确地计算,如误差小于2 ulps,又能可靠地计算,即在中间计算中没有下溢和溢出,而不采用高于本机精度的计算?

1 回复 | 直到 7 年前

njuffa 7 年前

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

下面,我展示了一个接受函数的设计的实现 float 参数。使用 浮动 这样可以更容易地评估准确性,因为这样可以对可能的测试用例进行更密集的采样。显然,总体的准确性将取决于实施的质量 logf() 和 logpf() 在数学图书馆。使用具有几乎正确四舍五入的函数的数学库(在 Logf() <0.524 ulp,最大误差 log1pf() <0.506 ulp),在 log_quotient() 是1.5 ulps。使用不同的库,并对函数进行精确的四舍五入实现(在 Logf() <0.851 ulp,最大误差 日志1pf() <0.874 ulp),在 对数商() 是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);
    }
}