一些数学库函数有存在的必要吗?
本文将给出一些在 C 数学库中乍一看似乎不必要存在的函数示例。
函数 log1p(x) = log(1 + x)
函数 log1p
计算 log(1 + x)
。这有多难实现呢?
1 | log(1 + x); |
好了。
等等。。。如果 x 非常小呢?例如,如果 x 是 $10^{-16}$,那么在典型系统上会出现 1 + x = 1
,因为机器精度缺乏足够的有效位来区分 1 + x
和 1(参见这里了解更多)。这意味着代码 log(1 + x)
会先计算 1 + x
,得到 1,然后计算 log(1)
,返回 0。但是 $log(1 + 10^{-16}) \approx 10^{-16}$
。这意味着绝对误差约 $10^{-16}$
,而 **相对误差为 100%**。对于那些比 $10^{-16}$
大但总体还是相当小的数来说,代码 log(1 + x)
可能并不完全不正确,但是相对误差还是大的不可接受。
幸运的是,这是一个很好解决的问题。对于小的 x 值,log(x + 1)
约等于 x。所以对于非常小的参数,只需要返回 x。对于大的参数,直接计算 log(1 + x)
。具体细节参见这里。
为什么这很重要?即使代码对于非零的答案直接返回零,其绝对误差也很小。这不是还可以吗?是的,它可能是可行的。这取决于你下一步做什么——如果将结果添加到大量数据中,那么答案中的相对误差无关紧要;但是如果将结果乘以大数,大的相对误差也会变成大的绝对误差。
函数 expm1(x) = exp(x) - 1
另一个似乎不必要的函数是 expm1
。这个函数计算 $e^x-1$
。为什么不写成这样?
1 | exp(x) - 1.0; |
对于适当的 x 来说这能很好的工作。但是,对于非常小的 x 值,exp(x)
接近于 1,这可能会导致太接近 1 以至于在机器精度上实际就是 1。在这种情况下,代码 exp(x) - 1
将返回 0 并产生 100% 的相对误差。和上面的一样,对于稍微大一点的 x 值,朴素代码不会完全错误,但可能不如所需的精确。解决方案是直接计算 exp(x) - 1
,而不是先计算 exp(x)
。exp(x)
的泰勒级数是 $1+x+\frac{x^2}{2}...$
,所以对于非常小的 x,exp(x) - 1
返回 x。对于稍大的 x,可以返回 $x + x^2/2$
。(具体细节参见这里)
函数 erf(x) 与 erfc(x)
C 数学库包含一对函数 erf
与 erfc
。结尾的 c 代表互补(complement),因为 erfc(x) = 1 - erf(x)
。 函数 erf(x)
被称为误差函数,实现并不容易。但是为什么还有一个单独的 erfc
函数?一旦你有 erf
的代码实现它不是很容易吗?对于某些 x 的值来说,这是真的。但是,如果 x 很大,则 erf(x)
接近 1,当结果很小且为正时,代码 1 - erf(x)
可能会返回 0。和上面的例子一样,朴素的实现会导致某些值精度的完全丧失,其他一些值精度也会部分丢失。
函数 $\Gamma(x)$
与 $log\Gamma(x)$
标准 C 数学库有两个与 gamma 函数相关的函数:返回 $\Gamma(x)$
的 tgamma
和返回 $log\Gamma(x)$
的 lgamma
。为什么有两个?为什么后者不能使用前者的对数?gamma 函数的增长非常快。对于一般大小的参数,其值超过计算机数值的容量。 有时你可能需要这些天文数字大小的值作为中间结果。你需要一个适当大小的数作为两个非常大的值的比值(这里有一个例子)。在这种情况下,只需要将 lgamma
值相减,而不是用 tgamma
的值相比。(为什么函数叫做 tgamma 而不是 gamma?请参阅这篇文章的最后一段与第一条评论)
总结
标准 C 数学库汲取了许多经验。对于一些参数来说有些函数是不必要的,但是,对于另一些参数来说它们是不可或缺的。