Stand-alone C++ code for log(1+x)

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 with an error roughly equal to p³/3. So if |p| is less than 10^-4, the error in approximating log(1+p) by p - p²/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 <sstream>
#include <iostream>
#include <stdexcept>

// 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.

Other versions of the same code: Python, C#

Stand-alone numerical code

 

Home