log (a/b)时的一个明显问题，其中a和b是给定精度(此处称为本机精度(的两个非零正有限浮点操作数，是商a/b可能无法表示为该精度的浮点数。此外，当源操作数的比率接近单位时，精度将丢失。这可以通过暂时切换到更高精度的计算来解决。但是这种更高的精度可能并不容易获得，例如，当本机精度double并且long double只是映射到double时。使用更高精度的计算也可能对性能产生非常显着的负面影响，例如在 GPU 上，float计算的吞吐量可能比double计算的吞吐量高出 32 倍。
人们可以决定使用对数的商规则来计算log (a/b)作为log(a) - log(b)，但是当a和b彼此接近时，这会使计算面临减法取消的风险，从而导致非常大的误差。
如何准确计算给定精度的两个浮点数的商的对数，例如误差小于 2 ulps，又稳健地计算，即在中间计算中没有下溢和溢出，而不诉诸高于本机精度的计算？
到目前为止，我确定的最佳方法是区分三种情况，这些情况基于较大源操作数的商除以较小的源操作数。这个比率告诉我们操作数之间的距离有多远。如果它太大以至于超过本机精度的最大可表示数，则必须使用商规则，并将结果计算为 log(a) - log(b) 。如果比率接近单位，则计算应利用函数log1p()来提高精度，将结果计算为log1p ((a - b) / b)。Sterbenz 引理表明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);
}
}