We all know that
\begin{equation}
\exp(x) = \sum_{n=0}^\infty \frac{x^n}{n!} = 1 + x + \frac12 x^2 + \dots
\end{equation}
implies that for $|x| \ll 1$, we have $\exp(x) \approx 1 + x$. This means that if we have to evaluate in floating point $\exp(x) -1$, for $|x| \ll 1$ catastrophic cancellation can occur.
This can be easily demonstrated in python:
>>> from math import (exp, expm1)
>>> x = 1e-8
>>> exp(x) - 1
9.99999993922529e-09
>>> expm1(x)
1.0000000050000001e-08
>>> x = 1e-22
>>> exp(x) - 1
0.0
>>> expm1(x)
1e-22
Exact values are
\begin{align}
\exp(10^{-8}) -1 &= 0.000000010000000050000000166666667083333334166666668 \dots \\
\exp(10^{-22})-1 &= 0.000000000000000000000100000000000000000000005000000 \dots
\end{align}
In general an "accurate" implementation of exp
and expm1
should be correct to no more than 1ULP (i.e. one unit of the last place). However, since attaining this accuracy results in "slow" code, sometimes a fast, less accurate implementation is available. For example in CUDA we have expf
and expm1f
, where f
stands for fast. According to the CUDA C programming guide, app. D the expf
has an error of 2ULP.
If you do not care about errors in the order of few ULPS, usually different implementations of the exponential function are equivalent, but beware that bugs may be hidden somewhere... (Remember the Pentium FDIV bug?)
So it is pretty clear that expm1
should be used to compute $\exp(x)-1$ for small $x$. Using it for general $x$ is not harmful, since expm1
is expected to be accurate over its full range:
>>> exp(200)-1 == exp(200) == expm1(200)
True
(In the above example $1$ is well below 1ULP of $\exp(200)$, so all three expression return exactly the same floating point number.)
A similar discussion holds for the inverse functions log
and log1p
since $\log(1+x) \approx x$ for $|x| \ll 1$.
log1p
you're referring to (especially how it's implemented, so we don't have to guess). $\endgroup$