Use C++11 standard library to generate random numbers

src/FGFDMExec.cpp

@@ -75,7 +75,8 @@ CLASS IMPLEMENTATION
 // Constructor
 
 FGFDMExec::FGFDMExec(FGPropertyManager* root, std::shared_ptr<unsigned int> fdmctr)
-  : RandomEngine(new default_random_engine), FDMctr(fdmctr)
+  : RandomSeed(0), RandomGenerator(make_shared<RandomNumberGenerator>(RandomSeed)),
+    FDMctr(fdmctr)
 {
   Frame           = 0;
   disperse        = 0;
@@ -86,7 +87,6 @@ FGFDMExec::FGFDMExec(FGPropertyManager* root, std::shared_ptr<unsigned int> fdmc
   IsChild = false;
   holding = false;
   Terminate = false;
-  RandomSeed = 0;
   HoldDown = false;
 
   IncrementThenHolding = false;  // increment then hold is off by default
@@ -1244,9 +1244,7 @@ void FGFDMExec::DoTrim(int mode)
 void FGFDMExec::SRand(int sr)
 {
   RandomSeed = sr;
-  gaussian_random_number_phase = 0;
-  RandomEngine->seed(sr);
-  srand(RandomSeed);
+  RandomGenerator->seed(RandomSeed);
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

src/FGFDMExec.h

@@ -42,7 +42,6 @@ INCLUDES
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
 #include <memory>
-#include <random>
 
 #include "models/FGPropagate.h"
 #include "models/FGOutput.h"
@@ -611,7 +610,7 @@ class JSBSIM_API FGFDMExec : public FGJSBBase
     TemplateFunctions[name] = std::make_shared<FGTemplateFunc>(this, el);
   }
 
-  auto GetRandomEngine(void) const { return RandomEngine; }
+  auto GetRandomGenerator(void) const { return RandomGenerator; }
 
 private:
   unsigned int Frame;
@@ -669,8 +668,8 @@ class JSBSIM_API FGFDMExec : public FGJSBBase
 
   bool HoldDown;
 
-  int RandomSeed;
-  std::shared_ptr<std::default_random_engine> RandomEngine;
+  unsigned int RandomSeed;
+  std::shared_ptr<RandomNumberGenerator> RandomGenerator;
 
   // The FDM counter is used to give each child FDM an unique ID. The root FDM
   // has the ID 0

src/FGJSBBase.cpp

@@ -77,8 +77,6 @@ CLASS IMPLEMENTATION
 const string FGJSBBase::needed_cfg_version = "2.0";
 const string FGJSBBase::JSBSim_version = JSBSIM_VERSION " " __DATE__ " " __TIME__ ;
 
-int FGJSBBase::gaussian_random_number_phase = 0;
-
 short FGJSBBase::debug_lvl  = 1;
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -109,34 +107,6 @@ string FGJSBBase::CreateIndexedPropertyName(const string& Property, int index)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-double FGJSBBase::GaussianRandomNumber(void)
-{
-  static double V1, V2, S;
-  double X;
-
-  if (gaussian_random_number_phase == 0) {
-    V1 = V2 = S = X = 0.0;
-
-    do {
-      double U1 = (double)rand() / RAND_MAX;
-      double U2 = (double)rand() / RAND_MAX;
-
-      V1 = 2 * U1 - 1;
-      V2 = 2 * U2 - 1;
-      S = V1 * V1 + V2 * V2;
-    } while(S >= 1 || S == 0);
-
-    X = V1 * sqrt(-2 * log(S) / S);
-  } else
-    X = V2 * sqrt(-2 * log(S) / S);
-
-  gaussian_random_number_phase = 1 - gaussian_random_number_phase;
-
-  return X;
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
 double FGJSBBase::PitotTotalPressure(double mach, double p)
 {
   if (mach < 0) return p;

src/FGJSBBase.h

@@ -43,6 +43,8 @@ INCLUDES
 #include <string>
 #include <cmath>
 #include <stdexcept>
+#include <random>
+#include <chrono>
 
 #include "JSBSim_API.h"
 #include "input_output/string_utilities.h"
@@ -62,6 +64,43 @@ class JSBSIM_API BaseException : public std::runtime_error {
     BaseException(const std::string& msg) : std::runtime_error(msg) {}
 };
 
+/**
+ * @brief Random number generator.
+ * This class encapsulates the C++11 random number generation classes for
+ * uniform and gaussian (aka normal) distributions as well as their seed.
+ * This class guarantees that whenever its seed is reset so are its uniform
+ * and normal random number generators.
+ */
+
+class JSBSIM_API RandomNumberGenerator {
+  public:
+    /// Default constructor using a seed based on the system clock.
+    RandomNumberGenerator(void) : uniform_random(-1.0, 1.0), normal_random(0.0, 1.0)
+    {
+      auto seed_value = std::chrono::system_clock::now().time_since_epoch().count();
+      generator.seed(static_cast<unsigned int>(seed_value));
+    }
+    /// Constructor allowing to specify a seed.
+    RandomNumberGenerator(unsigned int seed)
+      : generator(seed), uniform_random(-1.0, 1.0), normal_random(0.0, 1.0) {}
+    /// Specify a new seed and reinitialize the random generation process.
+    void seed(unsigned int value) {
+      generator.seed(value);
+      uniform_random.reset();
+      normal_random.reset();
+    }
+    /** Get a random number which probability of occurrence is uniformly
+     * distributed over the segment [-1;1( */
+    double GetUniformRandomNumber(void) { return uniform_random(generator); }
+    /** Get a random number which probability of occurrence is following Gauss
+     * normal distribution with a mean of 0.0 and a standard deviation of 1.0 */
+    double GetNormalRandomNumber(void) { return normal_random(generator); }
+  private:
+    std::default_random_engine generator;
+    std::uniform_real_distribution<double> uniform_random;
+    std::normal_distribution<double> normal_random;
+};
+
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 CLASS DOCUMENTATION
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
@@ -288,8 +327,6 @@ class JSBSIM_API FGJSBBase {
 
   static constexpr double sign(double num) {return num>=0.0?1.0:-1.0;}
 
-  static double GaussianRandomNumber(void);
-
 protected:
   static constexpr double radtodeg = 180. / M_PI;
   static constexpr double degtorad = M_PI / 180.;
@@ -318,8 +355,6 @@ class JSBSIM_API FGJSBBase {
 
   static std::string CreateIndexedPropertyName(const std::string& Property, int index);
 
-  static int gaussian_random_number_phase;
-
 public:
 /// Moments L, M, N
 enum {eL     = 1, eM,     eN    };

src/input_output/FGXMLElement.cpp

@@ -670,20 +670,20 @@ double Element::DisperseValue(Element *e, double val, const std::string& supplie
     double disp = e->GetAttributeValueAsNumber("dispersion");
     if (!supplied_units.empty()) disp *= convert[supplied_units][target_units];
     string attType = e->GetAttributeValue("type");
+    RandomNumberGenerator generator;
+
     if (attType == "gaussian" || attType == "gaussiansigned") {
-      double grn = FGJSBBase::GaussianRandomNumber();
-    if (attType == "gaussian") {
-      value = val + disp*grn;
-      } else { // Assume gaussiansigned
-        value = (val + disp*grn)*(fabs(grn)/grn);
-      }
+      double grn = generator.GetNormalRandomNumber();
+      if (attType == "gaussian")
+        value = val + disp*grn;
+      else // Assume gaussiansigned
+        value = (val + disp*grn)*FGJSBBase::sign(grn);
     } else if (attType == "uniform" || attType == "uniformsigned") {
-      double urn = ((((double)rand()/RAND_MAX)-0.5)*2.0);
-      if (attType == "uniform") {
-      value = val + disp * urn;
-      } else { // Assume uniformsigned
-        value = (val + disp * urn)*(fabs(urn)/urn);
-      }
+      double urn = generator.GetUniformRandomNumber();
+      if (attType == "uniform")
+        value = val + disp * urn;
+      else // Assume uniformsigned
+        value = (val + disp * urn)*FGJSBBase::sign(urn);
     } else {
       std::stringstream s;
       s << ReadFrom() << "Unknown dispersion type" << attType;

src/math/FGFunction.cpp

@@ -29,8 +29,6 @@ INCLUDES
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
 #include <iomanip>
-#include <random>
-#include <chrono>
 #include <memory>
 
 #include "simgear/misc/strutils.hxx"
@@ -292,17 +290,17 @@ void FGFunction::CheckOddOrEvenArguments(Element* el, OddEven odd_even)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-shared_ptr<default_random_engine> makeRandomEngine(Element *el, FGFDMExec* fdmex)
+shared_ptr<RandomNumberGenerator> makeRandomGenerator(Element *el, FGFDMExec* fdmex)
 {
   string seed_attr = el->GetAttributeValue("seed");
-  unsigned int seed;
   if (seed_attr.empty())
-    return fdmex->GetRandomEngine();
+    return fdmex->GetRandomGenerator();
   else if (seed_attr == "time_now")
-    seed = chrono::system_clock::now().time_since_epoch().count();
-  else
-    seed = atoi(seed_attr.c_str());
-  return make_shared<default_random_engine>(seed);
+    return make_shared<RandomNumberGenerator>();
+  else {
+    unsigned int seed = atoi(seed_attr.c_str());
+    return make_shared<RandomNumberGenerator>(seed);
+  }
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -606,10 +604,10 @@ void FGFunction::Load(Element* el, FGPropertyValue* var, FGFDMExec* fdmex,
         mean = atof(mean_attr.c_str());
       if (!stddev_attr.empty())
         stddev = atof(stddev_attr.c_str());
-      auto distribution = make_shared<normal_distribution<double>>(mean, stddev);
-      auto generator(makeRandomEngine(element, fdmex));
-      auto f = [generator, distribution]()->double {
-                 return (*distribution.get())(*generator);
+      auto generator(makeRandomGenerator(element, fdmex));
+      auto f = [generator, mean, stddev]()->double {
+                 double value = generator->GetNormalRandomNumber();
+                 return value*stddev + mean;
                };
       Parameters.push_back(new aFunc<decltype(f), 0>(f, PropertyManager, element,
                                                      Prefix));
@@ -622,10 +620,12 @@ void FGFunction::Load(Element* el, FGPropertyValue* var, FGFDMExec* fdmex,
         lower = atof(lower_attr.c_str());
       if (!upper_attr.empty())
         upper = atof(upper_attr.c_str());
-      auto distribution = make_shared<uniform_real_distribution<double>>(lower, upper);
-      auto generator(makeRandomEngine(element, fdmex));
-      auto f = [generator, distribution]()->double {
-                 return (*distribution.get())(*generator);
+      auto generator(makeRandomGenerator(element, fdmex));
+      double a = 0.5*(upper-lower);
+      double b = 0.5*(upper+lower);
+      auto f = [generator, a, b]()->double {
+                 double value = generator->GetUniformRandomNumber();
+                 return value*a + b;
                };
       Parameters.push_back(new aFunc<decltype(f), 0>(f, PropertyManager, element,
                                                      Prefix));

src/models/atmosphere/FGWinds.cpp

@@ -71,7 +71,8 @@ static inline double square_signed (double value)
 /// simply square a value
 constexpr double sqr(double x) { return x*x; }
 
-FGWinds::FGWinds(FGFDMExec* fdmex) : FGModel(fdmex)
+FGWinds::FGWinds(FGFDMExec* fdmex)
+  : FGModel(fdmex), generator(fdmex->GetRandomGenerator())
 {
   Name = "FGWinds";
 
@@ -213,7 +214,7 @@ void FGWinds::Turbulence(double h)
 
     double random = 0.0;
     if (target_time == 0.0) {
-      strength = random = 1 - 2.0*(double(rand())/double(RAND_MAX));
+      strength = random = generator->GetUniformRandomNumber();
       target_time = time + 0.71 + (random * 0.5);
     }
     if (time > target_time) {
@@ -299,10 +300,10 @@ void FGWinds::Turbulence(double h)
       tau_p = L_p/in.V, // eq. (9)
       tau_q = 4*b_w/M_PI/in.V, // eq. (13)
       tau_r =3*b_w/M_PI/in.V, // eq. (17)
-      nu_u = GaussianRandomNumber(),
-      nu_v = GaussianRandomNumber(),
-      nu_w = GaussianRandomNumber(),
-      nu_p = GaussianRandomNumber(),
+      nu_u = generator->GetNormalRandomNumber(),
+      nu_v = generator->GetNormalRandomNumber(),
+      nu_w = generator->GetNormalRandomNumber(),
+      nu_p = generator->GetNormalRandomNumber(),
       xi_u=0, xi_v=0, xi_w=0, xi_p=0, xi_q=0, xi_r=0;
 
     // values of turbulence NED velocities

src/models/atmosphere/FGWinds.h

@@ -394,6 +394,8 @@ class FGWinds : public FGModel {
   FGColumnVector3 vBurstGust;
   FGColumnVector3 vTurbulenceNED;
 
+  std::shared_ptr<RandomNumberGenerator> generator;
+
   void Turbulence(double h);
   void UpDownBurst();
 

src/models/flight_control/FGSensor.cpp

@@ -50,7 +50,8 @@ CLASS IMPLEMENTATION
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
 
-FGSensor::FGSensor(FGFCS* fcs, Element* element) : FGFCSComponent(fcs, element)
+FGSensor::FGSensor(FGFCS* fcs, Element* element)
+  : FGFCSComponent(fcs, element), generator(fcs->GetExec()->GetRandomGenerator())
 {
   // inputs are read from the base class constructor
 
@@ -182,11 +183,10 @@ void FGSensor::Noise(void)
 {
   double random_value=0.0;
 
-  if (DistributionType == eUniform) {
-    random_value = 2.0*(((double)rand()/(double)RAND_MAX) - 0.5);
-  } else {
-    random_value = GaussianRandomNumber();
-  }
+  if (DistributionType == eUniform)
+    random_value = generator->GetUniformRandomNumber();
+  else
+    random_value = generator->GetNormalRandomNumber();
 
   switch( NoiseType ) {
   case ePercent:

src/models/flight_control/FGSensor.h

@@ -178,6 +178,7 @@ class FGSensor  : public FGFCSComponent
   void bind(Element* el, FGPropertyManager* pm) override;
 
 private:
+  std::shared_ptr<RandomNumberGenerator> generator;
   void Debug(int from) override;
 };
 }

tests/unit_tests/FGJSBBaseTest.h

@@ -37,10 +37,6 @@ class FGJSBBaseTest : public CxxTest::TestSuite, public JSBSim::FGJSBBase
     Filter f0;
     Filter f(1.0, 1E-5);
     double x = f.execute(3.0);
-    // Called twice for 100% coverage
-    // Need to test that the numbers follow a Gaussian law ?
-    double ran0 = GaussianRandomNumber();
-    double ran1 = GaussianRandomNumber();
   }
 
   void testTemperatureConversion() {
@@ -68,4 +64,34 @@ class FGJSBBaseTest : public CxxTest::TestSuite, public JSBSim::FGJSBBase
     std::string version = GetVersion();
     disableHighLighting();
   }
+
+  void testRandomNumberGenerator() {
+    JSBSim::RandomNumberGenerator generator(17);
+
+    double u0 = generator.GetUniformRandomNumber();
+    double u1 = generator.GetUniformRandomNumber();
+    double u2 = generator.GetUniformRandomNumber();
+
+    double x0 = generator.GetNormalRandomNumber();
+    double x1 = generator.GetNormalRandomNumber();
+    double x2 = generator.GetNormalRandomNumber();
+
+    // Check that the seed() method correctly resets the random number generator
+    generator.seed(17);
+    double v0 = generator.GetUniformRandomNumber();
+    double v1 = generator.GetUniformRandomNumber();
+    double v2 = generator.GetUniformRandomNumber();
+
+    double y0 = generator.GetNormalRandomNumber();
+    double y1 = generator.GetNormalRandomNumber();
+    double y2 = generator.GetNormalRandomNumber();
+
+    TS_ASSERT_EQUALS(u0, v0);
+    TS_ASSERT_EQUALS(u1, v1);
+    TS_ASSERT_EQUALS(u2, v2);
+
+    TS_ASSERT_EQUALS(x0, y0);
+    TS_ASSERT_EQUALS(x1, y1);
+    TS_ASSERT_EQUALS(x2, y2);
+  }
 };