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 – p2/2 with an error roughly equal to p3/3. So if |p| is less than 10-4, the error in approximating log(1+p) by p – p2/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.