[R] problem with the precision of numbers

[(Ted Harding)]

On 25-Jan-10 00:12:29, William Dunlap wrote:
>> -----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

Bill, using your code strategy I agree with your mpfr result
(and with Maple) using 'bc:

s=0
j = -1
c = 1
for (i=0; i <= 80; i++ ){
        s = s + ((j^i) * c * (1 - (i + 1) * 1/200)^200)
        c = (c * (80 - i))/(i + 1)
    }
s
0.
0000000000 0000000006 6568524789 6620376996 7327577118
2135789699 5101929940 0816838328 9634271996 1257038416 
6039147205 2215051658 2822460868 2286772321 8214397124
2739979506 9906258378 1667597096 6899329803 4460676207

(time to look and see what I did wrong before ...)
Ted.

>> 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.
>> 
> 
> ______________________________________________
> 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.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 25-Jan-10                                       Time: 00:47:53
------------------------------ XFMail ------------------------------