vQ> if you really really need to have it done from within r,
vQ> you may want to use an external facility such as bc, the
vQ> 'basic calculator' [1,2]. for example, use the
vQ> (experimental!) r-bc:
vQ> source('http://r-bc.googlecode.com/svn/trunk/R/bc.R')
vQ> (you can also download the zipped package which will
vQ> install on windows, where you're likely not to have bc
vQ> yet; see http://code.google.com/p/r-bc/downloads/)
vQ> # an intuitive but slow approach implemented mostly in r
vQ> # (alternatively, you may want to have it recursive)
vQ> factorial.r = function(n) {
vQ> result = bc(1)
vQ> while (n > 1) {
vQ> result = result*n
vQ> n = n-1 }
vQ> result }
vQ> # an alternative, faster approach implemented mostly in bc
vQ> factorial.bc = function(n)
vQ> bc(sprintf('define fact(n) { if (n < 2) return 1; return n *
vQ> fact(n-1) }; fact(%d)', n))
vQ> library(rbenchmark)
vQ> benchmark(replications=10, columns=c('test', 'elapsed'),
vQ> r=factorial.r(500),
vQ> bc=factorial.bc(500))
vQ> # test elapsed
vQ> # 2 bc 0.101
vQ> # 1 r 34.181
vQ> this gives you factorials for arbitrary input, but note that the result
vQ> is not an integer, but an object of class 'bc' backed by a *character
vQ> string*:
vQ> result = factorial.bc(10^4)
vQ> is(result)
vQ> # "bc"
vQ> nchar(result)
vQ> # 35660
vQ> vQ
vQ> [1] http://www.gnu.org/software/bc/manual/html_mono/bc.html
vQ> [2] http://www.opengroup.org/onlinepubs/9699919799/utilities/bc.html
yet another alternative for arbitrary precision computing with R
is using using the GMP(http://www.gmp.org/) - based
MPFR(http://www.mpfr.org/) with the R package Rmpfr from
R-forge, http://r-forge.r-project.org/projects/rmpfr/
[ Unfortunately, it seems that there's no DLL version of the MPFR
library available for windows, and so the R-forge (and
Win-builder.r-project.org) maintainers currently do not build
Windows versions of the Rmpfr R package. ]
For the present case this brings no big advantage, as indeed the
factorial is trivial using multiprecision *integer* arithmetic;
The big advantage of MPFR and Rmpfr is the availability of many
transcendtal functions in multiprecision (arbitrary precision)
arithmetic.
First note
> factorial
function (x)
gamma(x + 1)
Then
> install.packages("Rmpfr", repos="http://R-Forge.R-project.org";)
> library(Rmpfr)
> gamma(as(1000,"mpfr"))
1 'mpfr' number of precision 128 bits
[1] 4.023872600770937735437024339230039857186e2564
> gamma(mpfr(1000, prec = 10000))
1 'mpfr' number of precision 10000 bits
[1]
4023872600770937735437024339230039857193748642107146325437999104299385123986290205920442084869694048004799886101971960586316668729948085589013238296699445909974245040870737599188236277271887325197795059509952761208749754624970436014182780946464962910563938874378864873371191810458257836478499770124766328898359557354325131853239584630755574091142624174743493475534286465766116677973966688202912073791438537195882498081268678383745597317461360853795345242215865932019280908782973084313928444032812315586110369768013573042161687476096758713483120254785893207671691324484262361314125087802080002616831510273418279777047846358681701643650241536913982812648102130927612448963599287051149649754199093422215668325720808213331861168115536158365469840467089756029009505376164758477284218896796462449451607653534081989013854424879849599533191017233555566021394503997362807501378376153071277619268490343526252000158885351473316117021039681759215109077880193931781141945452572238655414610628921879!
6022383897147608850627686296714667469756291123408243920816015378088989396451826324367161676217916890977991190375403127462228998800519544441428201218736174599264295658174662830295557029902432415318161721046583203678690611726015878352075151628422554026517048330422614397428693306169089796848259012545832716822645806652676995865268227280707578139185817888965220816434834482599326604336766017699961283186078838615027946595513115655203609398818061213855860030143569452722420634463179746059468257310379008402443243846565724501440282188525247093519062092902313649327349756551395872055965422874977401141334696271542284586237738753823048386568897646192738381490014076731044664025989949022222176590433990188601856652648506179970235619389701786004081188972991831102117122984590164192106888438712185564612496079872290851929681937238864261483965738229112312502418664935314397013742853192664987533721894069428143411852015801412334482801505139969429015348307764456909907315243327828826986460278986432113!
9083506217095002597389863554277196742822248757586765752344220207573630
569498825087968928162753848863396909959826280956121450994871701244516461260379029309120889086942028510640182154399457156805941872748998094254742173582401063677404595741785160829230135358081840096996372524230560855903700624271243416909004153690105933983835777939410970027753472000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
>
vQ> Murray Cooper wrote:
>> You don't say what the error was, for the R factorial function,
>> but it is probably irrelevant for your question.
>>
>> Factorials get to be big numbers rather quickly and unless you
>> are using a program that does arbitrary precission arithmetic
>> you will quickly exceed the precission limits, for storing a number.
>> If you have Maple, do 170! and count the number of digits in the
>> result. You will see what I mean.
>>
>> There are some tricks when working with large factorials, depending
>> on what you are doing with them. I'd first try the log factorial function
>> in R I think its called lfactorial. Just do a ?factorial and you'll find
>> documentation. If this doesn't work, for you, repost with a clear
>> description of what you're trying to do and someone may be able
>> to help.
>>
>> Murray M Cooper, Ph.D.
>> Richland Statistics
>> 9800 N 24th St
>> Richland, MI, USA 49083
>> Mail: richs...@earthlink.net
>>
>> ----- Original Message ----- From: "molinar" <sky...@hotmail.com>
>> To: <r-help@r-project.org>
>> Sent: Wednesday, April 22, 2009 3:21 PM
>> Subject: [R] large factorials
>>
>>
>>>
>>> I am working on a project that requires me to do very large factorial
>>> evaluations. On R the built in factorial function and the one I created
>>> both are not able to do factorials over 170. The first gives an
>>> error and
>>> mine return Inf.
>>>
>>> Is there a way to have R do these larger calculations (the calculator in
>>> accessories can do 10000 factorial and Maple can do even larger)
>>> --
>>> View this message in context:
>>> http://www.nabble.com/large-factorials-tp23175816p23175816.html
>>> Sent from the R help mailing list archive at Nabble.com.
>>>
>>> ______________________________________________
>>> 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.
vQ> ______________________________________________
vQ> R-help@r-project.org mailing list
vQ> https://stat.ethz.ch/mailman/listinfo/r-help
vQ> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
vQ> 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.
