On Tue, Jan 19, 2010 at 1:41 PM, Ted Harding <ted.hard...@manchester.ac.uk> wrote: > On 19-Jan-10 17:55:43, Ben Bolker wrote: >> kayj <kjaja27 <at> yahoo.com> writes: >>> Hi All, >>> >>> I was wodering if it is possible to increase the precision using R. >>> I ran the script below in R and MAPLE and I got different results >>> when k is large. >>> Any idea how to fix this problem? thanks for your help >>> >>> for (k in 0:2000){ >>> s=0 >>> for(i in 0:k){ >>> s=s+((-1)^i)*3456*(1+i*1/2000)^3000 >>> } >>> } >> >> (1) see >> http://wiki.r-project.org/rwiki/doku.php?id=misc:r_accuracy:high_precisi >> on_arithmetic >> >> (2) consider whether there is more accurate algorithm you >> could use. I don't recognize the series, but perhaps it >> has a closed form solution, maybe as a special function? >> How much accuracy do you really need in the solution? >> >> Ben Bolker > > I suspect this is an invented computation -- the "3456" strikes > me as "unlikely" (it reminds me of my habitual illustrative use > of set.seed(54321)). > > There is a definite problem with the development given by kayj. > When k=2000 and i=k, the formula requires evaluation of > > 3456*(2^3000) > > on a log10 scale this is > > log10(3456) + 3000*log10(2) = 906.6286 > > Since R "gives up" at 10^308.25471 = 1.79767e+308 > (10^308.25472 => Inf), this algorithm is going to be tricky to > evaluate! > > I don't know how well Rmpfr copes with very large numbers (the > available documentation seems cryptic). However, I can go along > with the recommendation in the URL the Ben gives, to use 'bc' > ("Berkeley Calculator"), available on unix[oid] systems since > a long time ago. That is an old friend of mine, and works well > (it can cope with exponents up to X^2147483647 in the version > I have). It can eat for breakfast the task of checking whether > Kate Bush can accurately sing pi to 117 significant figures: > > http://www.absolutelyrics.com/lyrics/view/kate_bush/pi > > (Try it in R). >
There is an R interface to bc here at http://r-bc.googlecode.com . Trying it for k up to 10: > source("http://r-bc.googlecode.com/svn/trunk/R/bc.R") > bc("for (k = 0; k <= 10; k = k + 1) { + s=0 + for (i = 0; i <= k; i = i + 1) { + s=s+((-1)^i)*3456*(1+i*1/2000)^3000 + } + } + s + ") [1] "8886117368.3070119572856212990071196502030186189331701144530548672570992204603757660023189324468582740298425344" ______________________________________________ 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.