[R] large factorials

[Ben Bolker]

J Dougherty wrote:
> 
> If you need an open source arbitrary precision calculator, you might want
> to 
> look at Octave which is OS and works similarly to Mathematica  - up to a
> point 
> books for Mathematica will be a significant help with Octave.
> 
> JDougherty
> 
> 

  ???  octave is a Matlab clone, not a Mathematica clone (I would be
interested in an open source Mathematica clone ...) ???

  octave seems to have the same (double-precision floating point)
limitations
as R:

octave:13> factorial(170)
ans =  7.2574e+306
octave:14> factorial(171)
ans = Inf

Another way you can do this is with the "Brobdingnag" package,
which represents large numbers by their logarithms ... 

> library(Brobdingnag)

> factorial(170)
[1] 7.257416e+306

> as.numeric(prod(as.brob(1:170)))
[1] 7.257416e+306

   according to bc,

define fact(n) { if (n<2) return 1; return n*fact(n-1) }
fact(170)
72574156153079989673967282111292631147169916812964513765435777989005\
61843401706157852350749242617459511490991237838520776666022565442753\
02532890077320751090240043028005829560396661259965825710439855829425\
75689663134396122625710949468067112055688804571933402126614528000000\
00000000000000000000000000000000000


> prod(as.brob(1:1000))
[1] +exp(5912.1)

with bc -l

l(fact(1000))
5912.12817848816334887812

and in base R
> lfactorial(1000)
[1] 5912.128

  so it's not too bad ...

 (someone want to post this thread to the wiki ... ???)

  All of this does suppose that you have a good answer to Murray
Cooper's question of why you really *need* to evaluate factorials
this large and can't get away with manipulating them on the
log scale ...

  Ben Bolker

-- 
View this message in context: http://www.nabble.com/large-factorials-tp23175816p23189628.html
Sent from the R help mailing list archive at Nabble.com.