numerically stable log1pexp calculation

What are good approximations for computing log1pexp for single precision and double precision floating point numbers?

Note: log1pexp(x) is log(1 + exp(x))

I have found few implementations of log1pexp for double precision but they don't provide an explanation on how they arrived at the approximations. Hence, I am not able to implement log1exp for single precision numbers (without converting to double precision intermediates of course).

Reference implementation: https://github.com/davisking/dlib/blob/master/dlib/dnn/utilities.h#L16-L29

Let $0

• the Taylor expansion of $y \mapsto \log(1 + y)$ converges absolutely at $e^x$.

Thus the relative error of $e^x$ from $\log(1 + e^x)$ is:

\begin{align}
\frac{|e^x - \log(1 + e^x)|}{\left|\log(1 + e^x)\right|}
&\leq 2 e^{-x} |e^x - \log(1 + e^x)| \\
&= 2 e^{-x} |e^x - (e^x - e^{2x}/2 + e^{3x}/3 - \cdots)| \\
&= 2 e^{-x} |e^{2x}/2 - e^{3x}/3 + e^{4x}/4 - \cdots| \\
&= 2 |e^x/2 - e^{2x}/3 + e^{3x}/4 - \cdots| \\
&

-
If $x \gg 0$, we can use the identity

\begin{equation}
\log(1 + e^x) = \log\bigl[e^x (1 + e^{-x})\bigr] = x + \log(1 + e^{-x})
\end{equation}

to rewrite $\log(1 + e^x)$ in terms of $\log(1 + y)$ for $0

-
If $x \ggg 0$, then $x + \log(1 + e^{-x})$ may simply be rounded to $x$. Specifically, the relative error of $x$ from $\log(1 + e^x) = x + \log(1 + e^{-x})$ is

\begin{equation}
\frac{\bigl|x - \bigl[x + \log(1 + e^{-x})\bigr]\bigr|}{\left|\log(1 + e^x)\right|}
= \frac{\log(1 + e^{-x})}{\log(1 + e^x)}

StackExchange Computer Science Q#110798, answer score: 6