diff --git a/applications/test/Random/Make/files b/applications/test/Random/Make/files
new file mode 100644
index 0000000000..ab0e0592c0
--- /dev/null
+++ b/applications/test/Random/Make/files
@@ -0,0 +1,3 @@
+Test-Random.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-Random
diff --git a/applications/test/Random/Make/options b/applications/test/Random/Make/options
new file mode 100644
index 0000000000..4c3dd783cb
--- /dev/null
+++ b/applications/test/Random/Make/options
@@ -0,0 +1,3 @@
+EXE_INC =
+
+EXE_LIBS =
diff --git a/applications/test/Random/Test-Random.C b/applications/test/Random/Test-Random.C
new file mode 100644
index 0000000000..b4bd2f2552
--- /dev/null
+++ b/applications/test/Random/Test-Random.C
@@ -0,0 +1,183 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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 .
+
+Application
+ Test-Random
+
+Description
+ Simple test for sequence of random numbers
+
+\*---------------------------------------------------------------------------*/
+
+#include "Rand48.H"
+#include "Random.H"
+
+#include
+#include
+#include
+
+#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
+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 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()*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;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OSspecific/POSIX/POSIX.C b/src/OSspecific/POSIX/POSIX.C
index c37e911aee..b66699bc18 100644
--- a/src/OSspecific/POSIX/POSIX.C
+++ b/src/OSspecific/POSIX/POSIX.C
@@ -66,14 +66,6 @@ Description
#include
#endif
-#ifdef USE_RANDOM
- #include
- #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& buffer)
-{
- #ifdef USE_RANDOM
- srandom_r((unsigned int)seed, reinterpret_cast(buffer.begin()));
- #else
- srand48_r(seed, reinterpret_cast(buffer.begin()));
- #endif
-}
-
-
-Foam::label Foam::osRandomInteger(List& buffer)
-{
- #ifdef USE_RANDOM
- int32_t result;
- random_r(reinterpret_cast(buffer.begin()), &result);
- return result;
- #else
- long result;
- lrand48_r(reinterpret_cast(buffer.begin()), &result);
- return result;
- #endif
-}
-
-
-Foam::scalar Foam::osRandomDouble(List& buffer)
-{
- #ifdef USE_RANDOM
- int32_t result;
- random_r(reinterpret_cast(buffer.begin()), &result);
- return (scalar)result/INT_MAX;
- #else
- double result;
- drand48_r(reinterpret_cast(buffer.begin()), &result);
- return result;
- #endif
-
-}
-
static Foam::DynamicList> threads_;
static Foam::DynamicList> mutexes_;
diff --git a/src/OpenFOAM/include/OSspecific.H b/src/OpenFOAM/include/OSspecific.H
index bdd87752c4..c5077f01a2 100644
--- a/src/OpenFOAM/include/OSspecific.H
+++ b/src/OpenFOAM/include/OSspecific.H
@@ -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& buffer);
-
-//- Return random integer (uniform distribution between 0 and 2^31)
-label osRandomInteger(List& buffer);
-
-//- Return random double precision (uniform distribution between 0 and 1)
-scalar osRandomDouble(List& buffer);
-
-
// Thread handling
//- Allocate a thread
diff --git a/src/OpenFOAM/primitives/random/Random/Rand48.H b/src/OpenFOAM/primitives/random/Random/Rand48.H
new file mode 100644
index 0000000000..fc0a2aec88
--- /dev/null
+++ b/src/OpenFOAM/primitives/random/Random/Rand48.H
@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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 .
+
+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
+#include
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+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(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(engine_() >> 17);
+ }
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+#endif
+
+// ************************************************************************* //
+
diff --git a/src/OpenFOAM/primitives/random/Random/Random.C b/src/OpenFOAM/primitives/random/Random/Random.C
index c17f5c2668..307fcb0eb1 100644
--- a/src/OpenFOAM/primitives/random/Random/Random.C
+++ b/src/OpenFOAM/primitives/random/Random/Random.C
@@ -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;
}
diff --git a/src/OpenFOAM/primitives/random/Random/Random.H b/src/OpenFOAM/primitives/random/Random/Random.H
index fb348fdafa..58d489b3e7 100644
--- a/src/OpenFOAM/primitives/random/Random/Random.H
+++ b/src/OpenFOAM/primitives/random/Random/Random.H
@@ -3,7 +3,7 @@
\\ / 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.
+ \\/ M anipulation | Copyright (C) 2017-2018 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -27,7 +27,6 @@ Class
Description
Random number generator.
-
SourceFiles
RandomI.H
Random.C
@@ -38,15 +37,18 @@ SourceFiles
#ifndef Random_H
#define Random_H
-#include "List.H"
+#include "Rand48.H"
+#include "label.H"
+#include "scalar.H"
+#include
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
-// Forward declaration of classes
-class Random;
+// Forward declarations
+template class UList;
/*---------------------------------------------------------------------------*\
Class Random Declaration
@@ -56,98 +58,99 @@ class Random
{
// Private data
- //- Buffer used by re-entrant random number generator
- List buffer_;
-
//- Initial random number seed
label seed_;
- //- Indicator, which tells if there is a stored gaussian sample
+ //- Random number generator on the int32 interval [0, 2^31)
+ Rand48 generator_;
+
+ //- Uniform distribution on the scalar interval [0, 1]
+ std::uniform_real_distribution uniform01_;
+
+ //- Is there a gaussian sample cached?
bool hasGaussSample_;
- //- Stored sample value
+ //- The cached gaussian sample value
scalar gaussSample_;
// Private Member Functions
- //- Returns the current sample
- scalar scalar01();
+ //- A uniformly distributed floating-point random number [0,1]
+ inline scalar scalar01();
public:
// Constructors
- //- Construct given seed and sample count
- Random(const label seed = 123456);
+ //- Construct with seed value
+ Random(const label seedValue = 123456);
- //- Copy constructor with optional reset of sampleI
+ //- Copy construct with optional reset of seed
Random(const Random& r, const bool reset = false);
// Destructor
- ~Random();
+ ~Random() = default;
- // Member functions
+ // Member Functions
- // Access
+ // Access
- //- Return const access to the initial random number seed
- inline label seed() const;
+ //- The initial random number seed that was used
+ inline label seed() const;
- //- Reset the random number generator seed
- void reset(const label seed);
+ //- Reset the random number generator seed.
+ inline void reset(const label seedValue);
- // Evaluation
+ // Random numbers
- // Random numbers
+ //- Return a random bit
+ inline int bit();
- //- Return a random bit
- int bit();
+ //- Return a sample whose components lie in the range 0-1
+ template
+ Type sample01();
- //- Return a sample whose components lie in the range 0-1
- template
- Type sample01();
+ //- Return a sample whose components are normally distributed
+ // with zero mean and unity variance N(0, 1)
+ template
+ Type GaussNormal();
- //- Return a sample whose components are normally distributed
- // with zero mean and unity variance N(0, 1)
- template
- Type GaussNormal();
+ //- Return a sample between start and end
+ template
+ Type position(const Type& start, const Type& end);
- //- Return a sample between start and end
- template
- Type position(const Type& start, const Type& end);
+ //- Randomise value in the range 0-1
+ template
+ void randomise01(Type& value);
- //- Randomise value in the range 0-1
- template
- void randomise01(Type& value);
-
- //- Shuffle the values in the list
- template
- void shuffle(UList& values);
+ //- Shuffle the values in the list
+ template
+ void shuffle(UList& values);
- // Global random numbers - consistent across all processors
+ // Global random numbers - consistent across all processors
- //- Return a sample whose components lie in the range 0-1
- template
- Type globalSample01();
+ //- Return a sample whose components lie in the range 0-1
+ template
+ Type globalSample01();
- //- Return a sample whose components are normally distributed
- // with zero mean and unity variance N(0, 1)
- template
- Type globalGaussNormal();
+ //- Return a sample whose components are normally distributed
+ // with zero mean and unity variance N(0, 1)
+ template
+ Type globalGaussNormal();
- //- Return a sample between start and end
- template
- Type globalPosition(const Type& start, const Type& end);
+ //- Return a sample between start and end
+ template
+ Type globalPosition(const Type& start, const Type& end);
- //- Randomise value in the range 0-1
- template
- void globalRandomise01(Type& value);
+ //- Randomise value in the range 0-1
+ template
+ void globalRandomise01(Type& value);
};
@@ -167,11 +170,7 @@ template<>
label Random::GaussNormal