Velocity Reviews - Computer Hardware Reviews

Velocity Reviews > Newsgroups > Programming > C Programming > construction of a uniform double RNG from a random bits RNG

Reply
Thread Tools

construction of a uniform double RNG from a random bits RNG

 
 
Francois Grieu
Guest
Posts: n/a
 
      04-07-2009
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
 
Reply With Quote
 
 
 
 
Francois Grieu
Guest
Posts: n/a
 
      04-07-2009
pete a écrit :
> Francois Grieu wrote:
>> 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
 
Reply With Quote
 
 
 
 
Francois Grieu
Guest
Posts: n/a
 
      04-07-2009
pete a écrit :
> Francois Grieu wrote:
>> 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
 
Reply With Quote
 
Ben Bacarisse
Guest
Posts: n/a
 
      04-07-2009
Francois Grieu <(E-Mail Removed)> writes:

> 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.

--
Ben.
 
Reply With Quote
 
Francois Grieu
Guest
Posts: n/a
 
      04-07-2009
Ben Bacarisse a écrit :
> Francois Grieu <(E-Mail Removed)> writes:
>
>> 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
 
Reply With Quote
 
robertwessel2@yahoo.com
Guest
Posts: n/a
 
      04-07-2009
On Apr 7, 5:41*am, Francois Grieu <(E-Mail Removed)> wrote:
> 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.
 
Reply With Quote
 
Francois Grieu
Guest
Posts: n/a
 
      04-07-2009
http://www.velocityreviews.com/forums/(E-Mail Removed) wrote :
> On Apr 7, 5:41 am, Francois Grieu <(E-Mail Removed)> wrote:
>> 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
 
Reply With Quote
 
Ben Bacarisse
Guest
Posts: n/a
 
      04-07-2009
Francois Grieu <(E-Mail Removed)> writes:

> Ben Bacarisse a écrit :
>> Francois Grieu <(E-Mail Removed)> writes:
>>
>>> 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).

--
Ben.
 
Reply With Quote
 
 
 
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are Off


Similar Threads
Thread Thread Starter Forum Replies Last Post
Need a good RNG and a LCG, both with a max period >= 31 bits Adem24 C++ 19 06-14-2008 08:57 PM
Need a good RNG and a LCG, both with a max period >= 31 bits Adem24 C Programming 18 06-14-2008 08:57 PM
Default construction versus construction with initial values Ook C++ 10 10-08-2005 09:00 PM
cannot convert parameter from 'double (double)' to 'double (__cdecl *)(double)' error Sydex C++ 12 02-17-2005 06:30 PM
8-Bits vs 12 or 16 bits/pixel; When does more than 8 bits count ? Al Dykes Digital Photography 3 12-29-2003 07:08 PM



Advertisments