看起来不必要的数学函数

微信公众号：Julia语言
每周一三五更新Julia语言；
每周二四六更新Python进阶；

Math library functions that seem unnecessary

https://www.johndcook.com/blog/2010/06/07/math-library-functions-that-seem-unnecessary/

This post will give several examples of functions include in the standard C math library that seem unnecessary at first glance.

Function log1p(x) = log(1 + x)

The function log1p computes log(1 + x). How hard could this be to implement?

julia> log(1+9)
2.302585092994046

julia> log1p(9)
2.302585092994046

But wait. What if x is very small? If x is 10^-16, for example, then on a typical system 1 + x = 1 because machine precision does not contain enough significant bits to distinguish 1 + x from 1. (For details, see Anatomy of a floating point number.) That means that the code log(1 + x) would first compute 1 + x, obtain 1, then compute log(1), and return 0. But log(1 + 10^-16) ≈ 10^-16. This means the absolute error is about 10-16 and the relative error is 100%. For values of x larger than 10-16 but still fairly small, the code log(1 + x) may not be completely inaccurate, but the relative error may still be unacceptably large.

Fortunately, this is an easy problem to fix. For small x, log(1 + x) is approximately x. So for very small arguments, just return x. For larger arguments, compute log(1 + x) directly. See details here.

Why does this matter? The absolute error is small, even if the code returns a zero for a non-zero answer. Isn’t that OK? Well, it might be. It depends on what you do next. If you add the result to a large number, then the relative error in the answer doesn’t matter. But if you multiply the result by a large number, your large relative error becomes a large absolute error as well.

Function expm1(x) = exp(x) – 1

Another function that may seem unnecessary is expm1. This function computes e^x – 1. Why not just write this?

julia> exp(8)-1
2979.9579870417283

julia> expm1(8)
2979.9579870417283

That’s fine for moderately large x. For very small values of x, exp(x) is close to 1, maybe so close to 1 that it actually equals 1 to machine precision. In that case, the code exp(x) - 1 will return 0 and result in 100% relative error. As before, for slightly larger values of xthe naive code will not be entirely inaccurate, but it may be less accurate than needed. The solution is to compute exp(x) – 1 directly without first computing exp(x). The Taylor series for exp(x) is 1 + x + x²/2 … So for very small x, we could just return x for exp(x) – 1. Or for slightly larger x, we could return x + x²/2. (See details here.)

欢迎关注微信公众账号Julia语言.jpg

点击阅读原文可查看历史文章