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

If x is very small, directly computing exp(x) - 1 can be inaccurate. Numerical libraries often include a function expm1 to compute this function. The need for such a function is easiest to see when x is extremely small. If x is small enough, exp(x) = 1 in machine arithmetic and so exp(x) - 1 returns 0 even though the correct result is positive. All precision is lost. If x 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 exp(x):

exp(x) = 1 + x + x2/2 + x36 + ...

If |x| < 10-5, the error in approximating exp(x) - 1 by x + x2/2 is on the order of 10-15 and so the relative error is on the order of 10-10 or better. If we compute exp(10-5) - 1 directly, the absolute error is about 10-16 and so the relative error is about 10-11. So by using the two-term Taylor approximation for |x| less than 10-5 and the direct method for |x| larger than 10-5, we obtain at least 10 significant figures for all inputs.

 

#include <cmath>
#include <iostream>

// Compute exp(x) - 1 without loss of precision for small values of x.
double expm1(double x)
{
	if (fabs(x) < 1e-5)
		return x + 0.5*x*x;
	else
		return exp(x) - 1.0;
}

void testExpm1()
{
    // Select a few input values
    double x[] = 
    {
        -1, 
        0.0, 
        1e-5 - 1e-8, 
        1e-5 + 1e-8,
		0.5
    };

    // Output computed by Mathematica
    // y = Exp[x] - 1
    double y[] = 
    { 
       -0.632120558828558, 
        0.0, 
        0.000009990049900216168, 
        0.00001001005010021717, 
        0.6487212707001282 
    };

	int numTests = sizeof(x)/sizeof(double);

    double maxError = 0.0;
    for (int i = 0; i < numTests; ++i)
    {
        double error = fabs(y[i] - expm1(x[i]));
        if (error > maxError)
            maxError = error;
    }

	std::cout << "Maximum error: " << maxError << "\n";
} 

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