Skip to content

large factorials

17 messages · molinar, Murray Cooper, John +7 more

#
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)
#
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
#
if you really really need to have it done from within r, you may want to
use an external facility such as bc, the 'basic calculator' [1,2].  for
example, use the (experimental!) r-bc:

    source('http://r-bc.googlecode.com/svn/trunk/R/bc.R')

(you can also download the zipped package which will install on windows,
where you're likely not to have bc yet; see
http://code.google.com/p/r-bc/downloads/)

    # an intuitive but slow approach implemented mostly in r
    # (alternatively, you may want to have it recursive)
    factorial.r = function(n) {
       result = bc(1)
       while (n > 1) {
          result = result*n
          n = n-1 }
       result }

    # an alternative, faster approach implemented mostly in bc
    factorial.bc = function(n)
       bc(sprintf('define fact(n) { if (n < 2) return 1; return n *
fact(n-1) }; fact(%d)', n))

    library(rbenchmark)
    benchmark(replications=10, columns=c('test', 'elapsed'),
       r=factorial.r(500),
       bc=factorial.bc(500))

    #   test elapsed
    # 2   bc   0.101
    # 1    r  34.181

this gives you factorials for arbitrary input, but note that the result
is not an integer, but an object of class 'bc' backed by a *character
string*:

    result = factorial.bc(10^4)
    is(result)
    # "bc"
    nchar(result)
    # 35660

vQ


[1] http://www.gnu.org/software/bc/manual/html_mono/bc.html
[2] http://www.opengroup.org/onlinepubs/9699919799/utilities/bc.html
Murray Cooper wrote:
#
On Wednesday 22 April 2009 12:21:41 pm molinar 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
#
J Dougherty wrote:
???  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 ...
[1] 7.257416e+306
[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
[1] +exp(5912.1)

with bc -l

l(fact(1000))
5912.12817848816334887812

and in base R
[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
#
J Dougherty wrote:
octave works similarly to *matlab* (to the point some call it a 'matlab
clone').  it's not so close to mathematica.

vQ
#
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
function (x) 
gamma(x + 1)

Then
1 'mpfr' number of precision  128   bits 
[1] 4.023872600770937735437024339230039857186e2564
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.
#
Also the R sympy package can handle this:
Loading required package: rJava
[1] "3628800"
+       bc=factorial.bc(500),
+       sympy = factorial.sympy(500))
   test elapsed
1    bc    2.17
2 sympy    0.09

See the rSymPy, r-bc and rbenchmark home pages:
http://rsympy.googlecode.com
http://r-bc.googlecode.com
http://rbenchmark.googlecode.com
On Wed, Apr 22, 2009 at 3:21 PM, molinar <sky2k2 at hotmail.com> wrote:
#
Thank you everyone all of your posts were very helpful.  I tried each one and
I think I have about 10 new packages installed.  The formula I need to
calculate did not involve any logarithms but infinite summations of
factorials, I'm sorry for not specifying.  I read some things about using
logarithms but I thought in my case I would have to do an e to the log and
by doing that R still gave me the same problems with numbers over 170.  

But I was able to get it to work by using the last post about the rsympy
packages.  

I tried downloading bc but I didn't know how to connect it to R, so R said
"could not find function bc".  

Thanks again for all of your help.
Samantha
Gabor Grothendieck wrote:

  
    
#
Here is what I did:
library(rSymPy)
factorial.sympy <- function(n) sympy(paste("factorial(", n, ")"))
factorial.sympy(171) 
[1]
"1241018070217667823424840524103103992616605577501693185388951803611996075221691752992751978120487585576464959501670387052809889858690710767331242032218484364310473577889968548278290754541561964852153468318044293239598173696899657235903947616152278558180061176365108428800000000000000000000000000000000000000000"
Which work perfectly. 

Here is one of my summation functions:

sum1 <- function(l,u,t,i,n,w) {
+ v <- 0
+ for (m in 0 :w) {
+ v1 <- ((u^(1/2))*(l^(1/2))*t)^(i-n+2*m)
+ v2 <- (factorial.sympy(i-n+m))*(factorial.sympy(m))
+ v3 <- v1/v2
+ v <- v+v3
+ }
+ return(v)
+ }

sum1(1,2,10,80,3,80)
Error in (factorial.sympy(i - n + m)) * (factorial.sympy(m)) : 
  non-numeric argument to binary operator

I'm not sure why it works when I do the factorial normally but when I call
my function it doesn't work?
molinar wrote:

  
    
#
Hi Samantha,

It is quite likely that you are not doing something right when you are
explicitly computing large factorials.  There is probably a good asymptotic
approximation that will simplify things for you.  In my experience, there is
seldom a need to explicitly compute factorials of integers.  What is the
real problem that you are trying to solve?

Ravi.


----------------------------------------------------------------------------
-------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan at jhmi.edu

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html



----------------------------------------------------------------------------
--------


-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of molinar
Sent: Thursday, April 23, 2009 9:42 AM
To: r-help at r-project.org
Subject: Re: [R] large factorials


Thank you everyone all of your posts were very helpful.  I tried each one
and I think I have about 10 new packages installed.  The formula I need to
calculate did not involve any logarithms but infinite summations of
factorials, I'm sorry for not specifying.  I read some things about using
logarithms but I thought in my case I would have to do an e to the log and
by doing that R still gave me the same problems with numbers over 170.  

But I was able to get it to work by using the last post about the rsympy
packages.  

I tried downloading bc but I didn't know how to connect it to R, so R said
"could not find function bc".  

Thanks again for all of your help.
Samantha
Gabor Grothendieck wrote:
--
View this message in context:
http://www.nabble.com/large-factorials-tp23175816p23197201.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.
#
sympy() returns a character string, not an R numeric -- it shouldn't
automatically return an R numeric since R can't represent all
the numbers that sympy can.

The development version of rSymPy has a second class which
produces objects of class c("Sym", "character") and those
can be manipulated with +, -, *, / producing other Sym
objects so try this:


library(rSymPy)

# next line pulls in code to handle Sym objects
source("http://rsympy.googlecode.com/svn/trunk/R/Sym.R")

# define factorial to return a Sym object
factorial.Sym <- function(n) Sym("factorial(", n, ")")

sum1 <- function(l,u,t,i,n,w) {
 v <- 0
 for (m in 0 :w) {
	 v1 <- ((u^(1/2))*(l^(1/2))*t)^(i-n+2*m)
	 v2 <- (factorial.Sym(i-n+m))*(factorial.Sym(m))
	 v3 <- v1/v2
	 v <- v+v3
 }
 sympy(v)
}

s <- sum1(1,2,10,80,3,80)
s # Sym object
as.numeric(s) # numeric
On Thu, Apr 23, 2009 at 10:00 AM, molinar <sky2k2 at hotmail.com> wrote:
#
The code in my prior post works (except one comment was wrong)
but try this instead.  The only change is the last line of the sum1
function.   This way it produces a Sym object rather than a character
string.

library(rSymPy)

# define factorial to return a Sym object
factorial.Sym <- function(n) Sym("factorial(", n, ")")

sum1 <- function(l,u,t,i,n,w) {
 v <- 0
 for (m in 0 :w) {
	 v1 <- ((u^(1/2))*(l^(1/2))*t)^(i-n+2*m)
	 v2 <- (factorial.Sym(i-n+m))*(factorial.Sym(m))
	 v3 <- v1/v2
	 v <- v+v3
 }
 Sym(sympy(v)) # send it to SymPy, make result a Sym obj
}

s <- sum1(1,2,10,80,3,80)
s # Sym
as.numeric(s) # numeric


On Thu, Apr 23, 2009 at 10:37 AM, Gabor Grothendieck
<ggrothendieck at gmail.com> wrote:
#
One more improvement.  Perhaps it would be best just to return a numeric
so sum1 inputs and outputs numerics:

library(rSymPy)

# define factorial to return a Sym object
factorial.Sym <- function(n) Sym("factorial(", n, ")")

sum1 <- function(l,u,t,i,n,w) {
 v <- 0
 for (m in 0 :w) {
        v1 <- ((u^(1/2))*(l^(1/2))*t)^(i-n+2*m)
        v2 <- (factorial.Sym(i-n+m))*(factorial.Sym(m))
        v3 <- v1/v2
        v <- v+v3
 }

 as.numeric(sympy(v)) # send it to SymPy, make result a Sym obj

}

s <- sum1(1,2,10,80,3,80)
s # numeric


On Thu, Apr 23, 2009 at 10:44 AM, Gabor Grothendieck
<ggrothendieck at gmail.com> wrote:
#
On Wed, Apr 22, 2009 at 08:26:51PM -0700, Ben Bolker wrote:

            
You might take a look at Sage.  It is not a mathematica clone, but 
open source mathematical software which has a development community 
similar to that of R.   

        http://www.sagemath.org/

albyn
#
Albyn Jones wrote:
... except for that ws has quite a different attitude than bdr.

vQ
3 days later