From bbolker at gmail.com  Tue Sep  2 23:43:17 2014
From: bbolker at gmail.com (Ben Bolker)
Date: Tue, 2 Sep 2014 17:43:17 -0400
Subject: [Rd] log(i, base=i) not giving 1
In-Reply-To: <21509.48311.450620.998931@stat.math.ethz.ch>
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>
	<21509.48311.450620.998931@stat.math.ethz.ch>
Message-ID: <540639F5.2080002@gmail.com>

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
>



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