Puissance du nombre proche de 1

Question

Je suppose qu'il existe une astuce standard que je n'ai pas pu trouver :Quoi qu'il en soit, je veux calculer une grande puissance d'un nombre très proche de 1 (pensez à 1-p où p<1e-17) de manière numériquement stable.1-p est tronqué à 1 sur mon système.

En utilisant le développement de Taylor du logarithme, j'obtiens les limites suivantes

Y a-t-il quelque chose de plus intelligent que je puisse faire ?

Answer 1

Vous pouvez calculer log(1+x) plus précisément pour |x| <= 1 en utilisant le log1p fonction.

Un exemple:

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

Et un autre:

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

Ici, cependant, nous sommes limités par la précision de double arithmétique FP basée sur le type (voir Ce que tout informaticien devrait savoir sur l'arithmétique à virgule flottante).

Une autre façon consiste à utiliser par ex.le Rmpfr (alias.Package fiable à virgule flottante à précision multiple) :

> 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

Le paquet utilise le gsl et mpfr bibliothèque pour implémenter des opérations FP de précision arbitraire (au prix d'une vitesse de calcul plus lente, bien sûr).