>>>>> peter dalgaard <pda...@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);
?
> On 02 Sep 2014, at 03:27 , Ben Bolker <bbol...@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....@cbs.dk Priv: pda...@gmail.com
> ______________________________________________
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
