How do I evaluate xe^x/(e^x-1) with numerical stability in Python?
<p>How do I evaluate xe^x/(e^x-1) with numerical stability around zero and when x is very positive or negative?  I have access to all the usual mathematical functions in <code>numpy</code> and <code>scipy</code>.</p>

<pre><code>def f(x):
  if abs(x) &gt; 0.1: return x*exp(x)/(exp(x)-1)
  else: return 1/(1.-x/2.+x**2/6.-x**3/24.)
</code></pre>

<p>The expansion in the last line can be extended in the obvious fashion if more precision is required, and can be made faster by re-phrasing.  As it stands, it errs by as much as 1e-6.</p>

python numerical-stability