[R] Rmpfr question

[Martin Maechler]

>>>>> Ulises M Alvarez <uma at fata.unam.mx>
>>>>>     on Sat, 16 Oct 2010 14:34:25 -0500 writes:

    > Hi: I'm trying to reproduce an arbitrary precision
    > constant from 'Why and How to Use Arbitrary Precision'
    > (Ghazi et al., COMPUTING IN SCIENCE & ENGINEERING May/June
    > 2010; http://perso.ens-lyon.fr/philippe.theveny/cise.pdf):

    > d = 173746a + 94228b − 78487c

    > where: a = sin(1022), b = log(17.1), and c = exp(0.42).
      	         ^^^^^^^^^
  from what you do below, this should be  sin(10^22)

    > Ghazi et al. report: d = −1.341818958e−12 whit IEEE-754
    > quadruple precision (113 bits).

    > I have tried to reproduce such result using the Rmpfr
    > library using:

    > a <- mpfr(sin(10^22), 230)

if you think a bit, it's clear that this cannot work!
sin(10^22) is computed with only double precision, hence is
inaccurate, only using 53 bits.

    > b <- mpfr(log(171/10), 230)

the same here, of course

    > c <- mpfr(exp(42/100), 230)

and here.
    > (d <- 173746*a + 94228*b - 78487*c)

    > 1 'mpfr' number of precision  230   bits
    > [1] 
    > 2.9904079212883516447618603706359863281250000000000000000000000000000000e-11

    > Which does not correspond to the value reported by Ghazi et al.

    > Could any one give me some some hint on how to obtain a value closer to 
    > the one of −1.341818958e−12 please.

I hope I'm not doing your home work here.
The correct solution of course is

> a <- sin(mpfr(10, 230)^22)
> b <- log(mpfr(171, 230)/10)
> c <- exp(mpfr(42, 230)/100)
> (d <- 173746*a + 94228*b - 78487*c)
1 'mpfr' number of precision  230   bits 
[1] -1.3418189578296195497042786842309588809452366232139762407816198747581278e-12
> 

Best regards,
Martin Maechler, ETH Zurich

    > (Yes, I have reproduce the 'd' value from Ghazi et al. using  'gcc 
    > -lmpfr -lgmp' on linux boxes).

    > #################
    > Session info:
    > R version 2.11.1 (2010-05-31)
    > x86_64-apple-darwin9.8.0

    > locale:
    > [1] C

    > attached base packages:
    > [1] stats     graphics  grDevices utils     datasets  methods   base

    > other attached packages:
    > [1] Rmpfr_0.2-3

    > loaded via a namespace (and not attached):
    > [1] tools_2.11.1

    > (Have also tried with R on Linux boxes without success)
    > -- 
    > Ulises M. Alvarez
    > http://www.fata.unam.mx/