How to calculate the Factorial of Numbers > 170

[DataMiner]

Hi

I use factorial() function, but it doesn't work for numbers more than 170. I also use lfactorial() but doesn't work too.

I put simple codes below :

> factorial(200)
[1] Inf
Warning message:
In factorial(200) : value out of range in 'gammafn'

> exp(lfactorial(200))
[1] Inf

> gamma(201)
[1] Inf
Warning message:
value out of range in 'gammafn'

[Dason]

What are you trying to do? Do you need that exact value or are you using it for further calculations? For example: if you're using binomial coefficients these have factorials in the definition but you can simplify the calculations so you don't need to compute the full factorials (otherwise it would be impossible to compute binomial probabilities for large sample sizes).

[TheEcologist]

The factorial of x, x! is given by the gamma function +1 ; gamma (x+1).

Therefore use ?gamma in r.

factorial of 100;

gamma(100+1)

>170 is likely too large for the standard gamma function so use the log gamma function (for example the factorial of 200):

lgamma(200+1) is 863.232, so the factorial of 200 is e^863.232. Does this help?

[TheEcologist]

Otherwise you can use this (from Martin Maechler´s reply);


library(Rmpfr)
gamma(as(1000,"mpfr"))
#1 'mpfr' number of precision 128 bits
# 4.023872600770937735437024339230039857186 e 2564

[Dason]

I didn't know about

library(Rmpfr)
gamma(as(1000,"mpfr"))

I was considering writing a little bit of code to do arbitrary length integer calculations. Then I thought it had probably been done before so I just stopped and since it wasn't pressing to any of my research I just forgot about it.

But I still say the question remains as to what the numbers are being used for because a lot of times you can reduce the calculations down a lot if lots of things are going to cancel.

[TheEcologist]

In the case of binomial probabilities much will indeed 'cancel out' and such algorithms already exist :Catherine Loader (2000). Fast and Accurate Computation of Binomial Probabilities. And therefore have already been implemented in functions as dbinom - I just automatically assumed DataMiner needs this for something that isn't already implemented... but you are correct.. it would be way better if DataMiner "Describes the goal and not only the step".

[Outlier]

Isn't there a "Stirling's approximation" to calculate large factorials?

[Dason]

Yes. And it's quite good for large n in a relative error sense. So depending on what you're doing it could be a decent approximation. But even using the approximation... 170! is too large for R to handle using the base numeric types it has built in so you'd either need an arbitrary length integer or you need to specify what you're doing and figure out a better way to do it than to do the full computations. So it all really depends on if the OP ever shows up again.

[Lazar]

Stirling<- function (n)((2*pi*n)^.5)*(n/exp(1))^n

I think this is the Stirling approximation.

[TheEcologist]

Dason, I love mpfr. It has proven very useful in the past ! (when correct rounding becomes important)

For those of you who dont know about it:
http://en.wikipedia.org/wiki/Arbitra...ion_arithmetic

and http://www.mpfr.org/

[Dason]

Oh I didn't realize it was actually a C library. That might come in even more useful than I thought before...

[Dragan]

This approximation is very accurate. For n=200 the absolute relative error is less than 2.107x10^(-18)

[TheEcologist]

This approximation is very accurate. For n=200 the absolute relative error is less than 2.107x10^(-18)

Relative is the correct word :-)

[Dason]

Exactly. Depending on what you're doing you might be more concerned that the absolute error gets QUITE large.

[Outlier]

While we're on the subject, how many significant figures should be carried with stats work? Is it different for p < .05 than for p < .001?