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
Is there anything smarter I can do?
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).
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With