13
$\begingroup$

The function $f \colon x \mapsto (e^x-1)/x$ has singularity near $x = 0$. That singularity can be lifted, though: for $x = 1$, one should have $f(x) = 1$, since $$ e^x = \sum_{k=0} \frac{x^k}{k!} $$ and thus $$ (e^x - 1)/x = \sum_{k=1} \frac{x^{k-1}}{k!} $$ However, the form $(e^x-1)/x$ is not only not defined at $x = 0$, it is also numerically unstable in the vicinity of that point; in order to evaluate $f(x)$ for very small $x$ numerically, one could use a Taylor expansion, i.e. a truncation of the afore-mentioned power series.

Q: Does the function $f$ have a name? In other words, is this a common problem?

Q: Is anyone aware of a C/C++ library that handles this situation nicely, i.e. uses the Taylor expansion of an appropriate degree near 0 and the other representation away from zero?

Improve this question
$\endgroup$
19
$\begingroup$

Possibly one could start with the function $\mathtt{expm1}$ which is part of the C99 standard, and calculates $e^x-1$ accurately near $x=0$.

Improve this answer
$\endgroup$
0
17
$\begingroup$

This is an instance of cancellation error. The C standard library (as of C99) includes a function called expm1 that avoids this problem. If you use expm1(x) / x instead of (exp(x) - 1.0) / x, you won't experience this issue (see graph below). <code>fabs(expm1(x) / x - (exp(x) - 1.0) / x)</code>

The details and solution of this particular problem are discussed at length in Section 1.14.1 of Accuracy and Stability of Numerical Algorithms. The same solution is also explained in page 19 of W. Kahan's paper titled How Futile are Mindless Assessments of Roundoff in Floating-Point Computation?. The actual implementation of expm1 in the GNU C library is different from the approach described in the references above and is thoroughly documented in the source code.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ Thank you, that's just what I needed! Unfortunately, I can only accept one answer... $\endgroup$ – anonymous Aug 31 '13 at 0:02
  • $\begingroup$ Of course! No problem :-) $\endgroup$ – Juan M. Bello-Rivas Aug 31 '13 at 0:04
3
$\begingroup$

To answer your first question, no, the function doesn't have a name (at least not one that is widely known).

As others have mentioned, the best way to compute the function is to treat several special cases. This is how any library would compute the function.

  1. Case 0: x = 0, return 1.
  2. Case 1: $|x|<\delta$, return $1+x/2$. Choose $\delta$ empirically for your numerical precision and data type. For double, it should be about 2e-8. For float, it's about 5e-4.
  3. Case else: return expm1(x)/x.

You can be more sophisticated and special case more things with truncated Taylor series, but it's probably not worth it. In fact, it's not entirely clear that case 1 needs to be handled separately, since as k20 pointed out, the cancellation is safe. However, handling it separately would let me feel more confident about it.

Improve this answer
$\endgroup$
2
$\begingroup$

I remember this question having been asked earlier on this site, and surprisingly the answer is that you only need to special-case exact equality to zero. The errors cancel out near zero. I don't have the link.

Yeah this answer was completely wrong. I'm not sure why it was upvoted so much, probably because it was stated so authoritatively. I found the link I had in mind. It was on the math stackexchange here, not on the scicomp stackexchange. The expm1-free error cancellation formula is given in the answer by J.M. and uses a u = exp(x) transformation.

Improve this answer
$\endgroup$
2
  • $\begingroup$ + Right. If $x$ = infinitessimal $dx$, $(e^{dx} - 1)/dx$ = $(1+dx - 1)/dx$ = $1$. $\endgroup$ – Mike Dunlavey Aug 30 '13 at 21:00
  • 1
    $\begingroup$ I don't think this works for $dx$ small enough since then in floating point accuracy, $1+dx=1$ and the whole thing blows up. $\endgroup$ – Wolfgang Bangerth Aug 30 '13 at 21:48
0
$\begingroup$

To answer the first question and provide a (probably numerically inefficient) method for the second, note that this is the inverse of the generating function of the Bernoulli numbers.

Improve this answer
$\endgroup$
4
  • $\begingroup$ That's an interesting connection, thanks for pointing that out. Unfortunately, I believe the triple sum will render this prohibitively expensive. Moreover, it is not immediately clear where to truncate each sum to obtain the desired accuracy. $\endgroup$ – anonymous Feb 17 '17 at 12:43
  • $\begingroup$ @anonymous: Which tripple sum do you mean? You don't need the Bernoulli polynomials, only the Bernoulli numbers, and you can list those in advance. But yes, it's still not to be better than the Taylor series. $\endgroup$ – Nikolaj-K Feb 19 '17 at 2:46
  • $\begingroup$ You can compute those in advance if it's clear that you only need a fixed finite number for any input only, though. $\endgroup$ – anonymous Feb 19 '17 at 10:25
  • $\begingroup$ @anonymous: Well yes, just like you'd list the Taylor coefficients in advance. $\endgroup$ – Nikolaj-K Feb 24 '17 at 21:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.