В R почему факториал (100) отображается по-разному с prod (1:100)?

Ответ 1

Ваш тест с помощью all.equal не дает ожидаемого результата. all.equal может сравнивать только два. Третий аргумент позиционируется на tolerance, что дает допуск операции сравнения. В вашем вызове all.equal вы даете ему допуск 100!, который определенно приводит к тому, что сравнение истинно для абсурдно разных значений:

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

Но даже если вы укажете только два аргумента, например

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

он все равно произведет TRUE, потому что допустимое значение по умолчанию .Machine$double.eps ^ 0.5, например. два операнда должны соответствовать примерно 8 цифрам, что определенно имеет место. С другой стороны, если вы установите допуск 0, то ни одна из трех возможных комбинаций не станет равной сравнению:

> 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"

Также обратите внимание, что только потому, что вы сказали R напечатать 200 значительных чисел, это не значит, что они все правильные. Действительно, 1/2 ^ 53 имеет около 53 десятичных цифр, но только первые 16 считаются значимыми.

Это также делает ваше сравнение с истинным значением ошибочным. Соблюдайте это. Конечными цифрами в том, что R дает вам для factorial(100), являются:

...01203456

Вы вычитаете из него n, где n - это "истинное" значение 100! поэтому он должен иметь 24 нуля в конце и, следовательно, разница также должна заканчиваться теми же цифрами, что и factorial(100). Но скорее это заканчивается:

...58520576

Это только показывает, что все эти цифры несущественны, и не следует действительно смотреть на их значение.

Требуется 525 бит бинарной точности, чтобы точно представлять 100! - что 10x точность double.

Ответ 2

Это должно быть сделано не с максимальным значением для double, а с его точностью.

100! имеет 158 значащих (десятичных) цифр. IEEE double (64 бит) имеет 52 бит пространства для хранения мантиссы, поэтому вы получаете ошибки округления после превышения примерно 16 десятичных цифр точности.

Кстати, 100! на самом деле, как вы подозревали,

93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

поэтому все вычисленные значения R неверны.

Теперь я не знаю R, но кажется, что all.equal() преобразует все три этих значения в float перед сравнением, и поэтому их отличия теряются.

Ответ 3

Я добавлю третий ответ, чтобы графически описать поведение, с которым вы сталкиваетесь. По сути, двойная точность для факториального расчета достаточна до 22!, затем она начинает все больше расходиться с реальным значением.

Вокруг 50! существует еще одно различие между двумя методами factorial (x) и prod (1: x), причем последнее дает, как вы указали, значения, более похожие на "реальный" фактор.

Код прилагается:

# Precision of factorial calculation (very important for the Fisher 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))

Ответ 4

Ну, вы можете сказать из тела factorial, что он вызывает gamma, который вызывает .Primitive("gamma"). Что выглядит .Primitive("gamma")? Как это.

Для больших входов поведение .Primitive("gamma") находится на строке 198 этого кода. Он вызывает

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

который приближенно.

Кстати, статья в Rmpfr использует factorial в качестве примера. Поэтому, если вы пытаетесь решить проблему, просто используйте библиотеку Rmpfr.