Potenza del numero vicino a 1

Question

Immagino che ci sia qualche trucco standard che non ero in grado di trovare: comunque voglio calcolare una grande potenza di un numero molto vicino a 1 (pensa 1-p dove P <1E-17) in aModa numericamente stabile.1-P viene troncato a 1 sul mio sistema.

Usando l'espansione Taylor del logaritmo ottengo i seguenti limiti

C'è qualcosa di più intelligente che posso fare?

Answer 1

È possibile calcolare il log(1+x) in modo più accurato per |x| <= 1 utilizzando la funzione log1p.

Un esempio:

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

.

e un altro:

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

.

Qui, tuttavia, siamo limitati con l'accuratezza dell'aritmetica FP di tipo double a base di tipo (vedere Cosa dovrebbe sapere ogni scienziato informatico sul punto flottante aritmetico ).

Un altro modo è quello di usare ad es.Il Rmpfr (a.k.a. Punto flottante di precisione multiplo Affidabile) Pacchetto:

> 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

.

Il pacchetto utilizza la libreria gsl e mpfr per implementare operazioni arbitrarie di precisione FP (al costo della velocità di calcolo più lenta, ovviamente).