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?
log(1 + x);
Done.
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 ex – 1. Why not just write this?
exp(x) - 1.0;
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 x the 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 + x2/2 … So for very small x, we could just return x for exp(x) – 1. Or for slightly larger x, we could return x + x2/2. (See details here.)
Functions erf(x) and erfc(x)
The C math library contains a pair of functions erf
and erfc
. The “c” stands for “complement” because erfc(x) = 1 – erf(x). The function erf(x) is known as the error function and is not trivial to implement. But why have a separate routine for erfc
? Isn’t it trivial to implement once you have code for erf
? For some values of x, yes, this is true. But if x is large, erf(x) is near 1, and the code 1 - erf(x)
may return 0 when the result should be small but positive. As in the examples above, the naive implementation results in a complete loss of precision for some values and a partial loss of precision for other values.
Functions Γ(x) and log Γ(x)
The standard C math library has two functions related to the gamma function: tgamma
that returns Γ(x) and lgamma
that return log Γ(x). Why have both? Why can’t the latter just use the log of the former? The gamma function grows extremely quickly. For moderately large arguments, its value exceeds the capacity of a computer number. Sometimes you need these astronomically large numbers as intermediate results. Maybe you need a moderate-sized number that is the ratio of two very large numbers. (See an example here.) In such cases, you need to subtract lgamma
values rather than take the ratio of tgamma
values. (Why is the function called “tgamma” and not just “gamma”? See the last paragraph and first comment on this post.)
Conclusion
The standard C math library distills a lot of experience. Some of the functions may seem unnecessary, and so they are for some arguments. But for other arguments these functions are indispensable.
More numerical posts