More on random numbers... (Proposed correction on the FAQ)

Discussion in 'C Programming' started by Army1987, Aug 4, 2007.

  1. Army1987

    Army1987 Guest

    In the FAQ I read:
    If you just need to do something with probability 1/N, you could use

    if(rand() < (RAND_MAX+1u) / N)

    This has an obvious bias problem which is easily fixed (not
    perfectly, but an improvement).

    There are RAND_MAX + 1 values RAND_MAX can return. So we want to
    pick (RAND_MAX + 1.0) / N of them. This is probably not an
    integer. So we might want to round it to nearest integer. If we
    are willing to use floating point, it turns out that if x is a
    positive real number, there are ceil(x) non-negative integers such
    as n < x. (For example, if x is 3 there are 3 such values: 0, 1,
    and 2. If x is 3.2 there are 4 such values: 0, 1, 2 and 3.) Now
    the nearest integer to x is ceil(x - 0.5). [1] So we might try to
    do
    if (rand() < (RAND_MAX + 1.0) / N - 0.5)
    or, for an arbitrary fraction p, not necessarily 1 / N:
    if (rand() < p * (RAND_MAX + 1.0) - 0.5)
    Now, if we don't want to use floating point, getting back to the
    case 1 / N, we can use
    if (rand() < (RAND_MAX + 1u + N / 2) / N)

    [1] The problem is when p * (RAND_MAX + 1.0) is exactly half an
    odd integer. When p is 1 / N this isn't a great concern, because,
    where RAND_MAX + 1.0 is a power of two, this is impossible, unless
    N is 2 * RAND_MAX, which can be handled in a more obvious way as
    if (rand() == 0 && rand() > RAND_MAX / 2).

    --
    Army1987 (Replace "NOSPAM" with "email")
    "Never attribute to malice that which can be adequately explained
    by stupidity." -- R. J. Hanlon (?)
     
    Army1987, Aug 4, 2007
    #1
    1. Advertisements

  2. Army1987

    Army1987 Guest

    On Sat, 04 Aug 2007 11:57:46 +0200, Army1987 wrote:

    > In the FAQ I read:
    > If you just need to do something with probability 1/N, you could use
    >
    > if(rand() < (RAND_MAX+1u) / N)
    >
    > This has an obvious bias problem which is easily fixed (not
    > perfectly, but an improvement).
    >
    > There are RAND_MAX + 1 values RAND_MAX can return. So we want to

    Ok. I meant ...values rand() can return.

    --
    Army1987 (Replace "NOSPAM" with "email")
    "Never attribute to malice that which can be adequately explained
    by stupidity." -- R. J. Hanlon (?)
     
    Army1987, Aug 4, 2007
    #2
    1. Advertisements

  3. Army1987

    Eric Sosman Guest

    Army1987 wrote:
    > In the FAQ I read:
    > If you just need to do something with probability 1/N, you could use
    >
    > if(rand() < (RAND_MAX+1u) / N)


    In fairness to the FAQ, it goes on to explain that this
    technique is only suitable for N much smaller than RAND_MAX.

    > This has an obvious bias problem which is easily fixed (not
    > perfectly, but an improvement).
    > [...]
    > if (rand() < (RAND_MAX + 1u + N / 2) / N)


    That is, rounding the quotient instead of truncating it.
    Yes, that could give a slightly more accurate result. How
    much more accurate? By less than one part in (RAND_MAX+1u),
    which could be as much as 0.0000305+ for the smallest legal
    RAND_MAX.

    But rounding can still leave you with an error of up to
    half that amount! If N is large enough that the truncation
    error of 0.0031% is significant, the rounding error of 0.0015%
    is probably *also* significant. Or to put it another way, if
    N is large enough to make rounding attractive, it is also
    large enough to make rounding ineffective; the improvement is
    only of interest when it's not good enough. For N that large,
    I think you should be using a rejection technique along the
    lines of the one illustrated in the FAQ.

    --
    Eric Sosman
    lid
     
    Eric Sosman, Aug 4, 2007
    #3
  4. In article <>,
    Army1987 <> wrote:

    >There are RAND_MAX + 1 values RAND_MAX can return.


    RAND_MAX is the maximum value that can be returned, but there is
    no certainty that every value from 0 to RAND_MAX is returnable.
    In particular, 0 was historically not returnable in some implementations.

    --
    There are some ideas so wrong that only a very intelligent person
    could believe in them. -- George Orwell
     
    Walter Roberson, Aug 4, 2007
    #4
  5. Army1987

    Army1987 Guest

    On Sat, 04 Aug 2007 18:36:56 +0000, Walter Roberson wrote:
    > RAND_MAX is the maximum value that can be returned, but there is
    > no certainty that every value from 0 to RAND_MAX is returnable.
    > In particular, 0 was historically not returnable in some implementations.


    Is there any way this can be known at compile-time or preprocessing
    time? I'd consider such an implementation to be broken [1], and I
    think it had better have rand() return _internal_rand() - 1 and
    #define RAND_MAX (_INTERNAL_RAND_MAX - 1).

    [1] An objection that could be made is that the Standard says that
    the range is 0 to RAND_MAX, but it doesn't say how bad they
    can be, so nothing stops an implementation never returning 0 to
    be nonconforming. But then, you could argue that not even an
    implementation never returning 1 is nonconforming, not even one
    never returning 2, and so on, until you could claim that
    #define RAND_MAX 32767
    int rand(void) { return RAND_MAX; }
    void srand(unsigned int __seed) { return; }
    is conforming.
    --
    Army1987 (Replace "NOSPAM" with "email")
    "Never attribute to malice that which can be adequately explained
    by stupidity." -- R. J. Hanlon (?)
     
    Army1987, Aug 5, 2007
    #5
    1. Advertisements

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. Lars-Erik Aabech
    Replies:
    8
    Views:
    1,019
    Lars-Erik Aabech
    Apr 28, 2005
  2. globalrev
    Replies:
    4
    Views:
    1,079
    Gabriel Genellina
    Apr 20, 2008
  3. Alex Untitled
    Replies:
    11
    Views:
    894
    Giampiero Zanchi
    Nov 16, 2009
  4. PerlFAQ Server

    FAQ 4.10 Why aren't my random numbers random?

    PerlFAQ Server, Feb 12, 2011, in forum: Perl Misc
    Replies:
    0
    Views:
    391
    PerlFAQ Server
    Feb 12, 2011
  5. PerlFAQ Server

    FAQ 4.10 Why aren't my random numbers random?

    PerlFAQ Server, Apr 27, 2011, in forum: Perl Misc
    Replies:
    0
    Views:
    400
    PerlFAQ Server
    Apr 27, 2011
  6. Richard Cornford

    FAQ update and proposed updates

    Richard Cornford, Nov 7, 2005, in forum: Javascript
    Replies:
    2
    Views:
    242
    Dr John Stockton
    Nov 9, 2005
  7. Randy Webb

    FAQ Proposed Anchor Names

    Randy Webb, Dec 19, 2006, in forum: Javascript
    Replies:
    5
    Views:
    304
    Randy Webb
    Dec 20, 2006
  8. VK
    Replies:
    15
    Views:
    1,919
    Dr J R Stockton
    May 2, 2010
Loading...