[R] large factorials

Martin Maechler maechler at stat.math.ethz.ch
Thu Apr 23 10:45:07 CEST 2009 

    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] 402387260077093773543702433923003985719374864210714632543799910429938512398629020592044208486969404800479988610197196058631666872994808558901323829669944590997424504087073759918823627727188732519779505950995276120874975462497043601418278094646496291056393887437886487337119181045825783647849977012476632889835955735432513185323958463075557409114262417474349347553428646576611667797396668820291207379143853719588249808126867838374559731746136085379534524221586593201928090878297308431392844403281231558611036976801357304216168747609675871348312025478589320767169132448426236131412508780208000261683151027341827977704784635868170164365024153691398281264810213092761244896359928705114964975419909342221566832572080821333186116811553615836546984046708975602900950537616475847728421889679646244945160765353408198901385442487984959953319101723355556602139450399736280750137837615307127761926849034352625200015888535147331611702103968175921510907788019393178114194545257223865541461062892187960223838971476088506276862967146674697562911234082439208160153780889893964518263243671616762179168909779911903754031274622289988005195444414282012187361745992642956581746628302955570299024324153181617210465832036786906117260158783520751516284225540265170483304226143974286933061690897968482590125458327168226458066526769958652682272807075781391858178889652208164348344825993266043367660176999612831860788386150279465955131156552036093988180612138558600301435694527224206344631797460594682573103790084024432438465657245014402821885252470935190620929023136493273497565513958720559654228749774011413346962715422845862377387538230483865688976461927383814900140767310446640259899490222221765904339901886018566526485061799702356193897017860040811889729918311021171229845901641921068884387121855646124960798722908519296819372388642614839657382291123125024186649353143970137428531926649875337218940694281434118520158014123344828015051399694290153483077644569099073152433278288269864602789864321139083506217095002597389863554277196742822248757586765752344220207573630569498825087968928162753848863396909959826280956121450994871701244516461260379029309120889086942028510640182154399457156805941872748998094254742173582401063677404595741785160829230135358081840096996372524230560855903700624271243416909004153690105933983835777939410970027753472000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.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: richstat at earthlink.net
    >> 
    >> ----- Original Message ----- From: "molinar" <sky2k2 at hotmail.com>
    >> To: <r-help at 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 at 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 at 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 at 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.

More information about the R-help mailing list