Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

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?

like image 897
Tobias Madsen Avatar asked Jun 03 '14 10:06

Tobias Madsen


1 Answers

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

like image 172
gagolews Avatar answered Sep 18 '22 22:09

gagolews