How can I compute (exp(t) - 1)/t in a numerically stable way?

Question

The expression (exp(t) - 1)/t converges to 1 as t tends to 0. However, when computed numerically, we get a different story:

In [19]: (exp(10**(-12)) - 1) * (10**12)                                        
Out[19]: 1.000088900582341

In [20]: (exp(10**(-13)) - 1) * (10**13)                                        
Out[20]: 0.9992007221626409

In [21]: (exp(10**(-14)) - 1) * (10**14)                                        
Out[21]: 0.9992007221626409

In [22]: (exp(10**(-15)) - 1) * (10**15)                                        
Out[22]: 1.1102230246251565

In [23]: (exp(10**(-16)) - 1) * (10**16)                                        
Out[23]: 0.0

Is there some way I can compute this expression without encountering these problems? I've thought of using a power series but I'm wary of implementing this myself as I'm not sure of implementation details like how many terms to use.

If it's relevant, I'm using Python with scipy and numpy.

Answer 1

The discussion in the comments about tiny values seems pointless. If t is so tiny that it causes underflow, then the expression is 1 "since a long time". Indeed the Taylor development is

1 + t/2 + t²/6 + t³/24...

and as soon as t < 1 ulp, the floating-point representation is exactly 1.

Above that, expm1(t)/t will do a good job.