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

[Martin Maechler]

>>>>> 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);

?



    > 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