[R] problem with the precision of numbers

[(Ted Harding)]

On 24-Jan-10 22:24:10, kayj wrote:
> 
> 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.6484825769379918039202329231788093343835337105287241996630873753794810
> 915791302829474613042869191092322033403649086087872119205043462728841279
> 067042348e-6
> 
> Maple result
> 6.656852479*10^(-20)
> 
> why are the two results different?
> I appreciate your help

The result from bc (run "raw" and not from R, and assuming I
programmed it correctly) is (with some formatting for clarity):

s
0.
0000000000 0000000000 0000000000 0000000000 0000000000
0000000000 0000000000 0000000000 0000049078 2995761647
7217765991 5850974377 7588454525 1985604552 6041517025
6783796473 3395739785 9774327566 5079076761 9753976534

i.e.
4.90782995761647.......... * 10^(-86)

The 'bc' program was:

scale=200
define f(x) {
  if(x <= 1) return (1);
  return ( f(x-1)*x );
}
define choose(n,m) {
  return ( f(n)/(f(m)*f(n-m)) )
}
s=0
j=-1
for(i=0;i<=200;i++){
  c=choose(200,i);
  s=s+((j^i)*c*(1-(i+1)*1/200)^200);
}
s

(and, despite the recursive definition, the loop only took
about 3 seconds).

So my attempt came out different from your mpfr and your Maple.
It might be worth someone checking my 'bc' code! And, just in
case anyone suspects the accuracy of numbers like f(200) = 200!,
in R I get:

   print(lgamma(201),15)
  # [1] 863.231987192405

and in bc:
  l(f(200))
  #     863.2319871924054734957......................

so there is agreement there.

Ted.

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