pow1p: Accurate (1+x)^y, notably for small |x|

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

Accurate (1+x)^y, notably for small |x|

Description

Compute (1+x)^y accurately, notably also for small |x|, where the naive formula suffers from cancellation, returning 1, often.

Usage

pow1p(x, y,
      pow = ((x + 1) - 1) == x || abs(x) > 0.5 || is.na(x))

Arguments

x, y	numeric or number-like; in the latter case, arithmetic incl. ^, comparison, exp, log1p, abs, and is.na methods must work.
pow	logical indicating if the “naive” / direct computation (1 + x)^y should be used (unless y is in 0:4, where the binomial is used, see ‘Details’). The current default is the one used in R's C-level function (but beware of compiler optimization there!).

Details

A pure R-implementation of R 4.4.0's new C-level pow1p() function which was introduced for more accurate dbinom_raw() computations.

Currently, we use the “exact” (nested) polynomial formula for y \in \{0,1,2,3,4\}.

MM is conjecturing that the default pow=FALSE for (most) x \le \frac 1 2 is sub-optimal.

Value

numeric or number-like, as x + y.

Author(s)

Originally proposed by Morten Welinder, see \Sexpr[results=rd]{tools:::Rd_expr_PR(18642)}; tweaked, notably for small integer y, by Martin Maechler.

See Also

^, log1p, dbinom_raw.

Examples

x <- 2^-(1:50)
y <- 99
f1 <- (1+x)^99
f2 <- exp(y * log1p(x))
fp <- pow1p(x, 99)
matplot(x, cbind(f1, f2, fp), type = "l", col = 2:4)
legend("top", legend = expression((1+x)^99, exp(99 * log1p(x)), pow1p(x, 99)),
       bty="n", col=2:4, lwd=2)
cbind(x, f1, f2, sfsmisc::relErrV(f2, f1))

DPQ documentation built on Sept. 11, 2024, 8:37 p.m.