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

In DPQ: Density, Probability, Quantile ('DPQ') Computations

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 (double precision) 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 \Sexpr[results=rd]{tools:::Rd_expr_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)

See Also

logcf, the auxiliary function, lgamma1p which calls log1pmx, log1p; also expm1x)() which computes expm1(x) - x accurately, whereas log1pmx(x) computes log1p(x) - x accurately

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)

DPQ documentation built on Nov. 3, 2024, 3 a.m.