Fix base::RandGenerator bug (it had non-uniform random distribution).

base/rand_util_unittest.cc

 // Copyright (c) 2011 The Chromium Authors. All rights reserved.
 // Use of this source code is governed by a BSD-style license that can be
 // found in the LICENSE file.
 
 #include "base/rand_util.h"
 
 #include <limits>
 
 #include "testing/gtest/include/gtest/gtest.h"
 
@@ skipping 43 matching lines @@
   EXPECT_NE(0, accumulator);
 }
 
 // Make sure that it is still appropriate to use RandGenerator in conjunction
 // with std::random_shuffle().
 TEST(RandUtilTest, RandGeneratorForRandomShuffle) {
   EXPECT_EQ(base::RandGenerator(1), 0U);
   EXPECT_LE(std::numeric_limits<ptrdiff_t>::max(),
             std::numeric_limits<int64>::max());
 }
+
+TEST(RandUtilTest, RandGeneratorIsUniform) {
+  // Verify that RandGenerator has a uniform distribution. This is a
+  // regression test that consistently failed when RandGenerator was
+  // implemented this way:
+  //
+  //   return base::RandUint64() % max;
+  //
+  // The worst case for such an implementation is e.g. a top of range
+  // that is 2/3rds of the way to MAX_UINT64, in which case the bottom

[jar (doing other things), 2011/05/28 17:26:40]

You got me to wondering about this ratio, and it wasn't clear that 2/3 was
optimally problematic.  I applied calculus (which I've only done about 4 times
in my life, after gettin' out'a school), and concluded that the maximum delta
between the broken-code-mean and the correct (range over 2 mean) is actually
achieved with a ratio R of:

3/4

(If I'm wrong, please correct me!  I'm notorious for clever ideas... and typos
in my arithmetic!).

As always, after grinding through the calculus, there is an easy way to see this
result (that the delta in the means is maximized with R=3/4)

We only consider R in the range of 1/2 to 1, so that we never have to consider
cases where the rand() values mapped more than twice into any return value.

The difference function is then a quadratic in R, where R is the ratio.  You can
pretty much (in retrospect) guess that since the calculation of the means is
done by summing all numbers in a range, and dividing by the count of numbers. 
We just have two different ranges that we need to sum (one after the modular
re-mapping).  Everything in that equation is based on R and (1-R), and the sum
of numbers from 1 to K is quadratic.

So you know it is *some* quadratic (i.e., a parabola)

If the ratio was 1/2, then there would be no difference in the mean (no bias
would be introduced).  Similarly, if the ratio were 1, there would be no bias. 
Hence we know that this parabola has a zero value at R=1/2, and R = 1.

Since parabolas are nice and symetric, the the zero-derivative (max) point is
half-way between those two equal values.

The precise equation (if I didn't screw up) for the difference in the means is:

MaxInt  * (R^2  -3R/2 + 1/2)

Again, feel free to correct my math.

+  // half of the range would be twice as likely to occur as the top
+  // half, assuming a naive modulus implementation of RandGenerator.
+  const uint64 kTopOfRange = (std::numeric_limits<uint64>::max() / 3L) * 2L;
+  const uint64 kExpectedAverage = kTopOfRange / 2L;
+  const uint64 kAllowedVariance = kExpectedAverage / 50L;  // +/- 2%
+  const int kMinAttempts = 1000;
+  const int kMaxAttempts = 1000000;
+
+  double cumulative_average = 0.0;
+  int count = 0;
+  while (count < kMaxAttempts) {
+    uint64 value = base::RandGenerator(kTopOfRange);
+    cumulative_average = (count * cumulative_average + value) / (count + 1);
+
+    // Don't quit too quickly for things to start converging, or we may have
+    // a false positive.
+    if (count > kMinAttempts &&
+        kExpectedAverage - kAllowedVariance < cumulative_average &&
+        cumulative_average < kExpectedAverage + kAllowedVariance) {
+      break;
+    }
+
+    ++count;
+  }
+
+  ASSERT_LT(count, kMaxAttempts) << "Expected average was " <<
+      kExpectedAverage << ", average ended at " << cumulative_average;
+}