Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 8 commits into
base: master
Choose a base branch
from
Open

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

Changes from 1 commit
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
daf21fc
Removed skipping of 0 (zero) in TRandom3.Rndm()
kjvbrt Feb 13, 2024
d869c06
Adding test
kjvbrt Feb 14, 2024
5db6b75
Small fixes
kjvbrt Feb 14, 2024
6ed8fba
Formatting
kjvbrt Feb 14, 2024
d98d6af
Rndm() now generates ]0, 1]
kjvbrt Feb 14, 2024
e69fa9d
Formatting
kjvbrt Feb 14, 2024
791c755
Runing only over two overflows
kjvbrt Feb 14, 2024
f501b84
Removed 0 skipping also in case of RndmArray
kjvbrt Feb 14, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
62 changes: 32 additions & 30 deletions math/mathcore/src/TRandom3.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand All @@ -74,7 +72,7 @@ ClassImp(TRandom3);

////////////////////////////////////////////////////////////////////////////////
/// Default constructor
/// If seed is 0, the seed is automatically computed via a TUUID object.
/// If seed is 0, the seed array is automatically computed via a TUUID object.
/// In this case the seed is guaranteed to be unique in space and time.

TRandom3::TRandom3(UInt_t seed)
Expand All @@ -85,21 +83,21 @@ TRandom3::TRandom3(UInt_t seed)
}

////////////////////////////////////////////////////////////////////////////////
/// Default destructor
/// Default destructor.

TRandom3::~TRandom3()
{
}

////////////////////////////////////////////////////////////////////////////////
/// Machine independent random number generator.
/// Produces uniformly-distributed floating points in (0,1)
/// \brief Machine independent random number generator.
///
/// \warning Can produce 0 (zero).
///
/// 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;
Expand All @@ -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;

Expand All @@ -132,21 +131,22 @@ 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) -> then returned value cannot be = 1.0
return static_cast<Double_t>(y * 2.3283064365386963e-10);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We can add the +1 here after converting to double to avoid overflows in the uint:

return (static_cast<Double_t>(y) + 1)* 2.3283064365386963e-10)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The overhead of adding a +1 should be minimal. We could test running tutorial/math/testRandom.C

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@kjvbrt could you take care?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I'm testing the implementation with +1 ATM.

Copy link
Contributor Author

@kjvbrt kjvbrt Feb 14, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is the result:

                   TRandom3    TRandom3 with +1       
Rndm..............   20.000      20.000                                               
RndmArray.........   10.000      15.000                                               
Gaus..............   60.000      60.000                                               
Rannor............   40.000      40.000                                               
Landau............   40.000      40.000                                               
Exponential.......   35.000      35.000
Binomial(5,0.5)...  125.000     130.000
Binomial(15,0.5)..  370.000     380.000
Poisson(3)........   90.000      95.000
Poisson(10).......  205.000     210.000
Poisson(70).......  295.000     300.000
Poisson(100)......  295.000     295.000
GausTF1...........  180.000     180.000
LandauTF1.........  175.000     175.000

}

////////////////////////////////////////////////////////////////////////////////
/// 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[
/// \warning Can produce 0 (zero).

void TRandom3::RndmArray(Int_t n, Float_t *array)
{
for(Int_t i=0; i<n; i++) array[i]=(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[
/// \warning Can produce 0 (zero).

void TRandom3::RndmArray(Int_t n, Double_t *array)
{
Expand Down Expand Up @@ -195,13 +195,15 @@ void TRandom3::RndmArray(Int_t n, Double_t *array)
}

////////////////////////////////////////////////////////////////////////////////
/// 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)
{
Expand Down