DPQ (version 0.5-3)

log1pmx: Accurate log(1+x) - x Computation

Description

Compute $$\log(1+x) - x$$ accurately also for small \(x\), i.e., \(|x| \ll 1\).

Since April 2021, the pure R code version log1pmx() also works for "mpfr" numbers (from package Rmpfr).

Usage

log1pmx (x, tol_logcf = 1e-14, eps2 = 0.01, minL1 = -0.79149064,
         trace.lcf = FALSE,
         logCF = if(is.numeric(x)) logcf else logcfR.)
log1pmxC(x)  # TODO in future: arguments (minL1, eps2, tol_logcf),
             # possibly with *different* defaults (!)

Value

a numeric vector (with the same attributes as x).

Arguments

x

numeric (or "mpfr" number) vector with values \(x > -1\).

tol_logcf

a non-negative number indicating the tolerance (maximal relative error) for the auxiliary logcf() function.

eps2

non-negative cutoff where the algorithm switches from a few terms, to using logcf() explicitly. Note that for more accurate mpfr-numbers the default eps = .01 is too large, even more so when the tolerance is lowered (from 1e-14).

minL1

negative cutoff, called minLog1Value in Morten Welinder's C code for log1pmx() in R/src/nmath/pgamma.c, hard coded there to -0.79149064 which seems not optimal for computation of log1pmx(), at least in some cases, and hence the default may be changed in the future.

trace.lcf

logical used in logcf(.., trace=trace.lcf).

logCF

the function to be used as logcf(). The default chooses the pure R logcfR() when x is not numeric, and chooses the C-based logcf() when is.numeric(x) is true.

Author

A translation of Morten Welinder's C code of Jan 2005, see R's bug issue tools:::Rd_expr_PR(7307), parametrized and tuned by Martin Maechler.

Details

In order to provide full accuracy, the computations happens differently in three regions for \(x\), $$m_l = \code{minL1} = -0.79149064$$ is the first cutpoint,

\(x < m_l\) or \(x > 1\):

use log1pmx(x) := log1p(x) - x,

\(|x| < \epsilon_2\):

use \(t((((2/9 * y + 2/7)y + 2/5)y + 2/3)y - x)\),

\(x \in [ml,1]\), and \(|x| \ge \epsilon_2\):

use \(t(2y logcf(y, 3, 2) - x)\),

where \(t := \frac{x}{2 + x}\), and \(y := t^2\).

Note that the formulas based on \(t\) are based on the (fast converging) formula $$\log(1+x) = 2\left(r + \frac{r^3}{3}+ \frac{r^5}{5} + \ldots\right),$$ where \(r := x/(x+2)\), see the reference.

log1pmxC() is an interface to R C API (Rmathlib) function.

References

Abramowitz, M. and Stegun, I. A. (1972) Handbook of Mathematical Functions. New York: Dover. https://en.wikipedia.org/wiki/Abramowitz_and_Stegun provides links to the full text which is in public domain.
Formula (4.1.29), p.68.

Martin Mächler (2021). log1pmx, ... Computing ... Probabilities in R. (DPQ package vignette)

See Also

logcf, the auxiliary function, lgamma1p which calls log1pmx, log1p

Examples

(doExtras <- DPQ:::doExtras()) # TRUE e.g. if interactive()
n1 <- if(doExtras) 1001 else 201
curve(log1pmx, -.9999, 7, n=n1); abline(h=0, v=-1:0, lty=3)
curve(log1pmx, -.1,  .1,  n=n1); abline(h=0, v=0, lty=3)
curve(log1pmx, -.01, .01, n=n1) -> l1xz2; abline(h=0, v=0, lty=3)
## C and R versions correspond closely:
with(l1xz2, stopifnot(all.equal(y, log1pmxC(x), tol = 1e-15)))

e <- if(doExtras) 2^-12 else 2^-8; by.p <- 1/(if(doExtras) 256 else 64)
xd <- c(seq(-1+e, 0+100*e, by=e), seq(by.p, 5, by=by.p)) # length 676 or 5476 if do.X.
plot(xd, log1pmx(xd), type="l", col=2, main = "log1pmx(x)")
abline(h=0, v=-1:0, lty=3)

## much more graphics etc in ../tests/dnbinom-tst.R  (and the vignette, see above)