DPQ (version 0.5-9)

pow1p: 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))

Value

numeric or number-like, as x + y.

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

Author

Originally proposed by Morten Welinder, see tools:::Rd_expr_PR(18642); tweaked, notably for small integer y, by Martin Maechler.

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.

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