>>>>> 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/