[Bug libstdc++/33815] New: tr1::uniform_int isn't uniform

[Bug libstdc++/33815] New: tr1::uniform_int isn't uniform

[jkherciueh at gmx dot net]

Consider:

#include <ctime>
#include <iostream>
#include <vector>
#include <tr1/random>

void test ( unsigned int num_buckets, unsigned int seed = 3141592 ) {
  std::cout << "using " << num_buckets << " buckets:\n";
  std::tr1::mt19937 eng(seed);
  const unsigned long range = eng.max() / 5.5;
  std::tr1::uniform_int<unsigned long> rnd( 0, range );
  std::vector<unsigned long> bucket ( num_buckets, 0 );
  for ( unsigned int i = 0; i < 1000000; ++i ) {
    ++bucket[ rnd(eng) / (range / num_buckets) ];
  }
  for ( unsigned int i = 0; i < num_buckets; ++i ) {
    std::cout << i << ": " << bucket[i] << '\n';
  }
  std::cout << '\n';
}

int main ( void ) {
  test( 2 );
  test( 3 );
  test( 4 );
  test( 6 );
  test( 10 );
  test( 20 );
  test( 30 );
  test( 40 );
}

It yields:

using 2 buckets.
0: 545389
1: 454611

using 3 buckets.
0: 363350
1: 333829
2: 302821

using 4 buckets.
0: 272507
1: 272882
2: 227489
3: 227122

using 6 buckets.
0: 181743
1: 181607
2: 182039
3: 151790
4: 151182
5: 151639

...


This doesn't look uniform. It looks more like the buckets in the
lower half get hit significantly more often.

Credit: This strange behavior had been reported by David Benbennick
to Boost long ago and was apparently never fixed.


-- 
           Summary: tr1::uniform_int isn't uniform
           Product: gcc
           Version: 4.2.1
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: libstdc++
        AssignedTo: unassigned at gcc dot gnu dot org
        ReportedBy: jkherciueh at gmx dot net
 GCC build triplet: i686-pc-linux
  GCC host triplet: i686-pc-linux
GCC target triplet: i686-pc-linux


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33815

[Bug libstdc++/33815] tr1::uniform_int isn't uniform

[pcarlini at suse dot de]

------- Comment #1 from pcarlini at suse dot de  2007-10-19 09:37 -------
Note that our implementation has nothing to do with the Boost one. The problem
happens when the range is big. Anyway, we know our implementation is naive in
this case, we'll look into it.


-- 

pcarlini at suse dot de changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|UNCONFIRMED                 |NEW
     Ever Confirmed|0                           |1
   Last reconfirmed|0000-00-00 00:00:00         |2007-10-19 09:37:59
               date|                            |


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33815

[Bug libstdc++/33815] tr1::uniform_int isn't uniform

[pcarlini at suse dot de]

------- Comment #2 from pcarlini at suse dot de  2007-10-19 09:40 -------
Actually, the problem happens only for some specific values of the range, like
eng.max() / 5.5, dividing by 2 or 10 is ok. Indeed, seems a binary arithmetic
problem.


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33815

[Bug libstdc++/33815] tr1::uniform_int isn't uniform

[pcarlini at suse dot de]

-- 

pcarlini at suse dot de changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
         AssignedTo|unassigned at gcc dot gnu   |pcarlini at suse dot de
                   |dot org                     |
             Status|NEW                         |ASSIGNED


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33815

[Bug libstdc++/33815] tr1::uniform_int isn't uniform

[jkherciueh at gmx dot net]

------- Comment #3 from jkherciueh at gmx dot net  2007-10-19 10:13 -------
> Actually, the problem happens only for some specific values of the range, like
> eng.max() / 5.5, dividing by 2 or 10 is ok. Indeed, seems a binary arithmetic
> problem.

I conjecture that the problem happens if and only if the range size does not
divide evenly into eng.max(). In fact, the fractional part determines which
buckets get hit more often. If you do 4.3 instead of 5.5, the first 3 out of 10
or the first 6 ot out 20 buckets get the load. If this conjecture is correct,
it would be more accurate to say that the algorithm performs as promised only
for a small set of range sizes (the divisors of eng.max()) and is biased for
all other values.


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33815

[Bug libstdc++/33815] tr1::uniform_int isn't uniform

[pcarlini at suse dot de]

------- Comment #4 from pcarlini at suse dot de  2007-10-19 10:46 -------
I'm going to commit to mainline a stop-gap solution. If you could confirm, as I
believe, that it's at least an improvement, we can have it in 4_2-branch too.


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33815

[Bug libstdc++/33815] tr1::uniform_int isn't uniform

[pcarlini at suse dot de]

------- Comment #5 from pcarlini at suse dot de  2007-10-19 10:49 -------
To wit, this kind of change, you can certainly apply it to the 4.2.1 headers
as-is:

Index: random
===================================================================
--- random      (revision 129456)
+++ random      (working copy)
@@ -1607,7 +1607,8 @@
         {
          typedef typename __gnu_cxx::__add_unsigned<typename
            _UniformRandomNumberGenerator::result_type>::__type __utype;
-         return result_type(__utype(__urng()) % (__max - __min + 1)) + __min;
+         return result_type((__max - __min + 1.0L) * __utype(__urng())
+                            / (__urng.max() + 1.0L)) + __min;
        }

       template<typename _UniformRandomNumberGenerator>


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33815

[Bug libstdc++/33815] tr1::uniform_int isn't uniform

[pcarlini at suse dot de]

------- Comment #6 from pcarlini at suse dot de  2007-10-19 11:26 -------
Actually, would be:

Index: random
===================================================================
--- random      (revision 129456)
+++ random      (working copy)
@@ -1607,7 +1607,8 @@
         {
          typedef typename __gnu_cxx::__add_unsigned<typename
            _UniformRandomNumberGenerator::result_type>::__type __utype;
-         return result_type(__utype(__urng()) % (__max - __min + 1)) + __min;
+         return result_type((__max - __min + 1.0L) * __utype(__urng())
+                            / (__utype(__urng.max()) + 1.0L)) + __min;
        }

       template<typename _UniformRandomNumberGenerator>


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33815

[Bug libstdc++/33815] tr1::uniform_int isn't uniform

[paolo at gcc dot gnu dot org]

------- Comment #7 from paolo at gcc dot gnu dot org  2007-10-19 17:36 -------
Subject: Bug 33815

Author: paolo
Date: Fri Oct 19 17:36:03 2007
New Revision: 129493

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=129493
Log:
2007-10-19  Paolo Carlini  <pcarlini@suse.de>

        PR libstdc++/33815
        * include/tr1_impl/random
        (uniform_int<>::_M_call(_UniformRandomNumberGenerator&, result_type,
        result_type, true_type)): Avoid the modulo (which uses the low-order
        bits).

Modified:
    trunk/libstdc++-v3/ChangeLog
    trunk/libstdc++-v3/include/tr1_impl/random


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33815

[Bug libstdc++/33815] tr1::uniform_int isn't uniform

[paolo at gcc dot gnu dot org]

------- Comment #8 from paolo at gcc dot gnu dot org  2007-10-20 10:02 -------
Subject: Bug 33815

Author: paolo
Date: Sat Oct 20 10:02:34 2007
New Revision: 129507

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=129507
Log:
2007-10-20  Paolo Carlini  <pcarlini@suse.de>

        * include/tr1/random
        (uniform_int<>::_M_call(_UniformRandomNumberGenerator&, result_type,
        result_type, true_type)): Fix small thinko.

2007-10-19  Paolo Carlini  <pcarlini@suse.de>

        PR libstdc++/33815
        * include/tr1/random
        (uniform_int<>::_M_call(_UniformRandomNumberGenerator&, result_type,
        result_type, true_type)): Avoid the modulo (which uses the low-order
        bits).

Modified:
    branches/gcc-4_2-branch/libstdc++-v3/include/tr1/random


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33815

[Bug libstdc++/33815] tr1::uniform_int isn't uniform

[paolo at gcc dot gnu dot org]

------- Comment #9 from paolo at gcc dot gnu dot org  2007-10-20 10:03 -------
Subject: Bug 33815

Author: paolo
Date: Sat Oct 20 10:03:10 2007
New Revision: 129508

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=129508
Log:
2007-10-20  Paolo Carlini  <pcarlini@suse.de>

        * include/tr1/random
        (uniform_int<>::_M_call(_UniformRandomNumberGenerator&, result_type,
        result_type, true_type)): Fix small thinko.

2007-10-19  Paolo Carlini  <pcarlini@suse.de>

        PR libstdc++/33815
        * include/tr1/random
        (uniform_int<>::_M_call(_UniformRandomNumberGenerator&, result_type,
        result_type, true_type)): Avoid the modulo (which uses the low-order
        bits).

Modified:
    branches/gcc-4_2-branch/libstdc++-v3/ChangeLog


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33815

[Bug libstdc++/33815] tr1::uniform_int isn't uniform

[pcarlini at suse dot de]

------- Comment #10 from pcarlini at suse dot de  2007-10-20 10:05 -------
Fixed for 4.2.3.


-- 

pcarlini at suse dot de changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|ASSIGNED                    |RESOLVED
         Resolution|                            |FIXED
   Target Milestone|---                         |4.2.3


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33815

[Bug libstdc++/33815] tr1::uniform_int isn't uniform

[jkherciueh at gmx dot net]

------- Comment #11 from jkherciueh at gmx dot net  2007-10-20 11:07 -------
(In reply to comment #6)
> Actually, would be:
> 
> Index: random
> ===================================================================
> --- random      (revision 129456)
> +++ random      (working copy)
> @@ -1607,7 +1607,8 @@
>          {
>           typedef typename __gnu_cxx::__add_unsigned<typename
>             _UniformRandomNumberGenerator::result_type>::__type __utype;
> -         return result_type(__utype(__urng()) % (__max - __min + 1)) + __min;
> +         return result_type((__max - __min + 1.0L) * __utype(__urng())
> +                            / (__utype(__urng.max()) + 1.0L)) + __min;
>         }
> 
>        template<typename _UniformRandomNumberGenerator>
> 

I just tested this. It does a much better job at obfuscating the bias, which is
still there (after all you are rescaling linearly). If you rescale a domain of
size 16 to a range of size 12, some buckets will get 2 and some will get 1.

The following program demonstrates this:

#include <ctime>
#include <cassert>
#include <iostream>
#include <vector>
#include <tr1/random>

void test ( unsigned int num_buckets, unsigned int seed = 3141592 ) {
  std::cout << "using " << num_buckets << " buckets:\n";
  std::tr1::mt19937 eng(seed);
  // we make sure that the range size is a multiple of
  // the number of buckets:
  const unsigned long target_range = ( eng.max() + 1.0L ) * 0.75;
  const unsigned long bucket_size = target_range / num_buckets;
  const unsigned long range = bucket_size * num_buckets - 1;
  assert( bucket_size * num_buckets == range + 1 );
  std::tr1::uniform_int<unsigned long> rnd( 0, range );
  std::vector<unsigned long> bucket ( num_buckets, 0 );
  for ( unsigned int i = 0; i < 1000000; ++i ) {
    ++bucket[ rnd(eng) % num_buckets ];
    //++bucket[ rnd(eng) / bucket_size ];
  }
  for ( unsigned int i = 0; i < num_buckets; ++i ) {
    std::cout << i << ": " << bucket[i] << '\n';
  }
  std::cout << '\n';
}

int main ( void ) {
  test( 2 );
  test( 3 );
  test( 4 );
  test( 5 );
  test( 6 );
  test( 7 );
}

Output:

0: 500742
1: 499258

using 3 buckets:
0: 499275
1: 250467
2: 250258

using 4 buckets:
0: 250024
1: 249533
2: 250718
3: 249725

using 5 buckets:
0: 199340
1: 200365
2: 199559
3: 200129
4: 200607

using 6 buckets:
0: 249861
1: 124864
2: 125278
3: 249414
4: 125603
5: 124980

using 7 buckets:
0: 142866
1: 142908
2: 142945
3: 142964
4: 142671
5: 142679
6: 142967

For 3 and 6 buckets, there is a 2:1 bias.


For a truly uniform distribution, one needs to discard values every once in a
while. Something like the shrink() method below:


#include <iostream>
#include <cstdlib>

template < typename Engine >
class remap_range {
public:

  typedef unsigned long long  value_type;

private:

  typedef Engine engine_type;

  engine_type  the_engine;
  value_type   the_maximum;


  value_type (remap_range::* op) ( void );


  // shrinking the range:
  // ====================
  value_type   bin_size;
  value_type   top_legal_value;

  void set_shrink_pars ( value_type engine_max, value_type range_max ) {
    bin_size = value_type( ( engine_max + 1.0L ) / ( range_max + 1.0L ) );
    value_type num_bins =
      value_type( ( the_engine.max() + 1.0L ) / bin_size );
    top_legal_value = value_type( 1.0L * num_bins * bin_size - 1 );
  }

  value_type shrink ( void ) {
    value_type random_value = the_engine();
    while ( top_legal_value < random_value ) {
      random_value = the_engine();
    }
    return ( random_value / bin_size );
  }

  // expanding the range:
  // ====================
  value_type expand ( void ) {
    return ( remap_to( the_maximum ) );
  }


  // remapping:
  // ==========
  value_type remap_to ( value_type range_max ) {
    if ( range_max <= the_engine.max() ) {
      set_shrink_pars( the_engine.max(), range_max );
      return ( shrink() );
    }
    value_type random_value = 0;
    do {
      random_value =
        value_type
        (
         ( the_engine.max() + 1.0L ) *
         remap_to( value_type( ( range_max + 1.0 )
                               / ( the_engine.max() + 1.0L ) ) )
         +
         the_engine()
         );
    } while ( random_value > range_max );
    return ( random_value );
  }


public:

  remap_range ( Engine eng, value_type max )
    : the_engine  ( eng )
    , the_maximum ( max )
  {
    if ( the_maximum <= the_engine.max() ) {
      set_shrink_pars( the_engine.max(), the_maximum );
      op = &remap_range::shrink;
    } else {
      op = &remap_range::expand;
    }
  }


  value_type max ( void ) const {
    return ( the_maximum );
  }

  value_type operator() ( void ) {
    return ( (this->*op)() );
  }

};



class crand {
public:

  typedef unsigned int value_type;

  crand ( unsigned int seed = 0 ) {
    std::srand( seed );
  }


  value_type max ( void ) const {
    return ( RAND_MAX );
  }

  value_type operator() ( void ) {
    return ( std::rand() );
  }

};


int main ( void ) {
  crand s;
  remap_range< crand > seven ( s, 7 );
  remap_range< remap_range< crand > > thausand ( seven, 20 );
  for ( unsigned int i = 0; i < 100000; ++i ) {
    std::cout << thausand() << '\n';
  }
}


Best

Kai-Uwe Bux


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33815

[Bug libstdc++/33815] tr1::uniform_int isn't uniform

[pcarlini at suse dot de]

------- Comment #12 from pcarlini at suse dot de  2007-10-20 11:27 -------
(In reply to comment #11)

> I just tested this. It does a much better job at obfuscating the bias, which is
> still there (after all you are rescaling linearly). If you rescale a domain of
> size 16 to a range of size 12, some buckets will get 2 and some will get 1.

Of course. What I did is just the naive rescaling, which you can find, for
example, in Numerical Recipes. It assumes, implicitly, that the source range is
much bigger than the target range.

In order to have a fully conforming implementation, we need something which
works for *any* source range and *any* target range. Are you willing to work on
that? Would be mainline-only of course, too large for 4.2.x. In that case you
have also to file a Copyright assignment. Please let me know.


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33815

[Bug libstdc++/33815] tr1::uniform_int isn't uniform

[jkherciueh at gmx dot net]

------- Comment #13 from jkherciueh at gmx dot net  2007-10-21 09:49 -------
(In reply to comment #12)
> In order to have a fully conforming implementation, we need something which
> works for *any* source range and *any* target range. Are you willing to
> work on that? Would be mainline-only of course, too large for 4.2.x. In
> that case you have also to file a Copyright assignment. Please let me know.

Sure thing.

I will need a little bit of time to familiarize myself with the current code.
Also, I will need a little bit of help in finding my way around. E.g., where do
the test-cases for tr1/random go?

Finally, I will need a little bit of guidance with regard to the process of
submitting patches. Is there a web-page for first-timers that describes the
submission process?


Best

Kai-Uwe Bux


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33815

[Bug libstdc++/33815] tr1::uniform_int isn't uniform

[pcarlini at suse dot de]

------- Comment #14 from pcarlini at suse dot de  2007-10-21 11:33 -------
(In reply to comment #13)
> Sure thing.

Excellent.

> I will need a little bit of time to familiarize myself with the current code.
> Also, I will need a little bit of help in finding my way around. E.g., where do
> the test-cases for tr1/random go?

Well in testsuite/tr1/5_numerical_facilities/random, no?

Do you have at hand the specifications? Would be N1836 for TR1 and the last
working paper N2369 for C++0x. Everything from:

 http://www.open-std.org/jtc1/sc22/wg21/docs/papers/

Frankly, I would suggest concentrating more on the C++0x specifications, for
the future: the random facilities will be rather different in C++0x! In TR1
some bits were simply broken in the design, like the idea of variate_generator,
which doesn't exist anymore in C++0x. For TR1 we have to live with that.

> Finally, I will need a little bit of guidance with regard to the process of
> submitting patches. Is there a web-page for first-timers that describes the
> submission process?

Of course you have to read:

  http://gcc.gnu.org/contribute.html

Otherwise, have a look to the patches posted on gcc-patches@ and libstdc++@.
But really, here we care about the algorithm, I can take care of the minor
issues.

Anyway, the MOST important thing at this stage is the Copyright Assignment!
Yours' in any case would be a non-trivial contribution, we need it. Please
start on that as soon as possible, it takes a lot of time. Again, see the above
link. In practice, you should now request the form to assignments@gnu.org.

Finally (I will be disconnected for some days), wanted to mention that, for the
present, I'm thinking about a relatively small change similar to what GSL is
doing (and your idea, I think, but I didn't really read it because I don't want
to plagiarize before your Assignment is in place and then proper credits):

Index: random
===================================================================
--- random      (revision 129509)
+++ random      (working copy)
@@ -1603,17 +1603,7 @@
       template<typename _UniformRandomNumberGenerator>
         result_type
         _M_call(_UniformRandomNumberGenerator& __urng,
-               result_type __min, result_type __max, true_type)
-        {
-         // XXX Must be fixed to also work when __urng.max() - __urng.min()
-         // is smaller than __max - __min.
-         typedef typename __gnu_cxx::__add_unsigned<typename
-           _UniformRandomNumberGenerator::result_type>::__type __utype;
-         return result_type((__max - __min + 1.0L)
-                            * (__utype(__urng()) - __utype(__urng.min()))
-                            / (__utype(__urng.max())
-                               - __utype(__urng.min()) + 1.0L)) + __min;
-       }
+               result_type __min, result_type __max, true_type);

       template<typename _UniformRandomNumberGenerator>
         result_type
Index: random.tcc
===================================================================
--- random.tcc  (revision 129509)
+++ random.tcc  (working copy)
@@ -750,6 +750,31 @@
     }


+  template<typename _IntType>
+    template<typename _UniformRandomNumberGenerator>
+      typename uniform_int<_IntType>::result_type
+      uniform_int<_IntType>::
+      _M_call(_UniformRandomNumberGenerator& __urng,
+             result_type __min, result_type __max, true_type)
+      {
+       // XXX Must be fixed to also work when __urng.max() - __urng.min()
+       // is smaller than __max - __min.
+       typedef typename __gnu_cxx::__add_unsigned<typename
+         _UniformRandomNumberGenerator::result_type>::__type __utype;
+
+       result_type __ret;
+
+       __utype __umin = __urng.min();
+       __utype __umax = __urng.max();
+       __utype __denom = (__umax - __umin) / (__max - __min + 1);
+
+       do
+         __ret = (__utype(__urng()) -  __umin) / __denom;
+       while (__ret > __max - __min);
+
+       return __ret + __min;
+      }
+
   template<typename _IntType, typename _CharT, typename _Traits>
     std::basic_ostream<_CharT, _Traits>&
     operator<<(std::basic_ostream<_CharT, _Traits>& __os,


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33815

[Bug libstdc++/33815] tr1::uniform_int isn't uniform

[paolo at gcc dot gnu dot org]

------- Comment #15 from paolo at gcc dot gnu dot org  2007-10-30 13:05 -------
Subject: Bug 33815

Author: paolo
Date: Tue Oct 30 13:05:26 2007
New Revision: 129769

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=129769
Log:
2007-10-19  Paolo Carlini  <pcarlini@suse.de>

        PR libstdc++/33815
        * include/tr1_impl/random
        (uniform_int<>::_M_call(_UniformRandomNumberGenerator&, result_type,
        result_type, true_type)): Avoid the modulo (which uses the low-order
        bits).

Modified:
    trunk/libstdc++-v3/ChangeLog
    trunk/libstdc++-v3/include/tr1_impl/random
    trunk/libstdc++-v3/include/tr1_impl/random.tcc


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33815