[R] accurate log(1+x) {was "(arbitrary) precision"}

[Martin Maechler]

>>>>> "Gabor" == Gabor Grothendieck <ggrothendieck at myway.com>
>>>>>     on Fri, 18 Feb 2005 02:55:15 +0000 (UTC) writes:

    Gabor> Michael Grottke <Michael.Grottke <at> duke.edu>
    Gabor> writes: : I am currently using R for fitting a model
    Gabor> to various data sets : (minimizing the negative
    Gabor> log-likelihood) and calculating a number of : metrics
    Gabor> based on the parameter estimates. Within these
    Gabor> calculations, I : have steps of the form : :
    Gabor> log(log(1+x)), : : where x can be very small (e.g.,
    Gabor> around 1e-16). Unfortunately, the : precision of
    Gabor> doubles does not suffice for coping with the
    Gabor> difference in : the orders of magnitude of 1 and x:
    Gabor> 1+x is rounded to 1.  : : One way for solving this
    Gabor> problem seems to be to use an arbitrary : precision
    Gabor> library implemented in C and call the respective
    Gabor> routines for : calculating the logarithm(s) from
    Gabor> within R.  : : My questions are as follows: : 1. Is
    Gabor> there any better/more direct way to solve the
    Gabor> problem?  : 2. Is there any arbitrary precision
    Gabor> library you can suggest in particular?  :

    Gabor> The approximation log(1+x) = x would be accuate to
    Gabor> several decimal places in your case so your
    Gabor> expression would reduce to log(log(1+x)) = log(x).

but I think the most efficient solution for this case is to use

 log(log1p(x))

The point of the  log1p(.) function is to return  log(1+x)
accurately when ``for any x'', i.e., most importantly for x << 1.

Martin Maechler, ETH Zurich.