[R] : Ramanujan and the accuracy of floating point computations - using Rmpfr in R

David Winsemius dwinsemius at comcast.net
Sat Jul 4 02:08:03 CEST 2015 


It doesn’t appear to me that mpfr was ever designed to handle expressions as the first argument.

— 
David

> On Jul 3, 2015, at 12:01 PM, Ravi Varadhan <ravi.varadhan at jhu.edu> wrote:
> 
> Thank you all.  I did think about declaring `pi' as a special constant, but for some reason didn't actually try it.  
> Would it be easy to have the mpfr() written such that its argument is automatically of extended precision? In other words, if I just called:  mpfr(exp(sqrt(163)*pi, 120), then all the constants, 163, pi, are automatically of 120 bits precision.
> 
> Is this easy to do?
> 
> Best,
> Ravi
> ________________________________________
> From: David Winsemius <dwinsemius at comcast.net>
> Sent: Friday, July 3, 2015 2:06 PM
> To: John Nash
> Cc: r-help; Ravi Varadhan
> Subject: Re: [R] : Ramanujan and the accuracy of floating point computations - using Rmpfr in R
> 
>> On Jul 3, 2015, at 8:08 AM, John Nash <John.Nash at uottawa.ca> wrote:
>> 
>> 
>> 
>> 
>> Third try -- I unsubscribed and re-subscribed. Sorry to Ravi for extra
>> traffic.
>> 
>> In case anyone gets a duplicate, R-help bounced my msg "from non-member" (I changed server, but not email yesterday,
>> so possibly something changed enough). Trying again.
>> 
>> JN
>> 
>> I got the Wolfram answer as follows:
>> 
>> library(Rmpfr)
>> n163 <- mpfr(163, 500)
>> n163
>> pi500 <- mpfr(pi, 500)
>> pi500
>> pitan <- mpfr(4, 500)*atan(mpfr(1,500))
>> pitan
>> pitan-pi500
>> r500 <- exp(sqrt(n163)*pitan)
>> r500
>> check <- "262537412640768743.99999999999925007259719818568887935385..."
>> savehistory("jnramanujan.R")
>> 
>> Note that I used 4*atan(1) to get pi.
> 
> RK got it right by following the example in the help page for mpfr:
> 
> Const("pi", 120)
> 
> The R `pi` constant is not recognized by mpfr as being anything other than another double .
> 
> 
> There are four special values that mpfr recognizes.
> 
> —
> Best;
> David
> 
> 
>> It seems that may be important,
>> rather than converting.
>> 
>> JN
>> 
>> On 15-07-03 06:00 AM, r-help-request at r-project.org wrote:
>> 
>>> Message: 40
>>> Date: Thu, 2 Jul 2015 22:38:45 +0000
>>> From: Ravi Varadhan <ravi.varadhan at jhu.edu>
>>> To: "'Richard M. Heiberger'" <rmh at temple.edu>, Aditya Singh
>>>     <aps6dl at yahoo.com>
>>> Cc: r-help <r-help at r-project.org>
>>> Subject: Re: [R] : Ramanujan and the accuracy of floating point
>>>     computations - using Rmpfr in R
>>> Message-ID:
>>>     <14ad39aaf6a542849bbf3f62a0c2f38f at DOM-EB1-2013.win.ad.jhu.edu>
>>> Content-Type: text/plain; charset="utf-8"
>>> 
>>> Hi Rich,
>>> 
>>> The Wolfram answer is correct.
>>> http://mathworld.wolfram.com/RamanujanConstant.html
>>> 
>>> There is no code for Wolfram alpha.  You just go to their web engine and plug in the expression and it will give you the answer.
>>> http://www.wolframalpha.com/
>>> 
>>> I am not sure that the precedence matters in Rmpfr.  Even if it does, the answer you get is still wrong as you showed.
>>> 
>>> Thanks,
>>> Ravi
>>> 
>>> -----Original Message-----
>>> From: Richard M. Heiberger [mailto:rmh at temple.edu]
>>> Sent: Thursday, July 02, 2015 6:30 PM
>>> To: Aditya Singh
>>> Cc: Ravi Varadhan; r-help
>>> Subject: Re: [R] : Ramanujan and the accuracy of floating point computations - using Rmpfr in R
>>> 
>>> There is a precedence error in your R attempt.  You need to convert
>>> 163 to 120 bits first, before taking
>>> its square root.
>>> 
>>>>> exp(sqrt(mpfr(163, 120)) * mpfr(pi, 120))
>>> 1 'mpfr' number of precision  120   bits
>>> [1] 262537412640768333.51635812597335712954
>>> 
>>> ## just the last four characters to the left of the decimal point.
>>>>> tmp <- c(baseR=8256, Wolfram=8744, Rmpfr=8333, wrongRmpfr=7837)
>>>>> tmp-tmp[2]
>>>    baseR    Wolfram      Rmpfr wrongRmpfr
>>>     -488          0       -411       -907
>>> You didn't give the Wolfram alpha code you used.  There is no way of verifying the correct value from your email.
>>> Please check that you didn't have a similar precedence error in that code.
>>> 
>>> Rich
>>> 
>>> 
>>> On Thu, Jul 2, 2015 at 2:02 PM, Aditya Singh via R-help <r-help at r-project.org> wrote:
>>>>> Ravi
>>>>> 
>>>>> I am a chemical engineer by training. Is there not something like law of corresponding states in numerical analysis?
>>>>> 
>>>>> Aditya
>>>>> 
>>>>> 
>>>>> 
>>>>> ------------------------------
>>>>> On Thu 2 Jul, 2015 7:28 AM PDT Ravi Varadhan wrote:
>>>>> 
>>>>>>> Hi,
>>>>>>> 
>>>>>>> Ramanujan supposedly discovered that the number, 163, has this interesting property that exp(sqrt(163)*pi), which is obviously a transcendental number, is real close to an integer (close to 10^(-12)).
>>>>>>> 
>>>>>>> If I compute this using the Wolfram alpha engine, I get:
>>>>>>> 262537412640768743.99999999999925007259719818568887935385...
>>>>>>> 
>>>>>>> When I do this in R 3.1.1 (64-bit windows), I get:
>>>>>>> 262537412640768256.0000
>>>>>>> 
>>>>>>> The absolute error between the exact and R's value is 488, with a relative error of about 1.9x10^(-15).
>>>>>>> 
>>>>>>> In order to replicate Wolfram Alpha, I tried doing this in "Rmfpr" but I am unable to get accurate results:
>>>>>>> 
>>>>>>> library(Rmpfr)
>>>>>>> 
>>>>>>> 
>>>>>>>>> exp(sqrt(163) * mpfr(pi, 120))
>>>>>>> 1 'mpfr' number of precision  120   bits
>>>>>>> 
>>>>>>> [1] 262537412640767837.08771354274620169031
>>>>>>> 
>>>>>>> The above answer is not only inaccurate, but it is actually worse than the answer using the usual double precision.  Any thoughts as to what I am doing wrong?
>>>>>>> 
>>>>>>> Thank you,
>>>>>>> Ravi
>>>>>>> 
>>>>>>> 
>>>>>>> 
>>>>>>>     [[alternative HTML version deleted]]
>>>>>>> 
>>>>>>> ______________________________________________
>>>>>>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>>>>>>> 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.
>>>>> ______________________________________________
>>>>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>>>>> 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.
>>> ------------------------------
>> 
>> 
>> ______________________________________________
>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> 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.
> 
> David Winsemius, MD
> Alameda, CA, USA
> 

David Winsemius, MD
Alameda, CA, USA

More information about the R-help mailing list