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

Question

computating一个明显的问题时log (a/b),其中a和b是给定的精度两个非零正有限浮点操作数(这里称为本机精度),是商a/b可能无法表示为在精度的浮点数.此外,当源操作数的比率接近于1时,精度将会丢失.

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

人们可以决定使用商法则对数计算的log (a/b)作为log(a) - log(b),但是这暴露了计算到的风险消减取消时a和b接近对方,产生非常大的误差.

如何精确计算给定精度的两个浮点数的商的对数,例如误差小于2 ulps并且鲁棒地计算,即在中间计算中没有下溢和溢出,而不需要求高于本机精度计算？

Answer 1

到目前为止,我已经确定的最佳方法区分了三种基于较大源操作数除以较小源操作数的商的情况.这个比率告诉我们操作数有多远.如果它超过本机精度的最大可表示数,则必须使用商规则,并将结果计算为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 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);
    }
}