QRandomGenerator: improve floating-point random generation

The previous version was good, just not optimal. Because the input was
an unsigned 64-bit number, compilers needed to generate extra code to
deal with HW instructions that only convert 64-bit signed input. And
that was useless because a double uniformly distributed from 0 to 1 can
only have 53 bits of randomness.

The previous implementation did exactly what the Microsoft libstdc++ and
libc++ implementations do. In my opinion, those implementations have an
imperfect distribution, which is corrected in this commit. In those, all
random input bigger than 0x20000000000000 has a different frequency
compared to input below that mark. For example, both 0x20000000000000
and 0x20000000000001 produce the same result (4.8828125e-4).

What's more, for the libc++ and MSVC implementations, input between
0xfffffffffffff001 and 0xffffffffffffffff results in 1.0 (probability 1
in 2⁵³), even though the Standard is very clear that the result should
be strictly less than 1. GCC 7's libstdc++ doesn't have this issue,
whereas the versions before would enter an infinite loop.

Change-Id: Ib17dde1a1dbb49a7bba8fffd14eced3c375dd2ec
Reviewed-by: Lars Knoll <lars.knoll@qt.io>
Reviewed-by: Edward Welbourne <edward.welbourne@qt.io>

src/corelib/global/qrandom.h

@@ -62,8 +62,16 @@ public:
     static Q_CORE_EXPORT quint64 generate64();
     static double generateDouble()
     {
-        // use get64() to get enough bits
-        return double(generate64()) / ((std::numeric_limits<quint64>::max)() + double(1.0));
+        // IEEE 754 double precision has:
+        //   1 bit      sign
+        //  10 bits     exponent
+        //  53 bits     mantissa
+        // In order for our result to be normalized in the range [0, 1), we
+        // need exactly 53 bits of random data. Use generate64() to get enough.
+        quint64 x = generate64();
+        quint64 limit = Q_UINT64_C(1) << std::numeric_limits<double>::digits;
+        x >>= std::numeric_limits<quint64>::digits - std::numeric_limits<double>::digits;
+        return double(x) / double(limit);
     }
 
     static qreal bounded(qreal sup)

tests/auto/corelib/global/qrandomgenerator/tst_qrandomgenerator.cpp

@@ -97,6 +97,8 @@ private slots:
 
     void generateReal_data() { generate32_data(); }
     void generateReal();
+    void qualityReal_data() { generate32_data(); }
+    void qualityReal();
 
     void seedStdRandomEngines();
     void stdUniformIntDistribution_data();
@@ -445,7 +447,7 @@ void tst_QRandomGenerator::generateReal()
     for (int i = 0; i < 4; ++i) {
         QVERIFY_3TIMES([] {
             qreal value = QRandomGenerator::generateDouble();
-            return value > 0 && value < 1 && value != RandomValueFP;
+            return value >= 0 && value < 1 && value != RandomValueFP;
         }());
     }
 
@@ -454,6 +456,52 @@ void tst_QRandomGenerator::generateReal()
         QVERIFY_3TIMES(QRandomGenerator::generateDouble() != QRandomGenerator::generateDouble());
 }
 
+void tst_QRandomGenerator::qualityReal()
+{
+    QFETCH(uint, control);
+    setRNGControl(control);
+
+    enum {
+        SampleSize = 160,
+
+        // Expected value: sample size times proportion of the range:
+        PerfectOctile = SampleSize / 8,
+        PerfectHalf = SampleSize / 2,
+
+        // Variance is (1 - proportion of range) * expected; sqrt() for standard deviations.
+        // Should usually be within twice that and almost never outside four times:
+        RangeHalf = 25,         // floor(4 * sqrt((1 - 0.5) * PerfectHalf))
+        RangeOctile = 16        // floor(4 * sqrt((1 - 0.125) * PerfectOctile))
+    };
+
+    double data[SampleSize];
+    std::generate(std::begin(data), std::end(data), &QRandomGenerator::generateDouble);
+
+    int aboveHalf = 0;
+    int belowOneEighth = 0;
+    int aboveSevenEighths = 0;
+    for (double x : data) {
+        aboveHalf += x >= 0.5;
+        belowOneEighth += x < 0.125;
+        aboveSevenEighths += x >= 0.875;
+
+        // these are strict requirements
+        QVERIFY(x >= 0);
+        QVERIFY(x < 1);
+    }
+
+    qInfo("Halfway distribution: %.1f - %.1f", 100. * aboveHalf / SampleSize, 100 - 100. * aboveHalf / SampleSize);
+    qInfo("%.1f below 1/8 (expected 12.5%% ideally)", 100. * belowOneEighth / SampleSize);
+    qInfo("%.1f above 7/8 (expected 12.5%% ideally)", 100. * aboveSevenEighths / SampleSize);
+
+    QVERIFY(aboveHalf < PerfectHalf + RangeHalf);
+    QVERIFY(aboveHalf > PerfectHalf - RangeHalf);
+    QVERIFY(aboveSevenEighths < PerfectOctile + RangeOctile);
+    QVERIFY(aboveSevenEighths > PerfectOctile - RangeOctile);
+    QVERIFY(belowOneEighth < PerfectOctile + RangeOctile);
+    QVERIFY(belowOneEighth > PerfectOctile - RangeOctile);
+}
+
 template <typename Engine> void seedStdRandomEngine()
 {
     QRandomGenerator rd;