On 02/09/2014 22:43, Ben Bolker wrote:
On 14-09-02 08:48 AM, Martin Maechler wrote:
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);
?
Should I submit a formal bug report about this, or should I assume it
has now been sufficiently noted by R-core?
It is under consideration by the author (not me).
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?
Yet again, it is better to write code which does not make false
assumptions about machine arithmetic. Fixing one very rare case at the
expense of speed for all others is not a good idea.
Studying http://www.validlab.com/goldberg/addendum.html is a good idea
for those unfamiliar with the pitfalls of i86 FPUs. It has an example
essentially the same as this one right at the top.
Also, Ben Bolker needs to change his compiler options to ones better
than the defaults in his unstated Linux distro. Why anyone is using
ix86 Linux nowadays is hard to understand when x86_64 is (on all but the
very smallest machines) faster and more reliable.
Ben Bolker
> 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
--
Brian D. Ripley, rip...@stats.ox.ac.uk
Emeritus Professor of Applied Statistics, University of Oxford
1 South Parks Road, Oxford OX1 3TG, UK
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
