Removed skipping of 0 (zero) in TRandom3.Rndm()

math/mathcore/src/TRandom3.cxx

@@ -7,27 +7,26 @@
 
 Random number generator class based on
   M. Matsumoto and T. Nishimura,
-  Mersenne Twister: A 623-diminsionally equidistributed
+  Mersenne Twister: A 623-dimensionally equidistributed
   uniform pseudorandom number generator
   ACM Transactions on Modeling and Computer Simulation,
   Vol. 8, No. 1, January 1998, pp 3--30.
 
 For more information see the Mersenne Twister homepage
   [http://www.math.keio.ac.jp/~matumoto/emt.html]
 
-Advantage:
+Advantages:
+- large period 2**19937 - 1
+- relatively fast (slightly slower than TRandom2 but much faster than TRandom1)
 
--  large period 2**19937 -1
--  relatively fast (slightly slower than TRandom2 but much faster than TRandom1)
-
-Drawback:
--  a relative large internal state of 624 integers 
-- generate only 32 random bits 
+Drawbacks:
+- a relative large internal state of 624 integers
+- generate only 32 random bits
 - not passing all the random generator tests. It fails some tests in TestU01
- (see [http://simul.iro.umontreal.ca/testu01/tu01.html])
+    (see [http://simul.iro.umontreal.ca/testu01/tu01.html])
 
-An altenativly excellent generator passing all tests of TestU01, having 61 random bits and 
-fast as Mersenne and Twister is MIXMAX (TRandomMixMax). 
+An alternatively excellent generator passing all tests of TestU01, having 61 random bits and
+being as fast as Mersenne and Twister is MIXMAX (TRandomMixMax).
 
 @ingroup Random
 
@@ -58,7 +57,6 @@ fast as Mersenne and Twister is MIXMAX (TRandomMixMax).
 #include "TRandom3.h"
 #include "TBuffer.h"
 #include "TRandom2.h"
-#include "TUUID.h"
 
 TRandom *gRandom = new TRandom3();
 #ifdef R__COMPLETE_MEM_TERMINATION
@@ -73,8 +71,10 @@ namespace {
 ClassImp(TRandom3);
 
 ////////////////////////////////////////////////////////////////////////////////
-/// Default constructor
-/// If seed is 0, the seed is automatically computed via a TUUID object.
+/// \brief Default constructor.
+///
+/// If seed is 0, the seed array is automatically computed via a TRandom2
+/// object, which internally uses TUUID.
 /// In this case the seed is guaranteed to be unique in space and time.
 
 TRandom3::TRandom3(UInt_t seed)
@@ -85,21 +85,19 @@ TRandom3::TRandom3(UInt_t seed)
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-/// Default destructor
+/// \brief Default destructor.
 
 TRandom3::~TRandom3()
 {
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-///  Machine independent random number generator.
-///  Produces uniformly-distributed floating points in (0,1)
+///  \brief Machine independent random number generator.
+///
+///  Produces uniformly-distributed floating points in ]0, 1].
 ///  Method: Mersenne Twister
-
 Double_t TRandom3::Rndm()
 {
-   UInt_t y;
-
    const Int_t  kM = 397;
    const Int_t  kN = 624;
    const UInt_t kTemperingMaskB =  0x9d2c5680;
@@ -108,6 +106,7 @@ Double_t TRandom3::Rndm()
    const UInt_t kLowerMask =       0x7fffffff;
    const UInt_t kMatrixA =         0x9908b0df;
 
+   UInt_t y;
    if (fCount624 >= kN) {
       Int_t i;
 
@@ -132,28 +131,25 @@ Double_t TRandom3::Rndm()
    y ^= ((y << 15) & kTemperingMaskC );
    y ^=  (y >> 18);
 
-   // 2.3283064365386963e-10 == 1./(max<UINt_t>+1)  -> then returned value cannot be = 1.0
-   if (y) return ( (Double_t) y * 2.3283064365386963e-10); // * Power(2,-32)
-   return Rndm();
+   // 2.3283064365386963e-10 == 1. / (max<UInt_t> + 1)
+   return (static_cast<Double_t>(y) + 1) * 2.3283064365386963e-10;
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-/// Return an array of n random numbers uniformly distributed in ]0,1]
+///  \brief Return an array of n random numbers uniformly distributed in ]0, 1].
 
 void TRandom3::RndmArray(Int_t n, Float_t *array)
 {
-  for(Int_t i=0; i<n; i++) array[i]=(Float_t)Rndm();
+   for (Int_t i = 0; i < n; i++) {
+      array[i] = static_cast<Float_t>(Rndm());
+   }
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-/// Return an array of n random numbers uniformly distributed in ]0,1]
+///  \brief Return an array of n random numbers uniformly distributed in ]0, 1].
 
 void TRandom3::RndmArray(Int_t n, Double_t *array)
 {
-   Int_t k = 0;
-
-   UInt_t y;
-
    const Int_t  kM = 397;
    const Int_t  kN = 624;
    const UInt_t kTemperingMaskB =  0x9d2c5680;
@@ -162,6 +158,8 @@ void TRandom3::RndmArray(Int_t n, Double_t *array)
    const UInt_t kLowerMask =       0x7fffffff;
    const UInt_t kMatrixA =         0x9908b0df;
 
+   Int_t k = 0;
+   UInt_t y;
    while (k < n) {
       if (fCount624 >= kN) {
          Int_t i;
@@ -187,21 +185,21 @@ void TRandom3::RndmArray(Int_t n, Double_t *array)
       y ^= ((y << 15) & kTemperingMaskC );
       y ^=  (y >> 18);
 
-      if (y) {
-         array[k] = Double_t( y * 2.3283064365386963e-10); // * Power(2,-32)
-         k++;
-      }
+      array[k] = (static_cast<Double_t>(y) + 1) * 2.3283064365386963e-10; // * Power(2,-32)
+      k++;
    }
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-///  Set the random generator sequence
-/// if seed is 0 (default value) a TUUID is generated and used to fill
-/// the first 8 integers of the seed array.
+/// \brief Set the random generator sequence.
+///
+/// If seed is 0 (default value) a TRandom2 (internally uses TUUID) is used to
+/// generate all 624 unsigned integers of the seed array.
 /// In this case the seed is guaranteed to be unique in space and time.
-/// Use upgraded seeding procedure to fix a known problem when seeding with values
-/// with many zero in the bit pattern (like 2**28).
-/// see http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
+///
+/// Upgraded seeding procedure is used to fix a known problem when seeding with
+/// values with many zero in the bit pattern (like 2**28), see
+/// http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
 
 void TRandom3::SetSeed(ULong_t seed)
 {
@@ -217,12 +215,12 @@ void TRandom3::SetSeed(ULong_t seed)
 
    } else {
 
-      // use TRandom2 (which is based on TUUId to generate the seed
-      // TRandom2 works fairly well  and has been tested against example
+      // use TRandom2 (which is based on TUUID to generate the seed.
+      // TRandom2 works fairly well and has been tested against example
       // layout in https://savannah.cern.ch/bugs/?99516
       TRandom2 r(0);
       for (Int_t i = 0; i< 624; i++) {
-         fMt[i]   = static_cast<UInt_t> (4294967296.*r.Rndm());
+         fMt[i] = static_cast<UInt_t>(4294967296. * r.Rndm());
       }
       // warm up the generator calling it 10 times
       for (Int_t i = 0; i < 10; ++i) Rndm();
@@ -232,7 +230,7 @@ void TRandom3::SetSeed(ULong_t seed)
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-/// Stream an object of class TRandom3.
+/// \brief Streamer for an object of class TRandom3.
 
 void TRandom3::Streamer(TBuffer &R__b)
 {

math/mathcore/test/CMakeLists.txt

@@ -96,6 +96,8 @@ ROOT_ADD_GTEST(testRootFinder testRootFinder.cxx  LIBRARIES ${Libraries})
 ROOT_ADD_GTEST(testKahan testKahan.cxx LIBRARIES Core MathCore)
 ROOT_ADD_GTEST(testDelaunay2D testDelaunay2D.cxx LIBRARIES Core MathCore)
 
+ROOT_ADD_GTEST(testTRandom3vsStdMt19937 testTRandom3vsStdMt19937.cxx LIBRARIES Core MathCore)
+
 if(clad)
   ROOT_ADD_GTEST(CladDerivatorTests CladDerivatorTests.cxx LIBRARIES Core MathCore)
 endif()

math/mathcore/test/testTRandom3vsStdMt19937.cxx

@@ -0,0 +1,31 @@
+// std
+#include <random>
+
+// ROOT
+#include "TRandom3.h"
+
+#include "gtest/gtest.h"
+
+TEST(TRandom3Test, TRandom3vsStdMt19937)
+{
+   unsigned int seed = 28u;
+
+   std::mt19937 gen_std(seed);
+   TRandom3 gen_root(seed);
+
+   // Test TRandom3 generator against other implementation of Mersenne Twister
+   // 32-bit algoritm.
+   //
+   // Test inspired by code of Giovanni Cerretani in this issue:
+   // https://its.cern.ch/jira/browse/ROOT-9733
+   // https://github.com/root-project/root/issues/14581
+   for (size_t generation = 0; generation < 2; ++generation) {
+      unsigned int iteration = 0;
+      while (++iteration) {
+         double rnd_std = (gen_std() + 1) * 2.3283064365386963e-10;
+         double rnd_root = gen_root.Rndm();
+
+         ASSERT_EQ(rnd_std, rnd_root);
+      }
+   }
+}