construction of a uniform double RNG from a random bits RNG

Discussion in 'C Programming' started by Francois Grieu, Apr 7, 2009.

  1. 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, Apr 7, 2009
    #1
    1. Advertising

  2. 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
    Francois Grieu, Apr 7, 2009
    #2
    1. Advertising

  3. 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
    Francois Grieu, Apr 7, 2009
    #3
  4. Francois Grieu <> 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.
    Ben Bacarisse, Apr 7, 2009
    #4
  5. Ben Bacarisse a écrit :
    > Francois Grieu <> 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
    Francois Grieu, Apr 7, 2009
    #5
  6. Francois Grieu

    Guest

    On Apr 7, 5:41 am, 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".



    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.
    , Apr 7, 2009
    #6
  7. wrote :
    > On Apr 7, 5:41 am, 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".

    >
    >
    > 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
    Francois Grieu, Apr 7, 2009
    #7
  8. Francois Grieu <> writes:

    > Ben Bacarisse a écrit :
    >> Francois Grieu <> 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.
    Ben Bacarisse, Apr 8, 2009
    #8
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Sydex
    Replies:
    12
    Views:
    6,452
    Victor Bazarov
    Feb 17, 2005
  2. aegis

    uniform random distribution

    aegis, Jan 30, 2005, in forum: C Programming
    Replies:
    7
    Views:
    354
    Julian V. Noble
    Jan 31, 2005
  3. Replies:
    8
    Views:
    431
  4. Adem24
    Replies:
    19
    Views:
    721
    Phil Carmody
    Jun 14, 2008
  5. Adem24
    Replies:
    18
    Views:
    663
    Phil Carmody
    Jun 14, 2008
Loading...

Share This Page