[R] problem with the precision of numbers

[Gabor Grothendieck]

On Tue, Jan 19, 2010 at 1:41 PM, Ted Harding
<Ted.Harding at 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"