\SweaveOpts{engine=R, keep.source=TRUE} \SweaveOpts{eps=FALSE, pdf=TRUE, width=9, height=6, strip.white=true} \setkeys{Gin}{width=\textwidth} <>= options(SweaveHooks= list(fig=function() par(mar=c(5.1, 4.1, 1.1, 2.1))), width = 75, digits = 7, # <-- here, keep R's default! prompt = "> ", continue=" ") try(Mlibrary(Rmpfr)) stopifnot(require("Rmpfr")) @ \subsection{Why $\log(1+x)$ is not good enough, but $\mathrm{log1p}(x)$ is} \begin{frame}[fragile]{Why $\log(1+x)$ is not good enough, but $\mathrm{log1p}(x)$ is} $1+x$ cannot be numerically accurate when $|x| \ll 1$. In double precision (53 bits $\approx 16$ digits) accuracy, $1+x$ sees'' only 2--3 digits of $x$ when $x = 10^{-14}$, <>= u <- 1 + (e <- 4e-13/9) ## then u - 1 == e mathematically: rbind(u-1 = u - 1, e) @ And the consequence for $\log(1 + x)$, <>= curve(abs(1 - log(1+x) / log1p(x)), 1e-17, .2, log = 'xy', main = "..", ..) sfsmisc :: eaxis(1); eaxis(2) <>= curve(abs(1 - log(1+x)/log1p(x)), 1e-17, .2, log='xy', n=1024, ylab='', main="| relative error of log(1+x) |", lwd = .25, col="red", axes=FALSE, frame=TRUE) sfsmisc :: eaxis(1); eaxis(2) @ \end{frame} \begin{frame}{Why $\mathrm{log1p}(x)$ beats $\log(1+x)$} \begin{center} \includegraphics[width=\textwidth]{log1p-ex-log1p-err-do} \end{center} Solution: Expand $\log(1+x)$ around $x=0$. Well known $\log(1+x) = x - x^2/2 + x^3/3 \pm \ldots = \sum_{n=1}^\infty (-1)^{n+1} \frac{x^n}{n},$ for $|x| < 1$. \medskip Fast version of this expansion: typically used in % the implementation of \Rfun{log1p}. \end{frame} %%% Local Variables: %%% TeX-command-default: "LaTeX PDF" %%% TeX-master: "Maechler-2014.tex" %%% End: