[R] problem with the precision of numbers

[William Dunlap]

> -----Original Message-----
> From: r-help-bounces at r-project.org 
> [mailto:r-help-bounces at r-project.org] On Behalf Of kayj
> Sent: Sunday, January 24, 2010 2:24 PM
> To: r-help at r-project.org
> Subject: [R] problem with the precision of numbers
> 
> 
> Hi All,
> 
> I was wondering if R can deal with high precsion numbers, below is an
> example that I tried on using R and Maple where I got 
> different results. I
> was not able to use the r-bc package in R, instead I used the Rmpfr
> package, below are both R and Maple results
> 
> 
> > library(Rmpfr) 
> > 
> > s<-mpfr(0,500)
> > j <- mpfr(-1,500)
> > 
> > for (i in 0:80){ 
> + p<-as.numeric(i)
> + c<-choose(80,i)
> + s=s+((j^i)*c*(1-(i+1)*1/200)^200)
> +  
> + }
> > s
> 1 'mpfr' number of precision  500   bits 
> [1]
> 4.648482576937991803920232923178809334383533710528724199663087
> 37537948109157913028294746130428691910923220334036490860878721
> 19205043462728841279067042348e-6

You are computing (1-(i+1)*1/200)^200 with ordinary
32-bit integers and 64-bit double precision numbers.
The same goes for choose(80,i).  Try something like

f <- function (precision) 
{
    require(Rmpfr)
    s <- mpfr(0, precision)
    j <- mpfr(-1, precision)
    c <- mpfr(1, precision)
    for (i in 0:80) {
        i <- mpfr(i, precision)
        s <- s + ((j^i) * c * (1 - (i + 1) * 1/200)^200)
        c <- (c * (80 - i))/(i + 1)
    }
    s
}
> print(f(precision=500), digits=25)
1 'mpfr' number of precision  500   bits 
[1] 6.656852478966203769967328e-20

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com 

> 
> 
> Maple result
> 
> 6.656852479*10^(-20)
> 
> 
> why are the two results different?
> 
> I appreciate your help
> 
> 
> 
> 
> -- 
> View this message in context: 
> http://n4.nabble.com/problem-with-the-precision-of-numbers-tp1
288905p1288905.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.
>