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

Gabor Grothendieck ggrothendieck at gmail.com
Sat Jul 4 16:42:43 CEST 2015 

You can do that with bc if you pass the entire expression to bc(...) in
quotes but in that case you will have to use bc notation, not R notation
so, for example, exp is e and atan is a.

> library(bc)
> bc("e(sqrt(163)*4*a(1))")
[1]
"262537412640768743.9999999999992500725971981856888793538563373369908627075374103782106479101186073116295306145602054347"


On Fri, Jul 3, 2015 at 3: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
>
>
> ______________________________________________
> 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.
>



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

	[[alternative HTML version deleted]]

More information about the R-help mailing list