What build of R can't calculate factorial(150)? I thought the max double of a 32 bit build would be on the order of 10^308 ~ 2^1024 (the material below seems to agree with me on this)
But yes, agreed with Ted: it's very helpful to think a bit about what you are calculating when doing these sorts of large number calculations before asking for the accurate calculation of huge numbers. Michael On Fri, Nov 18, 2011 at 4:31 PM, Ted Harding <ted.hard...@wlandres.net> wrote: > On 18-Nov-11 16:03:44, Gyanendra Pokharel wrote: >> Hi all, >> why factorial(150) shows the error out of range in 'gammafn'? >> I have to calculate the number of subset formed by 150 samples >> taking 10 at a time. How is this possible? >> best > > Because factorial(150) is nearly 10^263, which is far greater > than any number that R is capable of storing. > > I would, perhaps, suggest you work with lgamma(), the logarithm > (to base e) of the Gamma function. > > Thus factorial(150) = gamma(151), and > > lgamma(151) > # [1] 605.0201 > > so factorial(150) is close to e^605; to get it as a power of 10: > > lgamma(151)*log10(exp(1)) > # [1] 262.7569 > > > Hence the "nearly 10^263" above. > > If your question means "the number of ways of choosing 10 out > of 150, then, using lgamma(), its log to base e is: > > lgamma(151) - lgamma(11) - lgamma(141) > # [1] 34.6954 > > and its log to base 10 is: > > (lgamma(151) - lgamma(11) - lgamma(141))*log10(exp(1)) > # [1] 15.06802 > > so the result is about 10^15 which *is* within R's range. > > Hence: > > 10^( (lgamma(151) - lgamma(11) - lgamma(141))*log10(exp(1)) ) > # 1.169554e+15 > > You can see this (approximately) in an all-digits form by > using print(): > > X <- 10^( (lgamma(151) - lgamma(11) - lgamma(141))*log10(exp(1)) ) > print(X,18) > # [1] 1169554298222353 > > However, in many cases (such as your example) it will be feasible > to exploit the reduction in numbers of factors if you cancel out > the factors in the larger denominator factorial from the numerator > factorial, so that the number you want is > > 150*149*148*147*146*145*144*143*142*141/(10*9*8*7*6*5*4*3*2*1) > > = (150/10)*(149/9)*(148/8)* ... *(142/2)*(141/1) > > for which you can define a function nCr(n,r): > > nCr <- function(n,r){ > num <- (n:(n-r+1)) > den <- (r:1) > prod(num/den) > } > > and then: > > print(nCr(150,10),18) > # [1] 1169554298222310 > > which is not quite the same as the result 1169554298222353 > obtained above using lgamma(). This is because the previous > result is affected by rounding errors occuring in the > logarithms. The result obtained using the function nCr() is > exact. However, it can fail in its turn if you want one of > the larger possible results, such as nCr(1000,500): > > nCr(1500,750) > # [1] Inf > > since the result is too large for R to store (the largest on > my machine is about 2*10^308, see below)), and the products > in the function definition accumulate until this number is > exceeded. > > You can find out the limits of R storage on your machine > by using the command > > .Machine > > which gives a long list of the characteristics of the machine > you are using. In particular, you will see: > > $double.xmax > [1] 1.797693e+308 > > is the largest number that R will store; and: > > $double.eps > [1] 2.220446e-16 > > which is the smallest positive number; and: > > $double.neg.eps > [1] 1.110223e-16 > > which is the smallest negative number (note that the latter > is half the former). > > They can be individually accessed as: > > .Machine$double.xmax > # [1] 1.797693e+308 > > .Machine$double.eps > # [1] 2.220446e-16 > > .Machine$double.neg.eps > # [1] 1.110223e-16 > > Hoping this helps! > Ted. > > -------------------------------------------------------------------- > E-Mail: (Ted Harding) <ted.hard...@wlandres.net> > Fax-to-email: +44 (0)870 094 0861 > Date: 18-Nov-11 Time: 21:31:04 > ------------------------------ XFMail ------------------------------ > > ______________________________________________ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.