log2 buggy or is a real thing?

[Don Clugston]

On 04/04/12 13:40, bearophile wrote:
> Jonathan M Davis:
>
>> This progam:
>>
>> import std.math;
>> import std.stdio;
>> import std.typetuple;
>>
>> ulong log2(ulong n)
>> {
>>      return n == 1 ? 0
>>                    : 1 + log2(n / 2);
>> }
>>
>> void print(ulong value)
>> {
>>      writefln("%s: %s %s", value, log2(value), std.math.log2(value));
>> }
>>
>> void main()
>> {
>>      foreach(T; TypeTuple!(byte, ubyte, short, ushort, int, uint, long, ulong))
>>          print(T.max);
>> }
>>
>>
>> prints out
>>
>> 127: 6 6.98868
>> 255: 7 7.99435
>> 32767: 14 15
>> 65535: 15 16
>> 2147483647: 30 31
>> 4294967295: 31 32
>> 9223372036854775807: 62 63
>> 18446744073709551615: 63 64
>>
>>
>> So, the question is: Is std.math.log2 buggy, or is it just an issue with the
>> fact that std.math.log2 is using reals, or am I completely misunderstanding
>> something here?
>
> What is the problem you are perceiving?
>
> The values you see are badly rounded, this is why I said that the Python floating point printing function is better than the D one :-)
>
> I print the third result with more decimal digits using %1.20f, to avoid that bad rounding:
>
>
> import std.stdio, std.math, std.typetuple;
>
> ulong mlog2(in ulong n) pure nothrow {
>      return (n == 1) ? 0 : (1 + mlog2(n / 2));
> }
>
> void print(in ulong x) {
>      writefln("%s: %s %1.20f", x, mlog2(x), std.math.log2(x));
> }
>
> void main() {
>      foreach (T; TypeTuple!(byte, ubyte, short, ushort, int, uint, long, ulong))
>          print(T.max);
>      print(9223372036854775808UL);
>      real r = 9223372036854775808UL;
>      writefln("%1.20f", r);
> }
>
>
> The output is:
>
> 127: 6 6.98868468677216585300
> 255: 7 7.99435343685885793770
> 32767: 14 14.99995597176955822200
> 65535: 15 15.99997798605273604400
> 2147483647: 30 30.99999999932819277100
> 4294967295: 31 31.99999999966409638400
> 9223372036854775807: 62 63.00000000000000000100
> 18446744073709551615: 63 64.00000000000000000000
> 9223372036854775808: 63 63.00000000000000000100
> 9223372036854775807.80000000000000000000
>
> And it seems essentially correct, the only significant (but very little) problem I am seeing is
> log2(9223372036854775807) that returns a value>  63, while of course in truth it's strictly less than 63 (found with Wolfram Alpha):
> 62.999999999999999999843582690243412222355489972678233
> The cause of that small error is the limited integer representation precision of reals.

I don't think so. For 80-bit reals, every long can be represented 
exactly in an 80 bit real, as can every ulong from 0 up to and including 
ulong.max - 1. The only non-representable built-in integer is ulong.max, 
which (depending on rounding mode) gets rounded up to ulong.max+1.

The decimal floating point output for DMC, which DMD uses, is not very 
accurate. I suspect the value is actually <= 63.

> In Python there is a routine to compute precisely the logarithm of large integer numbers. It will be useful to have in D too, for BigInts too.
>
> Bye,
> bearophile