R: Accurate 'log(1+x) - x' Computation

log1pmx {DPQ}

R Documentation

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).

rlog1(x), provided mostly for reference and reproducibility, is used in TOMS Algorithm 708, see e.g. the reference of lgamma1p. and computes minus log1pmx(x), i.e., x - \log(1+x), using (argument reduction) and a rational approximation when x \in [-0.39, 0.57).

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 (!)
rlog1(x)

Arguments

`x`	numeric (or, for `log1pmx()` only, `"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. Also, for mpfr numbers, the default -0.79149064 may well be far from optimal.
`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.

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.

Value

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

Author(s)

A translation of Morten Welinder's C code of Jan 2005, see R's bug issue PR#7307, parametrized and tuned by Martin Maechler.

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)

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)

## --- Compare rexp1() with log1pmx() ----------------------------
x <- seq(-0.5, 5/8, by=1/256)
all.equal(log1pmx(x), -rlog1(x), tol = 0) # 2.838e-16 {|rel.error| <= 1.33e-15}
stopifnot(all.equal(log1pmx(x), -rlog1(x), tol = 1e-14))
## much more closely:
x <- c(-1+1e-9, -1+1/256, -(127:50)/128, (-199:295)/512, 74:196/128)
if(is.unsorted(x)) stop("x must be sorted for plots")
rlog1.x <- rlog1(x)
summary(relD <- sfsmisc::relErrV(log1pmx(x), -rlog1.x))
n.relD <- relD * 2^53
table(n.relD)
## 64-bit Linux F36 (gcc 12.2.1):
## -6  -5  -4  -3  -2  -1   0   2   4   6   8  10  12  14
##  2   3  13  24  79  93 259 120  48  22  14  15   5   1
stopifnot(-10 <= n.relD, n.relD <= 20) # above Lnx: [-6, 14]

if(requireNamespace("Rmpfr")) {
  relE <- Rmpfr::asNumeric(sfsmisc::relErrV(log1pmx(Rmpfr::mpfr(x,128)), -rlog1(x)))
  plot(x, pmax(2^-54, abs(relE)), log="y", type="l", main= "|rel.Err| of rlog1(x)")
  rl1.c <- c(-.39, 0.57, -.18, .18) # the cutoffs used inside rlog1()
  lc <- "gray"
  abline(v = rl1.c, col=lc, lty=2)
  axis(3, at=rl1.c, col=lc, cex.axis=3/4, mgp=c(2,.5,0))
  abline(h= (1:4)*2^-53,  lty=3, col = (cg <- adjustcolor(1, 1/4)))
  axis(4, at=(1:4)*2^-53, labels=expression(frac(epsilon[c],2), epsilon[c],
                                            frac(3,2)*epsilon[c], 2*epsilon[c]),
       cex.axis = 3/4, tcl=-1/4, las = 1, mgp=c(1.5,.5,0), col=cg)
  ## it seems the -.18 +.18 cutoffs should be slightly moved "outside"
}

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

[Package DPQ version 0.5-8 Index]

Accurate log(1+x) - x Computation