Search code examples
rfloating-pointfloating-point-precisionnumerical-stability

Power of number close to 1

I'm guessing there is some standard trick that I wasn't able to find: Anyway I want to compute a large power of a number very close to 1(think 1-p where p<1e-17) in a numerically stable fashion. 1-p is truncated to 1 on my system.

Using the taylor expansion of the logarithm I obtain the following bounds

formula

Is there anything smarter I can do?

Solution

  • You may calculate log(1+x) more accurately for |x| <= 1 by using the log1p function.

    An example:

    > p <- 1e-17
> log(1-p)
[1] 0
> log1p(-p)
[1] -1e-17

    And another one:

    > print((1+1e-17)^100, digits=22)
[1] 1
> print(exp(100*log1p(-1e-17)), digits=22)
[1] 0.9999999999999990007993

    Here, however, we're limited with the accuracy of double type-based FP arithmetic (see What Every Computer Scientist Should Know About Floating-Point Arithmetic).

    Another way is to use e.g. the Rmpfr (a.k.a. Multiple Precision Floating-Point Reliable) package:

    > options(digits=22)
> library(Rmpfr)
> .N <- function(.) mpfr(., precBits = 200) # see the package's vignette
> (1-.N(1e-20))^100
1 'mpfr' number of precision  200   bits 
[1] 0.99999999999999999900000000000000005534172854579042829381053529

    The package uses the gsl and mpfr library to implement arbitrary precision FP operations (at the cost of slower computation speed, of course).