What is a nice way to compute $f(x) = x / (\exp(x) - 1)$?

Question

I want it to be stable near $f(0) = 1$. Is there a nice function that does this already, like maybe a hyperbolic trig function or something like expm1, or should I just check if $x$ is near zero and then use a polynomial approximation?

Answer 1

If you don't want to use the expm1() function for some reason, one possibility, detailed in Higham's book, is to let $u=\exp\,x$ and then compute $\log\,u/(u-1)$. The trick is attributed to Velvel Kahan.

Answer 2

If your system provides expm1(x), they should have worried about errors and stability at least as well as you can. It is a better question if that is not available. Wolfram Alpha gives $1-\frac {x}2+\frac {x^2}{12}$ for the second order series around $0$, so you could check if $x$ is close to zero and use that.

Answer 3

You mention using hyperbolic functions. You might try $$ \frac{x}{\exp(x)-1}=\frac{x/2}{\exp(x/2)\sinh(x/2)} $$ This loses no precision if the $\sinh$ is computed to full precision by the underlying system.

Note that $\mathrm{expm1}(x)=2\exp(x/2)\sinh(x/2)$.

Example: 15 digit calculations

$x=.00001415926535897932$: $$ \begin{align} \frac{x}{\exp(x)-1} &=\frac{.00001415926535897932}{1.000014159365602-1}\\ &=\color{#C00000}{.9999929203}73447 \end{align} $$ $$ \begin{align} \frac{x/2}{\exp(x/2)\sinh(x/2)} &=\frac{.00000707963267948966}{1.000005000012500\cdot.00000707963267954879}\\ &=\color{#C00000}{.999992920384028} \end{align} $$

Answer 4

Consider the Bernoulli numbers, defined by the recursive formula:

$$B_0=1$$

$$\sum_{k<n} {n\choose k }B_k=0\text{ ; } n\geq 2$$

This gives the sequence:

$$\{B_n\}_{n\in \Bbb N}=\left\{ 1,-\frac 1 2,\frac 1 6 ,0,\frac 1 {30},0,\dots\right\}$$

It's generating function is

$$ \sum_{n=0}^\infty B_n \frac{x^n}{n!}=\frac{x}{e^x-1}$$

It's first few terms are

$$1-\frac x 2 +\frac {x^2}{12}-\frac{x^4}{720}+\cdots$$

The numbers' denominators grow pretty fast, so you should have no problem with convergence: in fact, the function is essentialy $=-x$ for large negative $x$ and $=0$ for large positive $x$, so a few terms should suffice the "near origin" approximation.

Answer 5

Just posting @user856's comment as an answer, since it's perfectly good; nothing complicated is needed here:

"Once you use expm1 to compute exp(x)−1, there's no further loss of significance in dividing x by it."

So:

if (x == 0)
  return 1;
else
  return x / expm1(x);

Note, however, that x/expm1(x) yields 1 for sufficiently small numbers that are still much larger than zero; this can be seen from its power series expansion (from worlframalpha.com):

$$1 - x/2 + x^2/12 - x^4/720 + O(x^6)$$

So the following will work too:

if (1 - x/2 == 1)
  return 1;
else
  return x / expm1(x);