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.