Commit b85d0b5c authored by Mark Olesen's avatar Mark Olesen
Browse files

ENH: provide Rand48 as generator in the expected C++11 form

- this removes an OS-specific dependency (eg, drand48_r is not POSIX)
  and allows easier use of other random number generators.

  The Rand48 generator has identical behaviour and period as the
  lrand48() library routine, but holds its own seed and state
  (which makes it re-entrant) and can be combined with other
  random distributions.

  However, when using the modified form to obtain scalar values
  they will not be identical to what drand48() yields.

  This is because drand48() uses the raw 48-bit values to directly
  set the mantissa of an IEEE double where as the newer distribution
  normalizes based on the 32-bit value.

STYLE: simplify code in Random::shuffle and use Swap
parent e4f5471e
Test-Random.C
EXE = $(FOAM_USER_APPBIN)/Test-Random
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2018 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
Test-Random
Description
Simple test for sequence of random numbers
\*---------------------------------------------------------------------------*/
#include "Rand48.H"
#include "Random.H"
#include <cstdlib>
#include <iostream>
#include <iomanip>
#define TEST_RAW_IEEE
// Construct a positive double with the 48 random bits distributed over
// its fractional part so the resulting FP number is [0.0,1.0).
//
// As per glibc erand48() implementation
#ifdef TEST_RAW_IEEE
#include <ieee754.h>
double randomFraction(const uint64_t bits)
{
// 48-bit value
unsigned short int xsubi[3];
xsubi[0] = (bits & 0xffff);
xsubi[1] = ((bits >> 16) & 0xffff);
xsubi[2] = ((bits >> 32) & 0xffff);
union ieee754_double temp;
temp.ieee.negative = 0;
temp.ieee.exponent = IEEE754_DOUBLE_BIAS;
temp.ieee.mantissa0 = (xsubi[2] << 4) | (xsubi[1] >> 12);
temp.ieee.mantissa1 = ((xsubi[1] & 0xfff) << 20) | (xsubi[0] << 4);
// The lower 4 bits of mantissa1 are always 0.
return temp.d - 1.0;
}
#endif
using namespace Foam;
// Output with cout instead of Info to keep values unsigned on output
using std::cout;
using std::setw;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
// Allow multiple passes to ensure reset is working properly
const label maxIter = 1;
std::default_random_engine deflt(123456);
std::mt19937 mtwist(123456);
std::uniform_real_distribution<scalar> uniform01;
{
Rand48 rnd(123456);
for (label iter=0; iter < maxIter; ++iter)
{
rnd.seed(123456);
::srand48(123456);
mtwist.seed(123456);
Info<< nl << "32-bit random with seed = 123456" << nl;
cout<< setw(12) << "Rand48()"
<< setw(12) << "lrand48()"
<< setw(12) << "mtwister"
<< setw(12) << "default" << nl;
for (int i=0; i<25; i++)
{
cout<< setw(12) << rnd()
<< setw(12) << long(::lrand48())
<< setw(12) << long(mtwist())
<< setw(12) << long(deflt()) << nl;
}
}
}
{
Random rnd(123456);
Rand48 manual(123456);
mtwist.seed(123456);
deflt.seed(123456);
// Two passes to ensure that reset is working properly
for (label iter=0; iter < maxIter; ++iter)
{
::srand48(123456);
rnd.reset(123456);
manual.seed(123456);
mtwist.seed(123456);
deflt.seed(123456);
cout<< nl << "Random (Rand48) with seed = " << rnd.seed()
<< " interval [0,1000]" << nl;
cout<< setw(12) << "Rand48()"
<< setw(12) << "drand48()";
#ifdef TEST_RAW_IEEE
cout<< setw(12) << "manual";
#endif
cout<< setw(12) << "mtwister";
cout<< nl;
for (int i=0; i<25; i++)
{
cout<< setw(12) << (rnd.sample01<scalar>()*1000)
<< setw(12) << (drand48()*1000);
#ifdef TEST_RAW_IEEE
cout<< setw(12) << (randomFraction(manual.raw())*1000);
#endif
cout<< setw(12) << (uniform01(mtwist)*1000);
cout<< setw(12) << (uniform01(deflt)*1000);
cout<< nl;
}
}
}
{
Rand48 rnd1(123456);
Rand48 rnd2(123456);
::srand48(123456);
rnd2.discard(10);
cout<< nl << "Rand48 - test with offset of 10" << nl;
cout<< setw(12) << "Rand48()"
<< setw(12) << "offset-10"
<< setw(12) << "lrand48()"
<< nl;
for (int i=0; i<25; i++)
{
cout<< setw(12) << (rnd1())
<< setw(12) << (rnd2())
<< setw(12) << long(::lrand48())
<< nl;
}
}
cout<< nl << "Done." << endl;
return 0;
}
// ************************************************************************* //
......@@ -66,14 +66,6 @@ Description
#include <link.h>
#endif
#ifdef USE_RANDOM
#include <climits>
#if INT_MAX != 2147483647
#error "INT_MAX != 2147483647"
#error "The random number generator may not work!"
#endif
#endif
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
......@@ -1666,84 +1658,6 @@ Foam::fileNameList Foam::dlLoaded()
return libs;
}
Foam::label Foam::osRandomBufferSize()
{
#ifdef USE_RANDOM
return sizeof(random_data);
#else
return sizeof(drand48_data);
#endif
}
void Foam::osRandomSeed(const label seed)
{
#ifdef USE_RANDOM
srandom((unsigned int)seed);
#else
srand48(seed);
#endif
}
Foam::label Foam::osRandomInteger()
{
#ifdef USE_RANDOM
return random();
#else
return lrand48();
#endif
}
Foam::scalar Foam::osRandomDouble()
{
#ifdef USE_RANDOM
return (scalar)random()/INT_MAX;
#else
return drand48();
#endif
}
void Foam::osRandomSeed(const label seed, List<char>& buffer)
{
#ifdef USE_RANDOM
srandom_r((unsigned int)seed, reinterpret_cast<random_data*>(buffer.begin()));
#else
srand48_r(seed, reinterpret_cast<drand48_data*>(buffer.begin()));
#endif
}
Foam::label Foam::osRandomInteger(List<char>& buffer)
{
#ifdef USE_RANDOM
int32_t result;
random_r(reinterpret_cast<random_data*>(buffer.begin()), &result);
return result;
#else
long result;
lrand48_r(reinterpret_cast<drand48_data*>(buffer.begin()), &result);
return result;
#endif
}
Foam::scalar Foam::osRandomDouble(List<char>& buffer)
{
#ifdef USE_RANDOM
int32_t result;
random_r(reinterpret_cast<random_data*>(buffer.begin()), &result);
return (scalar)result/INT_MAX;
#else
double result;
drand48_r(reinterpret_cast<drand48_data*>(buffer.begin()), &result);
return result;
#endif
}
static Foam::DynamicList<Foam::autoPtr<pthread_t>> threads_;
static Foam::DynamicList<Foam::autoPtr<pthread_mutex_t>> mutexes_;
......
......@@ -256,30 +256,6 @@ bool dlSymFound(void* handle, const std::string& symbol);
fileNameList dlLoaded();
// Low level random numbers. Use Random class instead.
//- Seed random number generator.
void osRandomSeed(const label seed);
//- Return random integer (uniform distribution between 0 and 2^31)
label osRandomInteger();
//- Return random double precision (uniform distribution between 0 and 1)
scalar osRandomDouble();
//- Return the size of the buffer for re-entrant random number generator
label osRandomBufferSize();
//- Seed random number generator.
void osRandomSeed(const label seed, List<char>& buffer);
//- Return random integer (uniform distribution between 0 and 2^31)
label osRandomInteger(List<char>& buffer);
//- Return random double precision (uniform distribution between 0 and 1)
scalar osRandomDouble(List<char>& buffer);
// Thread handling
//- Allocate a thread
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2018 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::Rand48
Description
A pseudo random number generator using the linear congruential algorithm
with the following parameters:
\verbatim
Xn+1 = (aXn + c) mod m, where n >= 0
a = 0x5DEECE66D,
c = 0xB
m = 2**48
\endverbatim
It delivers results identical to the \c lrand48() function that is
available on some systems.
\*---------------------------------------------------------------------------*/
#ifndef Rand48_H
#define Rand48_H
#include <cstdint>
#include <random>
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class Rand48 Declaration
\*---------------------------------------------------------------------------*/
class Rand48
{
//- It is a linear_congruential_engine
typedef std::linear_congruential_engine
<
uint64_t,
0x5DEECE66D,
0xB,
(uint64_t(1) << 48)
> engine;
// Private Data
//- The calculation engine
engine engine_;
// Private Member Functions
//- Convert integers
static constexpr uint64_t convert(uint32_t x) noexcept
{
return (static_cast<uint64_t>(x) << 16) | 0x330e;
}
public:
//- The type of the generated random value.
typedef uint32_t result_type;
//- The default seed
static constexpr result_type default_seed = 1u;
// Constructors
//- Construct with specified or default seed
Rand48(uint32_t val = default_seed)
:
engine_(convert(val))
{}
// Member Functions
//- The smallest value that the generator can produce
static constexpr uint32_t min() { return 0; }
//- The largest value that the generator can produce
static constexpr uint32_t max() { return 0x7FFFFFFF; }
//- Reseed the generator
void seed(uint32_t val = default_seed)
{
engine_.seed(convert(val));
}
//- Advance the state of the generator by \c z notches.
void discard(unsigned long long z)
{
engine_.discard(z);
}
//- Get the next random number in the sequence and return its raw
//- 48-bit value
uint64_t raw()
{
return engine_();
}
// Member Operators
//- Get the next random number in the sequence
uint32_t operator()()
{
return static_cast<uint32_t>(engine_() >> 17);
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
#endif
// ************************************************************************* //
......@@ -2,8 +2,8 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2017-2018 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -24,35 +24,25 @@ License
\*---------------------------------------------------------------------------*/
#include "Random.H"
#include "OSspecific.H"
#include "PstreamReduceOps.H"
// * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * //
Foam::scalar Foam::Random::scalar01()
{
return osRandomDouble(buffer_);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::Random::Random(const label seed)
Foam::Random::Random(const label seedValue)
:
buffer_(osRandomBufferSize()),
seed_(seed),
seed_(seedValue),
generator_(seed_),
uniform01_(),
hasGaussSample_(false),
gaussSample_(0)
{
// Initialise the random number generator
osRandomSeed(seed_, buffer_);
}
{}
Foam::Random::Random(const Random& r, const bool reset)
:
buffer_(r.buffer_),
seed_(r.seed_),
generator_(r.generator_),
uniform01_(),
hasGaussSample_(r.hasGaussSample_),
gaussSample_(r.gaussSample_)
{
......@@ -60,34 +50,13 @@ Foam::Random::Random(const Random& r, const bool reset)
{
hasGaussSample_ = false;
gaussSample_ = 0;
// Re-initialise the samples
osRandomSeed(seed_, buffer_);
generator_.seed(seed_);
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::Random::~Random()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::Random::reset(const label seed)
{
seed_ = seed;
osRandomSeed(seed_, buffer_);
}
int Foam::Random::bit()
{
return osRandomInteger(buffer_) & 1;
}
template<>
Foam::scalar Foam::Random::sample01()
{
......@@ -110,23 +79,24 @@ Foam::scalar Foam::Random::GaussNormal()
hasGaussSample_ = false;
return gaussSample_;
}
else
// Gaussian random number as per Knuth/Marsaglia.
// Input: two uniform random numbers, output: two Gaussian random numbers.
// cache one of the values for the next call.
scalar rsq, v1, v2;
do
{
scalar rsq, v1, v2;
do
{
v1 = 2*scalar01() - 1;
v2 = 2*scalar01() - 1;
rsq = sqr(v1) + sqr(v2);
} while (rsq >= 1 || rsq == 0);
v1 = 2*scalar01() - 1;
v2 = 2*scalar01() - 1;
rsq = sqr(v1) + sqr(v2);
} while (rsq >= 1 || rsq == 0);
scalar fac = sqrt(-2*log(rsq)/rsq);
const scalar fac = sqrt(-2*log(rsq)/rsq);
gaussSample_ = v1*fac;
hasGaussSample_ = true;
gaussSample_ = v1*fac;
hasGaussSample_ = true;
return v2*fac;
}
return v2*fac;
}
......@@ -174,16 +144,16 @@ Foam::scalar Foam::Random::globalSample01()
template<>
Foam::label Foam::Random::globalSample01()
{
scalar value = -GREAT;
label value = labelMin;
if (Pstream::master())
{
value = scalar01();
value = round(scalar01());
}
Pstream::scatter(value);
return round(value);
return value;
}
......
......@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox