0

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.

Jack M
  • 4,769
  • 6
  • 43
  • 67
  • Would using the [`decimal`](https://docs.python.org/3/library/decimal.html) module work for you? – Julia Mar 11 '21 at 20:31
  • 2
    Check if there is a function `expm1()`, which computes `exp(x)-1` in a numerically sound way, that is, without catastrophic cancellation near the origin. – njuffa Mar 11 '21 at 20:34
  • 1
    You could also make the same cancellation error twice and compute `(exp(t)-1)/((1+t)-1)`. // And yes, `numpy` has `expm1`. – Lutz Lehmann Mar 11 '21 at 20:50
  • @njuffa This seems to help partially (there is such a function in numpy), however, it seems in my case `t` is sometimes so small that `1 / t` sometimes returns `inf`, so I actually need a way of computing the expression "all at once", it seems. – Jack M Mar 11 '21 at 20:50
  • Then you need to replace the function with its Taylor-expansion if for instance `1+t*t==1`. If you do something a little more advanced like `(exp(t)-1-t)/(t**2)`, all the above tricks stop to work well. – Lutz Lehmann Mar 11 '21 at 20:53
  • @JackM `expm1(x)/x` works *perfectly* for tiny `x` when I try it in C, delivering a result of 1. Since 0/0 is undefined, it evaluates to NaN for x=0. `int main (void) { double x = 0x1.0p-1023; for (int i = 0; i < 53; i++) { printf ("x=%23.16e expm1(x)/x=%23.16e\n", x, expm1(x)/x); x = x * 0.5;} return 0; }` – njuffa Mar 11 '21 at 21:01
  • @njuffa I'm not sure how small that floating point is. For me I get `inf` when t is `3.16341187339184e-310`. – Jack M Mar 11 '21 at 21:32
  • @JackM Are you actually computing `expm1(x)/x`, or something else? For example, in floating-point arithmetic, `expm1(x)/x != expm1(x)*(1.0/x)`. In the latter, `(1.0/x)` can cause premature overflow. Maybe you can show your actual computation? I went by what the question in its present form is asking. `0x1.0p-1023` is about `1.1125e-308`, and the program halves it until `x` becomes zero. – njuffa Mar 11 '21 at 22:54
  • @njuffa I see, yes I'm using `(1/x)*expm1(x)` (I mean it's in Python using numpy, so I can be sure what's going on under the hood, but that's what my code looks like). Replacing it with `expm1(x)/x` seems to make an important difference. I'm still having problems but I think they're coming from other issues. – Jack M Mar 12 '21 at 08:32

1 Answers1

2

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.