https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65723

            Bug ID: 65723
           Summary: uniform_real_distribution isn't uniform.
           Product: gcc
           Version: unknown
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: libstdc++
          Assignee: unassigned at gcc dot gnu.org
          Reporter: art at blahonga dot org

uniform_real_distribution according to all documentation I can find, including
a c++ standard draft is supposed to return uniformly distributed numbers.

Here's a simple test program (yes, I know it's not correct c++, but gcc seems
to have no problems with hexadecimal floating point constants which aren't
allowed in c++):

$ cat urd.cxx
#include <random>
#include <iostream>

#include <stdio.h>
#include <assert.h>

int
main()
{
        std::default_random_engine gen;
        double from = 0x1p52, to = from + 2;
        std::uniform_real_distribution<double> dis(from, to);
        long long bucket[3] = { 0 };
        int i;

        for (i = 0; i < 1000000; i++) {
                double r = dis(gen);
                unsigned int b = (int)(r - from);
                if (b >= 3) {
                        printf("foo: %d %f\n", b, r);
                }
                assert(b < 3);
                bucket[b]++;
        }

        for (i = 0; i < 3; i++) {
                printf("%lld\n", bucket[i]);
        }


        return 0;
}
$ c++ -v
[... snip ...]
gcc version 4.9.2 (Ubuntu 4.9.2-0ubuntu1~14.04)
$ c++ -std=c++11 -o foo urd.cxx && ./foo
250862
498838
250300

Good news is that gcc is in good company, both llvm libc++ and visual studio
(whatever crazy thing it links with) have the exact same bug. The industry
standard of generating random floating point numbers seems to be to pretend
that floating point behaves like math we learned in school and hope for the
best. 

I've written a post about this here:
http://www.blahonga.org/2015/04/state-of-art-of-random-double-generation.html
which is too long to include in a bug report. It also contains links to my
attempts to at least get parts of this problem solved.

Also, since I'm here, the range implemented by you (and llvm and windows) is
supposed to be "a <= x < b" according to the standard. Yet all three implement
"a <= x <= b".

Reply via email to