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

________________________________________
From: David Winsemius <dwinsemius at comcast.net>
Sent: Saturday, July 4, 2015 1:12 PM
To: Duncan Murdoch
Cc: r-help; John Nash; Ravi Varadhan
Subject: Re: [R] : Ramanujan and the accuracy of floating point computations - using Rmpfr in R

On Jul 4, 2015, at 12:20 AM, Duncan Murdoch <murdoch.duncan at gmail.com> wrote:

On 04/07/2015 8:21 AM, David Winsemius wrote:
On Jul 3, 2015, at 11:05 PM, Duncan Murdoch <murdoch.duncan at gmail.com> wrote:

On 04/07/2015 3:45 AM, David Winsemius wrote:
On Jul 3, 2015, at 5:08 PM, David Winsemius <dwinsemius at comcast.net> wrote:



It doesn?t appear to me that mpfr was ever designed to handle expressions as the first argument.
This could be a start. Obviously one would wnat to include code to do other substitutions probably using the all.vars function to pull out the other ?constants? and ?numeric? values to make them of equivalent precision. I?m guessing you want to follow the parse-tree and then testing the numbers for integer-ness and then replacing by paste0( ?mpfr(?, val, ?L, ?, prec,?)? )

Pre <- function(expr, prec){ sexpr <- deparse(substitute(expr) )
Why deparse?  That's almost never a good idea.  I can't try your code (I
don't have mpfr available), but it would be much better to modify the
expression than the text representation of it.  For example, I think
your code would modify strings containing "pi", or variables with those
letters in them, etc.  If you used substitute(expr) without the
deparse(), you could replace the symbol "pi" with the call to the Const
function, and be more robust.

Really? I did try. I was  fairly sure that someone could do better but I don?t see an open path along the lines you suggest. I?m pretty sure I tried `substitute(expr, list(pi= pi))` when `expr` had been the formal argument and got disappointed because there is no `pi` in the expression `expr`. I _thought_ the problem was that `substitute` does not evaluate its first argument, but I do admit to be pretty much of a klutz with this sort of programming. I don?t think you need to have mpfr installed in order to demonstrate this.
The substitute() function really does two different things.
substitute(expr) (with no second argument) grabs the underlying
expression out of a promise.  substitute(expr, list(pi = pi)) tries to
make the substitution in the expression "expr", so it doesn't see "pi?.
Thank you. That was really helpful. I hope it ?sticks? to  sufficiently durable set of neurons.
This should work:

do.call(substitute, list(expr = substitute(expr), env=list(pi =
Const(pi, 120))))

(but I can't evaluate the Const function to test it).
The expression `pi` as the argument to Const only needed to be character()-ized and then evaluation on that result succeeded:

library(mpfr)
#  Pre <- function(expr, prec){ do.call(substitute,
                                        list(expr = substitute(expr),
                                             env=list(pi = Const(pi, prec)))) }

Pre(exp(pi),120)
Error in Const(pi, prec) :
  'name' must be one of 'pi', 'gamma', 'catalan', 'log2'


 Pre <- function(expr, prec){ do.call(substitute,
                                     list(expr = substitute(expr),
                                          env=list(pi = Const('pi', prec)))) }

Pre(exp(pi),120)
exp(list(<S4 object of class "mpfr1">))
eval( Pre(exp(pi),120) )
1 'mpfr' number of precision  120   bits
[1] 23.140692632779269005729086367948547394

So the next step for delivering Ravi?s request would be to add the rest of the ?Const? options:

Pre <- function(expr, prec){ do.call(substitute,
                                     list(expr = substitute(expr),
                                         env=list(pi = Const('pi', prec),
                                         catalan=Const('catalan', prec),
                                         log2=Const('log2', prec),
                                         gamma=Const('gamma', prec)))) }


eval(Pre(exp(gamma)))
Error in stopifnot(is.numeric(prec)) :
  argument "prec" is missing, with no default
eval(Pre(exp(gamma), 120))
1 'mpfr' number of precision  120   bits
[1] 1.781072417990197985236504103107179549

This works:

do.call(substitute, list(expr = substitute(expr), env = list(pi =
quote(Const(pi, 120)))))

because it never evaluates the Const function, it just returns some code
that you could modify more if you liked.

Duncan Murdoch
David Winsemius, MD
Alameda, CA, USA