If *p* is very small, directly computing log(1+*p*) can be inaccurate. Numerical libraries often include a function `log1p`

to compute this function. The need for such a function is easiest to see when p is extremely small. If *p* is small enough, 1 + *p* = 1 in machine arithmetic and so `log(1+p)`

returns `log(1)`

which equals zero. All precision is lost. If *p* is small but not so extremely small, direct computation still loses precision, just not as much.

We can avoid the loss of precision by using a Taylor series to evaluate log(1 + *p*). For small *p*, log(1+*p*) ≈ *p* – *p*^{2}/2 with an error roughly equal to *p*^{3}/3. So if |*p*| is less than 10^{-4}, the error in approximating log(1+*p*) by *p* – *p*^{2}/2 is less than 10^{-12}. So we could implement the function `LogOnePlusX`

as follows.

The following code first appeared in A literate program to compute the inverse of the normal CDF.

#include <cmath> #include <sstreamh> #include <iostreamh> #include <stdexcepth> // compute log(1+x) without losing precision for small values of x double LogOnePlusX(double x) { if (x <= -1.0) { std::stringstream os; os << "Invalid input argument (" << x << "); must be greater than -1.0"; throw std::invalid_argument( os.str() ); } if (fabs(x) > 1e-4) { // x is large enough that the obvious evaluation is OK return log(1.0 + x); } // Use Taylor approx. log(1 + x) = x - x^2/2 with error roughly x^3/3 // Since |x| < 10^-4, |x|^3 < 10^-12, relative error less than 10^-8 return (-0.5*x + 1.0)*x; }

This code is in the public domain. Do whatever you want to with it, no strings attached.