[Newbies] Re: What's wrong with this statement?

[nicolas cellier]

 <johnps11 <at> bigpond.com> writes:

> 
> The issue is almost certainly in atRandom:
> 
> Running
> 
> 30  timesRepeat: [ | a |
>     a:= (2 raisedTo: 57) atRandom.
>     Transcript show: a ; show:  ' ' ; show: a even ; cr ].
> 
> Gives me 30 falses.
> 
> The odds of that are less than one in a billion (assuming uniformly
> distributed integers).
> 
> This doesn't happen for 2^56 or lower, and does for all 2^n, n>= 57.
> 
> My guess is that it's either a hardware issue or something about Random. 
> I'm not clever enough to understand PRNGs, so I'll leave it others to work
> out what the answer is, although I suspect that for more than 56 bits you
> need a PRNG that uses LargeIntegers  and not Floats.
> 

Nothing surprising once again.
Float have a 53 bits integer mantissa m and and a biased exponent e:
    f := m * (2 raisedTo: e).
    m highBit = 53.
The PRNG has a uniform distribution on [0,1).
f < 1 means 100% of floats have e < -53.
Let us divide the uniform intervals in two recursively:

50% have e = -54 (0.5<=f<1)
50% have e < -54 (f<0.5)

25% have e = -55 (0.25<=f<0.5)
25% have e < -55 (f<0.25)
75% have e >=-55 (f>=0.25)

12.5% have e = -56 (0.125<=f<0.25)
12.5% have e < -56 (f<0.125)
87.5% have e >=-56 (f>=0.125)

6.25%  have e = -57 (0.0625<=f<0.125)
6.25%  have e < -57 (f<0.0625)
92.75% have e >=-57 (f>=0.0625)

3.125%  have e = -58 (0.03125<=f<0.0625)
3.125%  have e < -58 (f<0.03125)
96.875% have e >=-58 (f>=0.03125)

So when you multiply f by (2 raisedTo: 57),
that means (m bitShift: e+57),
87.5% have (e+57>0),
thus are shifted at least one bit to the left,
and you get 87.5% of even numbers.
since algorithm truncate then add 1, you then get 87.5% of odd results

I suspect that due to division by 16r7FFFFFFF, last bit of mantissa m is often
zero. In which case, you get 92.75% of odd results.

Or maybe last two bits, in which case you have 96.875% of odd results.

   ((1 to: 10000) count: [:i | (2 raisedTo: 57) atRandom odd])  / 10000.0 .
-> 0.9692
Yes, probably last two bits of these floats are set to 0...

The limit is by no way (2 raisedTo: 57).

The true limit at which you won't get any even number ?
well, smallest generated pseudo random float is 1/(2^31-1)

A little higher than 1.0 timesTwoPoser: -31.
   smallest := ((1.0 timesTwoPower: 31) - 1) reciprocal.
   smallest significandAsInteger printStringBase: 2.
-> '10000000000000000000000000000001000000000000000000000'

The biased exponent is then (-31-52=-83).
   (smallest significandAsInteger asFloat timesTwoPower: smallest exponent-52)
    = smallest
-> true

So, (2 raisedTo: 84) atRandom is guaranteed to be odd.

We conjectured that last two bits are zero,
maybe (2 raisedTo: 82) atRandom is already 100% odd?
   (1 to: 10000000 by: 11) detect: [:e | ((e * smallest) significandAsInteger
     bitAnd: 2r11) > 0] ifNone: [nil]
->1048587
Hmm the conjecture is false...

Floats are tricky, but once you understand they are just rational with
denominator being a power of two and inexact operations truncating the mantissa
at 53 bits, the mystery vanishes a bit.

Nicolas