If *x* is very small, directly computing `exp(`

can be inaccurate. Numerical libraries often include a function *x*) - 1`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(`

returns 0 even though the correct result is positive. All precision is lost. If *x*) - 1*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* + *x*^{2}/2 + *x*^{3}6 + …

If |*x*| < 10^{-5}, the error in approximating exp(*x*) – 1 by *x* + *x*^{2}/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

#include

// 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"; } [/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#