rexpm1 {DPQ}

TOMS 708 Approximation REXP(x) of expm1(x) = exp(x) - 1

Description

Originally REXP(), now rexpm1() is a numeric (double precision) approximation of exp(x) - 1, notably for small |x| \ll 1 where direct evaluation looses accuracy through cancellation.

Fully accurate computations of exp(x) - 1 are now known as expm1(x) and have been provided by math libraries (for C, C++, ..) and R, (and are typically more accurate than rexp1()).

The rexpm1() approximation was developed by Didonato & Morris (1986) and uses a minimax rational approximation for |x| <= 0.15; the authors say “accurate to within 2 units of the 14th significant digit” (top of p.379).

Usage

rexpm1(x)

Arguments

x

a numeric vector.

Value

a numeric vector (or array) as x,

Author(s)

Martin Maechler, for the C to R *vectorized* translation.

References

Didonato, A.R. and Morris, A.H. (1986) Computation of the Incomplete Gamma Function Ratios and their Inverse. ACM Trans. on Math. Softw. 12, 377–393, doi:10.1145/22721.23109; The above is the “flesh” of ‘TOMS 654’:

Didonato, A.R. and Morris, A.H. (1987) Algorithm 654: FORTRAN subroutines for Compute the Incomplete Gamma Function Ratios and their Inverse. ACM Transactions on Mathematical Software 13, 318–319, doi:10.1145/29380.214348.

See Also

pbeta, where the C version of rexpm1() has been used in several places, notably in the original TOMS 708 algorithm.

Examples

x <- seq(-3/4, 3/4, by=1/1024)
plot(x,     rexpm1(x)/expm1(x) - 1, type="l", main = "Error wrt expm1()")
abline(h = (-8:8)*2^-53, lty=1:2, col=adjustcolor("gray", 1/2))
cb2 <- adjustcolor("blue", 1/2)
do.15 <- function(col = cb2) {
    abline(v = 0.15*(-1:1), lty=3, lwd=c(3,1,3), col=col)
    axis(1, at=c(-.15, .15), col=cb2, col.axis=cb2)
}
do.15()

 op <- par(mar = par("mar") + c(0,0,0,2))
plot(x, abs(rexpm1(x)/expm1(x) - 1),type="l", log = 'y',
     main = "*Relative* Error wrt expm1() [log scale]")#, yaxt="n"
abline(h = (1:9)*2^-53, lty=2, col=adjustcolor("gray", 1/2))
axis(4, at = (1:9)*2^-53, las = 1, labels =
     expression(2^-53, 2^-52, 3 %*% 2^-53, 2^-51, 5 %*% 2^-53,
                6 %*% 2^-53, 7 %*% 2^-53, 2^-50, 9 %*% 2^-53))
do.15()
 par(op)

## "True" Accuracy comparison of  rexpm1() with [OS mathlib based] expm1():
if(require("Rmpfr")) withAutoprint({
  xM <- mpfr(x, 128); Xexpm1 <- expm1(xM)
  REr1 <- asNumeric(rexpm1(x)/Xexpm1 - 1)
  REe1 <- asNumeric(expm1(x) /Xexpm1 - 1)
  absC <- function(E) pmax(2^-55, abs(E))

  plot(x, absC(REr1), type= "l", log="y",
       main = "|rel.Error|  of exp(x)-1 computations wrt 128-bit MPFR ")
  lines(x, absC(REe1), col = (c2 <- adjustcolor(2, 3/4)))
  abline(h = (1:9)*2^-53, lty=2, col=adjustcolor("gray60", 1/2))
  do.15()
  axis(4, mgp=c(2,1/4,0),tcl=-1/8, at=2^-(53:51), labels=expression(2^-53, 2^-52, 2^-51), las=1)
  legend("topleft", c("rexpm1(x)", " expm1(x)"), lwd=2, col=c("black", c2),
         bg = "gray90", box.lwd=.1)

})