[R] large factorials

[Gabor Grothendieck]

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:
> 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:
>>
>> 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:
>>>
>>> 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:
>>>>
>>>> Also the R sympy package can handle this:
>>>>
>>>>> library(rSymPy)
>>>> Loading required package: rJava
>>>>
>>>>> factorial.sympy <- function(n) sympy(paste("factorial(", n, ")"))
>>>>
>>>>> # note that first time sympy is called it loads java, jython and sympy
>>>>> # but on subsequent calls its faster.  So make a dummy call first.
>>>>> factorial.sympy(10)
>>>> [1] "3628800"
>>>>
>>>>> # code from earlier post defining factorial.bc to be inserted here
>>>>
>>>>>    benchmark(replications=10, columns=c('test', 'elapsed'),
>>>> +       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:
>>>>>
>>>>> 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.
>>>>
>>>>
>>>
>>>
>>
>> --
>> View this message in context: http://www.nabble.com/large-factorials-tp23175816p23197344.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.
>>
>