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

#include

#include

#include

// 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;
}
[/code]
This code is in the public domain. Do whatever you want to with it, no strings attached.
Other versions of the same code: Python, C#