From ripley at stats.ox.ac.uk  Tue Sep  2 14:21:30 2014
From: ripley at stats.ox.ac.uk (Prof Brian Ripley)
Date: Tue, 2 Sep 2014 13:21:30 +0100
Subject: [Rd] ggplot2/plyr interaction with latest R-devel?
In-Reply-To: <D800863A-31C6-4361-BB5C-036AD5D7C642@gmail.com>
References: <5404BE1D.2020609@mcmaster.ca>	<CAN_1p9ygzA+Kbb45QKqgASWE9vg11DV0=tLruSsMxWHULfYoCA@mail.gmail.com>	<CABghstQyz36aboSegmTbKbETGnncK4oc8rc68d6uVcMbpkKOQw@mail.gmail.com>	<CABghstQQdS1zgEBjW_cSYAY2M-M3cAQ_oUoUQRY6X80yPWLZQg@mail.gmail.com>
	<D800863A-31C6-4361-BB5C-036AD5D7C642@gmail.com>
Message-ID: <5405B64A.3010408@stats.ox.ac.uk>

On 02/09/2014 12:43, peter dalgaard wrote:
> Impressive. Never ceases to amaze me what computers can do these days. ;-)
>
> It's even more impressive given that we have
>
> static double logbase(double x, double base)
> {
> #ifdef HAVE_LOG10
>      if(base == 10) return x > 0 ? log10(x) : x < 0 ? R_NaN : R_NegInf;
> #endif
> #ifdef HAVE_LOG2
>      if(base == 2) return x > 0 ? log2(x) : x < 0 ? R_NaN : R_NegInf;
> #endif
>      return R_log(x) / R_log(base);
> }
>
> which, except possibly for base-10 and base-2, would seem to be quite hard to convince to return anything other than 1 if x == base....

I suspect PD knows (or has known) why, but for the sake of those who are 
not much bitten by ix86 compilers ....

It's the curse of extended-precision arithmetic (and not enough 
registers).  It does

compute in an FPU register, then store z1 = R_log(x)
compute z2 = R_log(base) in the FPU
load z1 to an FPU register.
compute z1/z2

(at least at -O2 with the version of gcc whose assembler output I looked 
at).  So z1 is has a 64-bit mantissa and z2 a 53-bit one.

x86_64 has more registers, and so avoids the store.  This can be 
circumvented on i686 by compiling with -msse2, which for some reason 
Linux distros do not make the default.

floor(log2(x)/3) is more reliable (and this is the underlying reason why 
we have log10 and log2).

>
> -pd
>
>
> On 02 Sep 2014, at 03:27 , Ben Bolker <bbolker at gmail.com> wrote:
>
>> log(8, base=8L)-1
>> log(8, base=8)-1
>> logvals <- setNames(log(2:25,base=2:25)-1,2:25)
>> logvals[logvals!=0]  ## 5,8,14,18,19,25 all == .Machine$double.eps/2
>


-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Emeritus Professor of Applied Statistics, University of Oxford
1 South Parks Road, Oxford OX1 3TG, UK



[Rd] ggplot2/plyr interaction with latest R-devel?
Prof Brian Ripley