[Rd] log(i, base=i) not giving 1

[Ben Bolker]

On 14-09-02 08:48 AM, Martin Maechler wrote:
>>>>>> peter dalgaard <pdalgd at gmail.com>
>>>>>>     on Tue, 2 Sep 2014 13:43:21 +0200 writes:
> 
>     > Impressive. Never ceases to amaze me what computers can do these days. ;-)
> 
> Indeed,
> 
>     > 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....
> 
> Amazingly indeed, it does:  From the few platforms I can try
> here, I only see the problem 
> on 32 bit Linux, both an (old) ubuntu 12.04.5 and Fedora 19.
> 
>> i <- 2:99; i[log(i, base=i) != 1]
> [1]  5  8 14 18 19 25 58 60 64 65 66 67 75 80 84 86 91
> 
> so '8' is not so special (as Ben suggested) and also not the
> only power of two for which this happens :
> 
>> i <- 2^(1:20); i[log(i, base=i) != 1]
> [1] 8    64   128  4096  8192 16384
> 
> Interestingly, it does not happen on 32-bit Windows, nor on any
> 64-bit version of R I've tried.
> Yes, it is amazing that with computer arithmetic, we can't
> even guarantee that   U / U = 1  exactly.
> 
> Is it basically free to check for x == base in there, so we
> could change to
> 
>  return (x == base) ? 1. : R_log(x) / R_log(base);
> 
> ?

  Should I submit a formal bug report about this, or should I assume it
has now been sufficiently noted by R-core?

  Opinions about whether it is better to fix in logbase (according to
Martin's suggestion) or in version.R (e.g. BDR's suggestion of using
log2(x)/3 or some other way that makes the function less sensitive) or both?

  Ben Bolker


> 
> 
> 
>     > 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
> 
>     > -- 
>     > Peter Dalgaard, Professor,
>     > Center for Statistics, Copenhagen Business School
>     > Solbjerg Plads 3, 2000 Frederiksberg, Denmark
>     > Phone: (+45)38153501
>     > Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com
> 
>     > ______________________________________________
>     > R-devel at r-project.org mailing list
>     > https://stat.ethz.ch/mailman/listinfo/r-devel
>