construction of a uniform double RNG from a random bits RNG

[Francois Grieu]

Hi,

suppose I have an unsigned 32-bit type tu32 and a robust
32-bit Random Number Generator
tu32 rng32(void);
producing output uniformly distributed on [0 .. 0xFFFFFFFF].

I want to build a robust double RNG with output uniformly
distributed on ]0..1[ (or similar but known, e.g.
[DBL_MIN..1-DBL_EPSILON] ), for some sound definition of
uniformly distributed matching what double allows.
double rngdouble(void);
In particular, I do care that even a subsample of the output
restricted to some interval like [0 .. 1e-7] appear uniform
rather than "grainy".

Is there a classic good technique?

Francois Grieu

 

[Francois Grieu]

pete a écrit :

Hi,

suppose I have an unsigned 32-bit type tu32 and a robust
32-bit Random Number Generator
tu32 rng32(void);
producing output uniformly distributed on [0 .. 0xFFFFFFFF].

I want to build a robust double RNG with output uniformly
distributed on ]0..1[ (or similar but known, e.g.
[DBL_MIN..1-DBL_EPSILON] ), for some sound definition of
uniformly distributed matching what double allows.
double rngdouble(void);
In particular, I do care that even a subsample of the output
restricted to some interval like [0 .. 1e-7] appear uniform
rather than "grainy".

Is there a classic good technique?

http://c-faq.com/lib/rand48.html

From this I infer:

#define NEEDEDBITS 96
double rngdouble(void)
{
double x = 0;
double denom = 1;
int b;

tu32 r;
while (0==(r = rng32()))
denom *= 4294967296.; // 2^32
for(b = NEEDEDBITS; b>0; b -= 32)
{
denom *= 4294967296.; // 2^32
x += rng32()/denom;
}
return x;
}

I'm uncertain on the effect of rounding; does it introduce
a bias? Can 0 and 1 be reached ?

François Grieu

 

[Francois Grieu]

pete a écrit :

Hi,

suppose I have an unsigned 32-bit type tu32 and a robust
32-bit Random Number Generator
tu32 rng32(void);
producing output uniformly distributed on [0 .. 0xFFFFFFFF].

I want to build a robust double RNG with output uniformly
distributed on ]0..1[ (or similar but known, e.g.
[DBL_MIN..1-DBL_EPSILON] ), for some sound definition of
uniformly distributed matching what double allows.
double rngdouble(void);
In particular, I do care that even a subsample of the output
restricted to some interval like [0 .. 1e-7] appear uniform
rather than "grainy".

Is there a classic good technique?

http://c-faq.com/lib/rand48.html

From this I infer (hopefully corrected, and simplified)

#define NEEDEDBITS 96
double rngdouble(void)
{
double x = 0;
double denom = 1;
int b;
tu32 r;
for(b = 0; b<NEEDEDBITS
{
x += (r = rng32())/(denom *= 4294967296.);
if (b||r)
b += 32;
}
return x;
}

I'm uncertain on the effect of rounding; does it introduce
a bias? Can 0 and 1 be reached ?

François Grieu

 

[Ben Bacarisse]

suppose I have an unsigned 32-bit type tu32 and a robust
32-bit Random Number Generator
tu32 rng32(void);
producing output uniformly distributed on [0 .. 0xFFFFFFFF].

I want to build a robust double RNG with output uniformly
distributed on ]0..1[ (or similar but known, e.g.
[DBL_MIN..1-DBL_EPSILON] ), for some sound definition of
uniformly distributed matching what double allows.
double rngdouble(void);
In particular, I do care that even a subsample of the output
restricted to some interval like [0 .. 1e-7] appear uniform
rather than "grainy".

Is there a classic good technique?

If you don't want 0 or 1, and you are not too fussy about exactly how
close you get (provided the bounds are known) then I would just return

(rng32() + 1.0)/(0xffffffff + 2.0)

or

rng32()/(0xffffffff + 1.0) + DBL_EPSILON

There will be expressions that get closer to bounds ideal range or
(0..1) but you did not seem to be too concerned about that.

 

[Francois Grieu]

Ben Bacarisse a écrit :

suppose I have an unsigned 32-bit type tu32 and a robust
32-bit Random Number Generator
tu32 rng32(void);
producing output uniformly distributed on [0 .. 0xFFFFFFFF].

I want to build a robust double RNG with output uniformly
distributed on ]0..1[ (or similar but known, e.g.
[DBL_MIN..1-DBL_EPSILON] ), for some sound definition of
uniformly distributed matching what double allows.
double rngdouble(void);
In particular, I do care that even a subsample of the output
restricted to some interval like [0 .. 1e-7] appear uniform
rather than "grainy".

Is there a classic good technique?

If you don't want 0 or 1, and you are not too fussy about exactly how
close you get (provided the bounds are known) then I would just return

(rng32() + 1.0)/(0xffffffff + 2.0)

or

rng32()/(0xffffffff + 1.0) + DBL_EPSILON

Fine as far as the range goes, but these do not match my criteria:
when you consider the subset of the output which is no more than
1e-7, it is very sparse (few hundred values). A mere density plot
of -log(rngdouble()) for a few billion outputs makes it visible.

Francois Grieu

 

[robertwessel2]

Hi,

suppose I have an unsigned 32-bit type  tu32  and a robust
32-bit Random Number Generator
   tu32 rng32(void);
producing output uniformly distributed on [0 .. 0xFFFFFFFF].

I want to build a robust double RNG with output uniformly
distributed on ]0..1[ (or similar but known, e.g.
[DBL_MIN..1-DBL_EPSILON] ), for some sound definition of
uniformly distributed matching what  double  allows.
   double rngdouble(void);
In particular, I do care that even a subsample of the output
restricted to some interval like [0 .. 1e-7] appear uniform
rather than "grainy".

I want to point out that the values you can actually put in a double
are not uniformly distributed, so you need to clarify exactly what you
mean by that. But let's say you want to generate numbers from
[0..1,000,000) at .001 steps (acknowledging that this will not
generate exact floating point values for the usual reasons). Generate
a good quality random integer of at least 30 bits, discard any values
larger than or equal to the biggest multiple of 1,000,000,000 smaller
than the maximum random number your source will generate (for example
if your source produces a 32 bit unsigned, discard any values larger
than 3,999,999,999 (and if you discard, you need to generate a new
value). Mod that by 1,000,000,000, and scale the result to your
output range.

 

[Francois Grieu]

Hi,

suppose I have an unsigned 32-bit type tu32 and a robust
32-bit Random Number Generator
tu32 rng32(void);
producing output uniformly distributed on [0 .. 0xFFFFFFFF].

I want to build a robust double RNG with output uniformly
distributed on ]0..1[ (or similar but known, e.g.
[DBL_MIN..1-DBL_EPSILON] ), for some sound definition of
uniformly distributed matching what double allows.
double rngdouble(void);
In particular, I do care that even a subsample of the output
restricted to some interval like [0 .. 1e-7] appear uniform
rather than "grainy".

I want to point out that the values you can actually put in a double
are not uniformly distributed, so you need to clarify exactly what you
mean by that.

I simply want the result to behave as close to a uniform random
variable over the interval 0..1 as allowed by double; or nearly as
close.

But let's say you want to generate numbers from
[0..1,000,000) at .001 steps (acknowledging that this will not
generate exact floating point values for the usual reasons).

Rather, say that I am interested in simulating the quantity
-log(v) where v is a uniform random variable over 0..1
[This models in particular the delay between events occurring
randomly if there is one event per time unit on average]

> Generate good quality random integer of at least 30 bits, discard
> any values larger than or equal to the biggest multiple of
> 1,000,000,000 smaller than the maximum random number your source
> will generate (for example if your source produces a 32 bit unsigned,
> discard any values larger than 3,999,999,999 (and if you discard,
> you need to generate a new value). Mod that by 1,000,000,000,
> and scale the result to your output range.

If I make an rngdouble() according to this principle and plots
the value -log(rngdouble()) some 10e11 times on a linear scale from
0 to 25, there will be a distinctive pattern on the upper end,
which is due only to a defect of the implementation of rngdouble(),
and not at all on a limitation of double.

Francois Grieu

 

[Ben Bacarisse]

Ben Bacarisse a écrit :

suppose I have an unsigned 32-bit type tu32 and a robust
32-bit Random Number Generator
tu32 rng32(void);
producing output uniformly distributed on [0 .. 0xFFFFFFFF].

I want to build a robust double RNG with output uniformly
distributed on ]0..1[ (or similar but known, e.g.
[DBL_MIN..1-DBL_EPSILON] ), for some sound definition of
uniformly distributed matching what double allows.
double rngdouble(void);
In particular, I do care that even a subsample of the output
restricted to some interval like [0 .. 1e-7] appear uniform
rather than "grainy".

Is there a classic good technique?

If you don't want 0 or 1, and you are not too fussy about exactly how
close you get (provided the bounds are known) then I would just return

(rng32() + 1.0)/(0xffffffff + 2.0)

or

rng32()/(0xffffffff + 1.0) + DBL_EPSILON

Fine as far as the range goes, but these do not match my criteria:
when you consider the subset of the output which is no more than
1e-7, it is very sparse (few hundred values). A mere density plot
of -log(rngdouble()) for a few billion outputs makes it visible.

With 32 bits, the above is about as good as it gets so you have to
decide what you want to sacrifice. If you are prepared to compromise
period and independence by using two calls to your PRNG to build a
more dense floating point number, this may work for you (if you have a
64 bit integer type).

unsigned long long rng64(void)
{
return ((unsigned long long)rng32() << 32) + rng32();
}

double drng(void)
{
return rng64() / 18446744073709551616.0 + DBL_EPSILON;
}

This halves the effective period of your generator and there is a
small chance you might get odd effects coming from the fact that the
two numbers used to make the double are algorithmically related,
albeit in a complex way.

Alternatively, if you know you are going to be sampling in a narrow
range, use the 32 bits to get the smoothest double in that narrower
range (but you will have thought of that already).