numero - funciones r en r ejemplos

En R, ¿por qué el factorial(100) se muestra de manera diferente a prod(1: 100)?

Question

En RI estoy encontrando un comportamiento extraño que no puedo explicar y espero que alguien aquí pueda. Creo que el valor de 100! es este gran numero

Unas pocas líneas de la consola que muestran el comportamiento esperado ...

>factorial( 10 )
[1] 3628800
>prod( 1:10 )
[1] 3628800
> prod( as.double(1:10) )
[1] 3628800
> cumprod( 1:10 )
[1]       1       2       6      24     120     720    5040   40320  362880 3628800

Sin embargo, cuando intento 100! Obtengo (observe cómo los números resultantes comienzan a diferir en aproximadamente 14 dígitos):

> options(scipen=200) #set so the whole number shows in the output
> factorial(100)
[1] 93326215443942248650123855988187884417589065162466533279019703073787172439798159584162769794613566466294295348586598751018383869128892469242002299597101203456
> prod(1:100)
[1] 93326215443944102188325606108575267240944254854960571509166910400407995064242937148632694030450512898042989296944474898258737204311236641477561877016501813248
> prod( as.double(1:100) )
[1] 93326215443944150965646704795953882578400970373184098831012889540582227238570431295066113089288327277825849664006524270554535976289719382852181865895959724032
> all.equal( prod(1:100) , factorial(100) , prod( as.double(1:100) ) )
[1] TRUE

Si hago algunas pruebas contra una variable establecida en el número ''conocido'' de 100! Entonces veo lo siguiente:

# This is (as far as I know) the ''true'' value of 100!
> n<- as.double(93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000)
> factorial(100) - n
[1] -1902315522848807765998160811905210717565551993186466795054798772271710903343294674760811531554315419925519536152107160826913610179566298858520576
> prod(1:100) - n
[1] -48777321098687378615337456715518223527321845979140174232174327494146433419058837814379782860367062049372295798771978482741374619988879457910784
> prod(as.double(1:100)) - n
[1] 0

El resultado final se evalúa a cero, pero el número devuelto para prod( as.double( 1:100 ) ) no se muestra como esperaría, a pesar de que evalúa correctamente prod( as.double( 1:100 ) ) - n donde n es una variable establecida en el valor de 100 !.

¿Alguien puede explicarme este comportamiento por favor? No debería estar relacionado con el desbordamiento, etc., que yo sepa, ya que estoy usando un sistema x64. Versión y información de la máquina a continuación:

> .Machine$double.xmax
[1] 1.798e+308
> str( R.Version() )
List of 14
 $ platform      : chr "x86_64-apple-darwin9.8.0"
 $ arch          : chr "x86_64"
 $ os            : chr "darwin9.8.0"
 $ system        : chr "x86_64, darwin9.8.0"
 $ status        : chr ""
 $ major         : chr "2"
 $ minor         : chr "15.2"
 $ year          : chr "2012"
 $ month         : chr "10"
 $ day           : chr "26"
 $ svn rev       : chr "61015"
 $ language      : chr "R"
 $ version.string: chr "R version 2.15.2 (2012-10-26)"
 $ nickname      : chr "Trick or Treat"

¿Puede alguien explicarme esto? No dudo que R haga todo correctamente y esto es muy probable que esté relacionado con el uso. Podría indicar que ya que prod( as.double( 1:100 ) ) - n evalúa correctamente lo que me molesta, pero estoy haciendo Problema de Project Euler 20, por lo que necesitaba que se mostraran los dígitos correctos.

Gracias

Answer 1

Agregaré una tercera respuesta solo para describir gráficamente el comportamiento que está encontrando. Esencialmente, la precisión doble para el cálculo factorial es suficiente hasta 22, luego comienza a separarse cada vez más del valor real.

¡Alrededor de los 50 !, hay una distinción adicional entre los dos métodos factorial (x) y prod (1: x), con este último que produce, como usted indicó, valores más similares al factor "real".

Código adjunto:

# Precision of factorial calculation (very important for the Fisher''s Exact Test)
library(gmp)
perfectprecision<-list()
singleprecision<-c()
doubleprecision<-c()
for (x in 1:100){
    perfectprecision[x][[1]]<-factorialZ(x)
    singleprecision<-c(singleprecision,factorial(x))
    doubleprecision<-c(doubleprecision,prod(1:x))
}


plot(0,col="white",xlim=c(1,100),ylim=c(0,log10(abs(doubleprecision[100]-singleprecision[100])+1)),
        ,ylab="Log10 Absolute Difference from Big Integer",xlab="x!")
for(x in 1:100) {
    points(x,log10(abs(perfectprecision[x][[1]]-singleprecision[x])+1),pch=16,col="blue")
    points(x,log10(abs(perfectprecision[x][[1]]-doubleprecision[x])+1),pch=20,col="red")
}
legend("topleft",col=c("blue","red"),legend=c("factorial(x)","prod(1:x)"),pch=c(16,20))

Answer 2

Bien, puedes decir del cuerpo de factorial que llama gamma , que llama .Primitive("gamma") . ¿Qué .Primitive("gamma") tiene .Primitive("gamma") ? Como este

Para entradas grandes, el .Primitive("gamma") está en la línea 198 de ese código. Esta llamando

exp((y - 0.5) * log(y) - y + M_LN_SQRT_2PI +
            ((2*y == (int)2*y)? stirlerr(y) : lgammacor(y)));

que es sólo una aproximación .

Por cierto, el artículo sobre Rmpfr usa factorial como ejemplo. Entonces, si estás tratando de resolver el problema, "solo usa la biblioteca Rmpfr ".

Answer 3

Esto tiene que ver no con el valor máximo para un double sino con su precisión.

100! Tiene 158 dígitos (decimales) significativos. IEEE double s (64 bits) tiene 52 bits de espacio de almacenamiento para la mantisa, por lo que se obtienen errores de redondeo después de que se han excedido unos 16 dígitos decimales de precisión.

Por cierto, 100! es de hecho, como sospechabas,

93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

por lo que todos los valores R calculados son incorrectos.

Ahora no sé R, pero parece que all.equal() convierte los tres de esos valores en float antes de compararlos, por lo que sus diferencias se pierden.

Answer 4

Tu prueba con all.equal no produce lo que esperas. all.equal solo puede comparar dos valores. El tercer argumento se hace coincidir en la posición con la tolerance , lo que da la tolerancia de la operación de comparación. ¡En tu invocación a all.equal , dale una tolerancia de 100! lo que definitivamente lleva a que la comparación sea verdadera para valores absurdamente diferentes:

> all.equal( 0, 1000000000, prod(as.double(1:100)) )
[1] TRUE

Pero incluso si le das dos argumentos solamente, por ejemplo

all.equal( prod(1:100), factorial(100) )

todavía produciría TRUE porque la tolerancia predeterminada es .Machine$double.eps ^ 0.5 , por ejemplo, los dos operandos tienen que coincidir con aproximadamente 8 dígitos, lo que definitivamente es el caso. Por otro lado, si establece la tolerancia en 0 , entonces ninguna de las tres combinaciones posibles emerge igual a la comparación:

> all.equal( prod(1:100), factorial(100), tolerance=0.0 )
[1] "Mean relative difference: 1.986085e-14"
> all.equal( prod(1:100), prod( as.double(1:100) ), tolerance=0.0 )
[1] "Mean relative difference: 5.22654e-16"
> all.equal( prod(as.double(1:100)), factorial(100), tolerance=0.0 )
[1] "Mean relative difference: 2.038351e-14"

También tenga en cuenta que solo porque le ha dicho a R que imprima 200 números significativos no significa que todos sean correctos. De hecho, 1/2 ^ 53 tiene aproximadamente 53 dígitos decimales, pero solo los primeros 16 se consideran significativos.

Esto también hace que su comparación con el valor "verdadero" sea defectuosa. Observe esto. Los dígitos finales en lo que R te da para factorial(100) son:

...01203456

¡Resta n de él, donde n es el valor "verdadero" de 100! por lo tanto, debería tener 24 ceros al final y, por lo tanto, la diferencia también debería terminar con los mismos dígitos que el factorial(100) . Pero más bien termina con:

...58520576

Esto solo muestra que todos esos dígitos no son significativos y uno no debería realmente mirar su valor.

¡Se necesitan 525 bits de precisión binaria para representar exactamente 100! - Eso es 10x la precisión del double .