Re: Random number generation quirk

[Philipp Kutin]

Thinking about this some more. First, another correction:

On Sat, Mar 1, 2014 at 11:38 AM, Philipp Kutin <address@hidden> wrote:
> With epsf := eps('single') being 2^-23, the predecessor of 1 is
> 1-2^-23.

That should be "the predecessor of 1 is 1 - 2^-24".

If you want to do away with the coarseness, you may think along the
lines of, well, let's just do the intermediate calculation in double
precision, and only at the end demote to float. We repeat this as long
as this result equals one of the boundary values. Here's what you
might come up with then:

float randu32(void)  // (0, 1)
{
    static const double C = 1.0 / (65536.0 * 65536.0);  // 2^-32
    float f;

    do {
        f = C * (double)randi32();
    } while (f == 0.0 || f == 1.0);

    return f;
}

As an optimization, since f == 0.0 IFF randi32() returns 0, we can do
the f == 0.0 comparison in integer beforehand. How does this function
compare with the one posted earlier? First, it's clear that it samples
more values (without telling which ones these are). However, comparing
the converted result is slower than doing it on the integer early.
These are some timings for a whole 0 .. 2^32-1 cycle:

first version: 12.8 s
above version: 23.7 s
above version with early i==0 check: 20.7 s

But besides these efficiency issues, there's another numeric one that
relates to the FP rounding in the case of ties.

octave:37> for i=0:10, fprintf('%ld %ld\n', 2^24+i, single(2^24+i)); end
16777216 16777216
16777217 16777216  % round down
16777218 16777218
16777219 16777220  % round up
16777220 16777220
16777221 16777220  % round down
16777222 16777222
16777223 16777224  % round up
16777224 16777224
16777225 16777224  % round down
16777226 16777226

Because of this rounding definition, we don't sample the FP numbers in
an interval [2^(-n-1) 2^-n) (n >= 0) uniformly if, denoting the
difference between two consecutive numbers in this interval as d(n),
numbers "finer" than d(n) are used as intermediate values.

So much for "fun with floating point" for now. I'm attaching a C
source file with the three variants of randu32(), plus some simple
checking in main().

--Philipp

randfloat.c
Description: Text Data