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.hard...@manchester.ac.uk> Fax-to-email: +44 (0)870 094 0861 Date: 24-Jan-10 Time: 23:54:54 ------------------------------ XFMail ------------------------------ ______________________________________________ R-help@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.