[R] problem with the precision of numbers

[(Ted Harding)]

On 19-Jan-10 18:48:47, Gabor Grothendieck wrote:
> 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_preci
>>> si
>>> 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.307011957285621299007119650203018618933170114453054867257099
> 2204603757660023189324468582740298425344"

Excellent reource! Thanks for pointing it out, Gabor.
Now for Kate Bush.

First, KB:
==========
Sweet and gentle sensitive man
With an obsessive nature and deep fascination
For numbers
And a complete infatuation with the calculation
Of PI

Oh he love, he love, he love
He does love his numbers
And they run, they run, they run him
In a great big circle
In a circle of infinity

3.1415926535 897932
3846 264 338 3279

Oh he love, he love, he love
He does love his numbers
And they run, they run, they run him
In a great big circle
In a circle of infinity
But he must, he must, he must
Put a number to it

50288419 716939937510
582319749 44 59230781
6406286208 821 4808651 32

Oh he love, he love, he love
He does love his numbers
And they run, they run, they run him
In a great big circle
In a circle of infinity

82306647 0938446095 505 8223?
========================================
KB she say:
3.
14159 26535 89793  23846    20
26433 83279 50288  41971    40
69399 37510 582*31* 97494   60
45923 07816 40628  6208||8  80
21480 86513 28230  66470   100
93844 60955 05822  3       116+1 = 117


Next, bc:
========
  source("http://r-bc.googlecode.com/svn/trunk/R/bc.R")
  bc("scale=200
  + 4*a(1)
  + ")
  [1]
"
3.
14159  26535 89793   23846
26433  83279 50288   41971
69399  37510 582*0*9 74944
59230  78164 06286   208|99
86280  34825 34211   70679
|82148 08651 32823   06647
09384  46095 50582   23

172
53594 08128 48111  74502
84102 70193 85211  05559
64462 29489 54930  38196
"
[edited for layout and indicators]

So KB replaces a "0" at decimal place 54 by "31",
and omits "99 86280 34825 34211  70679".

Maybe there were overwhelming poetic reasons for this.
Ted.

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