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

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?

• Once you use expm1 to compute $\exp(x)-1$, there's no further loss of significance in dividing $x$ by it, is there? – Rahul Aug 21 '12 at 22:13
• @Rahul: You are suggesting to use x/expm1(x) but with a single extra check for 0/0 when x is exactly 0? – martin Aug 21 '12 at 22:29
• Suggesting, yes, (though to be clear, not with a great deal of authority). – Rahul Aug 21 '12 at 23:06
• Why not Bernoulli numbers? – Pedro Tamaroff Aug 22 '12 at 0:46
• $$f(x)={e^{-x/2}\over{\rm sinc}(i x/2)}\ .$$ – Christian Blatter Aug 22 '12 at 10:09

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.

• This is pretty cool and in fact it's what I meant by "check if x is near zero and then use a polynomial approximation." – martin Aug 22 '12 at 1:57
• @kxx, you might even consider building a Padé approximant from the series expansion Peter has shown you. For instance, the series expansion of $$\frac{(x-4)^2+4}{\frac{(x+3)^2}{3}+17}$$ matches the original function up to the fourth series term. – J. M. is a poor mathematician Aug 22 '12 at 11:34
• Dear Pedro - Please forgive this if it's a dopy question. Would you please show how to derive your recursive formula - $\sum_{k<n} =0$ from the generating function relation, which is my typical starting point. Thanks. Regards, – user12802 Jun 10 '15 at 21:32

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.

• I am totally OK with using expm1(). – martin Aug 21 '12 at 23:44
• Does this give any better precision? $u-1=\exp(x)-1$ would seem to lose precision near $x=0$. The computation of $\exp(x)$ is the place where precision is lost. Am I missing something? – robjohn Aug 22 '12 at 0:54
• @robjohn: Neither $u-1$ nor $\log u$ is very accurate for $u \sim 1$. But if you follow the link, it asserts the errors compensate for each other. – Excluded and Offended Aug 22 '12 at 1:54
• @rob, yes, there is an example there in Higham's book that shows how the effect of the errors in computing the numerator and denominator cancel out. I'll try to find the original ref by Kahan, and will link to it in the answer if I manage to do find the ref. – J. M. is a poor mathematician Aug 22 '12 at 11:36
• @J.M.: Yes, I believe this shows it: $$\frac{\log(\exp(x)+\delta)}{\exp(x)+\delta-1} -\frac{\log(\exp(x)}{\exp(x)-1}\sim\frac12\delta$$ for small $x$ and $\delta$. Very nice. – robjohn Aug 22 '12 at 12:24

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}

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.

• My system provides expm1. So you are seconding Rahul's suggestion to use x/expm1(x) with the check for x exactly equal to zero? – martin Aug 21 '12 at 23:46
• @kxx: that's what I would do. You could try some cases, comparing against Alpha (which will give arbitrary precision), for confirmation. – Ross Millikan Aug 21 '12 at 23:57

You can construct Lagrange interpolation polynomial, of degree $n$, given $n$ a positive integer. You should first extract values of $f (x) = \frac {x} {e^x - 1}$ at any $n + 1$ points you like, say, $x_0, x_1, \cdots, x_n$. Then $$L (x) := \sum_{0 \leqslant k \leqslant n} f (x_k) \prod_{0 \leqslant j \leqslant n} \frac {x - x_j} {x_k - x_j},$$ where $j \neq k$, is the Lagrange polynomial that interpolates $f (x)$ at the points $$(x_0, f (x_0)), (x_1, f (x_1)), \cdots, (x_n, f (x_n)).$$

But even if Lagrange polynomial is a very good approximation of the function that it interpolates, its analytic behaviour at turn points may radically differ from that of the original function. That is why, you may want to interpolate $f (x)$ with such a polynomial, that not only its values at given points overlap with those of the function, but also the values of any of its derivatives at the given points overlap with those of $f (x)$. Such is Taylor polynomial, i.e., Taylor series that is truncated at any given place (in our case, after the $(n + 1)$st term): $$\sum_{0 \leqslant k \leqslant n} \frac {f^{(k)}(a)}{k!} \, (x-a)^{k}.$$

My suggestion may look rather naive to the people here, but interpolation polynomials have long been used in computer algebra systems classically.