Can this function be rewritten to improve numerical stability?

We once had this (or at least something very similar) in my first numerical analysis class:

There it was proposed to calculate $f$ in terms of $y = \exp(-ux)$ instead. I.e. first set $y$ and then calculate $$f(x) = \frac{1-e^{-ux}}{u} = -x\frac{1-y}{\log(y)}$$

I have never understood why this should have any effect on the calculation, however... (It may only have an effect, if we make some assumptions on how the values get calculated inside the computer).

Edit: J.M. gave a very nice link in the comments, where it is explained why the above works!

If you make some numerical tests, let me know what you get (I'd be interested whether this actually does what they claimed it would)! ;)

The reason you get error is due to the machine precision of floating point arithmetic.

What is happening is not so much that $1-e^{ux}$ has issues with machine precision, but rather that you are dividing by something that is quite small, as well.

Essentially, if $e^{ux}$ is less than approximately $2.2 \times 10^{-16}$ (the machine $\varepsilon$) for double-precision arithmetic, then the total quantity $1-e^{ux}$ will return exactly 1. However, this is not so much an issue... 1 times something is very nearly the same thing as 0.99999999999999... times something.

The problem is then dividing by a very small number: $1/u$. If your $u < 2.2\times 10^{-16}$, then $1/u$ is very large, and the numerator is unable to "cancel" this effect. Furthermore, and I don't recall what the OFL limit for doubles are off the top of my head, if $1/u$ is larger than something like $2^{32}$, you're out of luck anyways.

One thing that you can do is compute $\log f(x)$ first: $\log f(x) = \log(1-e^{ux})-\log(u)$. You will still have the machine epsilon issue (ie, small values of $u$ will essentially make the log return $0$, but this is going to be nearly true for values in a neighborhood anyways), but your denominator will be more manageable.

Edit: Of course, one other method is to use the trusty ol' Taylor expansion for $e^{ut}$. Then, you have $f(x) = \frac{1}{u}+\sum_i \frac{(ux)^i}{i!u}$. The numerator is small, and the denominator is small, so they will to some degree cancel each other out. Compute each term in the sum independently, store them, and then add.

As J.M. pointed out in a comment to another answer, some platforms and programming languages have a primitive that computes $2^x-1$ or $e^x-1$ in a single step with without precision loss for small $x$. For example, in C99 <math.h> there is the expm1() function, which ought to map to the F2XM1 instruction on the x86 platform (via a scaling of $x$ by $\frac{1}{\log2}$).