Skip to content

GitLab

Commit 8283936b authored by Mark Fishcler's avatar Mark Fishcler
Browse files

Code to handle the case where RAND_MAX is not one of the two "expected" 

values of 2**15-1 or 2**31-1.
parent e6cb23b6
Hide whitespace changes
Inline Side-by-side
Random/Random/RandEngine.h
View file @ 8283936b
// $Id: RandEngine.h,v 1.3 2003/10/23 21:29:51 garren Exp $
// $Id: RandEngine.h,v 1.3.2.1 2004/04/29 18:43:45 fischler Exp $
// -*- C++ -*-
//
// -----------------------------------------------------------------------
......@@ -25,6 +25,8 @@
// engine counter: 15th Feb 1998
// Ken Smith - Added conversion operators: 6th Aug 1998
// replace mx by mantissa_bit_32
// M Fischler - Inserted warnings about the fact that the quality of rand()
// is quite poor.
// =======================================================================
#ifndef RandEngine_h
......@@ -54,6 +56,10 @@ public:
// It returns a pseudo random number between 0 and 1,
// according to the standard stdlib random function rand()
// but excluding the end points.
//
// WARNING: rand() is quite a weak generator on most systems, <
// will not pass several randomness tests, and does not give a <
// reproducible sequence of numbers. <
void flatArray (const int size, double* vect);
// Fills the array "vect" of specified size with flat random values.
......@@ -67,10 +73,14 @@ public:
void saveStatus( const char filename[] = "Rand.conf" ) const;
// Saves on file Rand.conf the current engine status.
// WARNING: This is non-functional, as rand() on various systems will <
// not give reproducible streams. <
void restoreStatus( const char filename[] = "Rand.conf" );
// Reads from file Rand.conf the last saved engine status
// and restores it.
// WARNING: This is non-functional, as rand() on various systems will <
// not give reproducible streams. <
void showStatus() const;
// Dumps the engine status on the screen.
......
Random/src/RandEngine.cc
View file @ 8283936b
// $Id: RandEngine.cc,v 1.4.2.1 2004/04/28 16:04:57 garren Exp $
// $Id: RandEngine.cc,v 1.4.2.2 2004/04/29 18:43:45 fischler Exp $
// -*- C++ -*-
//
// -----------------------------------------------------------------------
......@@ -22,6 +22,15 @@
// obtained by concatenation: 15th Feb 1998
// Ken Smith - Added conversion operators: 6th Aug 1998
// J. Marraffino - Remove dependence on hepString class 13 May 1999
// M. Fischler - Rapaired bug that in flat() that relied on rand() to
// deliver 15-bit results. This bug was causing flat()
// on Linux systems to yield randoms with mean of 5/8(!)
// - Modified algorithm such that on systems with 31-bit rand()
// it will call rand() only once instead of twice. Feb 2004
// M. Fischler - Modified the general-case template for RandEngineBuilder
// such that when RAND_MAX is an unexpected value the routine
// will still deliver a sensible flat() random.
//
// =======================================================================
#include "CLHEP/Random/defs.h"
......@@ -33,6 +42,21 @@
namespace CLHEP {
#ifdef NOTDEF
// The way to test for proper behavior of the RandEngineBuilder
// for arbitrary RAND_MAX, on a system where RAND_MAX is some
// fixed specialized value and rand() behaves accordingly, is
// to set up a fake RAND_MAX and a fake version of rand()
// by enabling this block.
#undef RAND_MAX
#define RAND_MAX 9382956
#include "CLHEP/Random/MTwistEngine.h"
#include "CLHEP/Random/RandFlat.h"
MTwistEngine * fakeFlat = new MTwistEngine;
RandFlat rflat (fakeFlat, 0, RAND_MAX+1);
int rand() { return (int)rflat(); }
#endif
static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
// number of instances with automatic seed selection
......@@ -166,19 +190,140 @@ void RandEngine::showStatus() const
std::cout << "----------------------------------------" << std::endl;
}
double RandEngine::flat()
{
// rand() is such a horrible generator, it only generates 15 bits.
// The quick fix here to get it to at LEAST 32 bits is to grab two values,
// and concatenate them, but that still leaves two bits empty and these
// we grab from the number of times flat() has been called, i.e. seq.
return double( (unsigned int)( (rand() << 17) | // bits 31-17
(seq++ & 0x3 << 15) | // bits 16,15
rand() // bits 14-0
) * mantissa_bit_32
);
}
// ====================================================
// Implementation of flat() (along with needed helpers)
// ====================================================
// Here we set up such that **at compile time**, the compiler decides based on
// RAND_MAX how to form a random double with 32 leading random bits, using
// one or two calls to rand(). Some of the lowest order bits of 32 are allowed
// to be as weak as mere XORs of some higher bits, but not to be always fixed.
//
// The method decision is made at compile time, rather than using a run-time
// if on the value of RAND_MAX. Although it is possible to cope with arbitrary
// values of RAND_MAX of the form 2**N-1, with the same efficiency, the
// template techniques needed would look mysterious and perhaps over-stress
// some older compilers. We therefore only treat RAND_MAX = 2**15-1 (as on
// most older systems) and 2**31-1 (as on the later Linux systems) as special
// and super-efficient cases. We detect any different value, and use an
// algorithm which is correct even if RAND_MAX is not one less than a power
// of 2.
template <int> struct RandEngineBuilder { // RAND_MAX any arbitrary value
static unsigned int thirtyTwoRandomBits() {
static bool prepared = false;
static unsigned int iT;
static unsigned int iK;
static unsigned int iS;
static int iN;
static double fS;
static double fT;
if ( (RAND_MAX >> 31) > 0 )
{
// Here, we know that integer arithmetic is 64 bits.
if ( !prepared ) {
iS = RAND_MAX + 1;
iK = 1;
// int StoK = S;
int StoK = iS;
if ( (RAND_MAX >> 32) == 0) {
iK = 2;
// StoK = S*S;
StoK = iS*iS;
}
int m;
for ( m = 0; m < 64; ++m ) {
StoK >>= 1;
if (StoK == 0) break;
}
iT = 1 << m;
prepared = true;
}
int v = 0;
do {
for ( int i = 0; i < iK; ++i ) {
v = iS*v+rand();
}
} while (v < iT);
return v & 0xFFFFFFFF;
}
else if ( (RAND_MAX >> 26) == 0 )
{
// Here, we know it is safe to work in terms of doubles without loss
// of precision, but we have no idea how many randoms we will need to
// generate 32 bits.
if ( !prepared ) {
fS = RAND_MAX + 1;
double twoTo32 = ldexp(1.0,32);
double StoK = fS;
for ( iK = 1; StoK < twoTo32; StoK *= fS, iK++ ) { }
int m;
fT = 1.0;
for ( m = 0; m < 64; ++m ) {
StoK *= .5;
if (StoK < 1.0) break;
fT *= 2.0;
}
prepared = true;
}
double v = 0;
do {
for ( int i = 0; i < iK; ++i ) {
v = fS*v+rand();
}
} while (v < fT);
return ((unsigned int)v) & 0xFFFFFFFF;
}
else
{
// Here, we know that 16 random bits are available from each of
// two random numbers.
if ( !prepared ) {
iS = RAND_MAX + 1;
int SshiftN = iS;
for (iN = 0; SshiftN > 1; SshiftN >>= 1, iN++) { }
iN -= 17;
prepared = true;
}
unsigned int x1, x2;
do { x1 = rand(); } while (x1 < (1<<16) );
do { x2 = rand(); } while (x2 < (1<<16) );
return x1 | (x2 << 16);
}
}
};
template <> struct RandEngineBuilder<2147483647> { // RAND_MAX = 2**31 - 1
inline static unsigned int thirtyTwoRandomBits() {
unsigned int x = rand() << 1; // bits 31-1
x ^= ( (x>>23) ^ (x>>7) ) ^1; // bit 0 (weakly pseudo-random)
return x & 0xFFFFFFFF; // mask in case int is 64 bits
}
};
template <> struct RandEngineBuilder<32767> { // RAND_MAX = 2**15 - 1
inline static unsigned int thirtyTwoRandomBits() {
unsigned int x = rand() << 17; // bits 31-17
x ^= rand() << 2; // bits 16-2
x ^= ( (x>>23) ^ (x>>7) ) ^3; // bits 1-0 (weakly pseudo-random)
return x & 0xFFFFFFFF; // mask in case int is 64 bits
}
};
double RandEngine::flat()
{
double r;
do { r = RandEngineBuilder<RAND_MAX>::thirtyTwoRandomBits();
} while ( r == 0 );
return r/4294967296.0;
}
void RandEngine::flatArray(const int size, double* vect)
{
......@@ -189,7 +334,7 @@ void RandEngine::flatArray(const int size, double* vect)
}
RandEngine::operator unsigned int() {
return (unsigned int)( (rand() << 17) | ((seq++&0x3) << 15) | rand() );
return RandEngineBuilder<RAND_MAX>::thirtyTwoRandomBits();
}
std::ostream & operator << ( std::ostream& os, const RandEngine& e )
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment