numerically stable log1pexp calculation
Let $0 &lt; \varepsilon \lll 1$ be the relative error bound of the floating-point system—$2^{-53}$ in IEEE 754 binary64 arithmetic.

First, the naive formula log1p(exp(x)) always gives a good approximation unless exp(x) overflows:  If exp(x) computes $(1 + \delta_1) e^x$, and if log1p(y) computes $(1 + \delta_2) \log(1 + y)$, where $|\delta_1|, |\delta_2| \leq \varepsilon$, then log1p(exp(x)) yields

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

where $|\delta_1'|$ is bounded by a constant multiple of $\delta_1$ since $(1 + \delta_1) e^x$ is positive and $y \mapsto \log(1 + y)$ has small condition number on the positive real axis.  Thus the relative error is bounded by $|\delta_2| + |\delta_1'| + |\delta_2 \delta_1'|$ which in turn is bounded by small polynomial function of $\varepsilon$.


  Exercise: Compute a bound on $|\delta_1'/\delta_1|$.


But you can do better, with less computation (only one transcendental evaluation, not two), in some special cases:


If $x \lll 0$, then $e^x$ is near 0 so $\log(1 + e^x) \approx e^x$.  Specifically, $0 &lt; e^x &lt; 1$, so:


$\log(1 + e^x) \geq e^x/2$, since $y \mapsto \log(1 + y)$ is a concave function with derivative $y \mapsto 1/(1 + y)$ bounded below by $1/2$ for $y &lt; 1$; and
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|}
     &amp;\leq 2 e^{-x} |e^x - \log(1 + e^x)| \\
     &amp;= 2 e^{-x} |e^x - (e^x - e^{2x}/2 + e^{3x}/3 - \cdots)| \\
     &amp;= 2 e^{-x} |e^{2x}/2 - e^{3x}/3 + e^{4x}/4 - \cdots| \\
     &amp;= 2 |e^x/2 - e^{2x}/3 + e^{3x}/4 - \cdots| \\
     &amp;&lt; 2 |e^x/2| \\
     &amp;= e^x.
   \end{align}

For $x \leq \log \varepsilon$, this relative error is bounded by $\varepsilon$, which justifies the formula exp(x).

In IEEE 754 binary64 arithmetic, $\log \varepsilon \approx -37$.
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 &lt; y &lt; 1$, with $y = e^{-x}$, as in case (1).  The value of $\log(1 + e^x)$ lies between $x + e^{-x}$ and $x + e^{-x} - e^{-2x}/2$, so the relative error of $x + e^{-x}$ from $\log(1 + e^x)$ is below the relative error of $x + e^{-x}$ from $x + e^{-x} - e^{-2x}/2$, which is

\begin{equation}
     \frac{|x + e^{-x} - (x + e^{-x} - e^{-2x}/2)|}{|x + e^{-x} - e^{-2x}/2|}
      = \frac{e^{-2x}/2}{|x + e^{-x} - e^{-2x}/2|}
      \leq e^{-2x}/2.
   \end{equation}

For $x \geq \log(1/\sqrt{2\varepsilon})$, this is bounded by $\varepsilon$, which justifies the formula x + exp(-x).

In IEEE 754 binary64 arithmetic, $\log(1/\sqrt{2\varepsilon}) \approx 18$.
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)}
     &lt; \log(1 + e^{-x}),
   \end{equation}

as long as $e^x \geq e - 1$ so that the denominator $\log(1 + e^x)$ is at least $1$.

For $x \geq \log[1/(e^\varepsilon - 1)]$, then this error is below $\varepsilon$, which justifies the formula x.

In IEEE 754 binary64 arithmetic, $\log[1/(e^\varepsilon - 1)] \approx 37$.


  Exercise: Find a nice formula for a bound on $x$ that corresponds to the cutoff $33.3$ in the source you cited, by taking the denominator $\log(1 + e^x)$ into account.