--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/core/model/random-variable-stream.cc Tue Jul 10 21:31:47 2012 -0700
@@ -0,0 +1,189 @@
+/* -*- Mode:C++; c-file-style:"gnu"; indent-tabs-mode:nil; -*- */
+/*
+ * Copyright (c) 2006 Georgia Tech Research Corporation
+ * Copyright (c) 2011 Mathieu Lacage
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License version 2 as
+ * published by the Free Software Foundation;
+ *
+ * This program 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 this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ * Authors: Rajib Bhattacharjea<raj.b@gatech.edu>
+ * Hadi Arbabi<marbabi@cs.odu.edu>
+ * Mathieu Lacage <mathieu.lacage@gmail.com>
+ *
+ * Modified by Mitch Watrous <watrous@u.washington.edu>
+ *
+ */
+#include "random-variable-stream.h"
+#include "assert.h"
+#include "boolean.h"
+#include "double.h"
+#include "integer.h"
+#include "string.h"
+#include "pointer.h"
+#include "log.h"
+#include "rng-stream.h"
+#include "rng-seed-manager.h"
+#include <cmath>
+#include <iostream>
+
+NS_LOG_COMPONENT_DEFINE ("RandomVariableStream");
+
+namespace ns3 {
+
+NS_OBJECT_ENSURE_REGISTERED (RandomVariableStream);
+
+TypeId
+RandomVariableStream::GetTypeId (void)
+{
+ static TypeId tid = TypeId ("ns3::RandomVariableStream")
+ .SetParent<Object> ()
+ .AddAttribute("Stream",
+ "The stream number for this RNG stream. -1 means \"allocate a stream automatically\". "
+ "Note that if -1 is set, Get will return -1 so that it is not possible to know which "
+ "value was automatically allocated.",
+ IntegerValue(-1),
+ MakeIntegerAccessor(&RandomVariableStream::SetStream,
+ &RandomVariableStream::GetStream),
+ MakeIntegerChecker<int64_t>())
+ .AddAttribute("Antithetic", "Set this RNG stream to generate antithetic values",
+ BooleanValue (false),
+ MakeBooleanAccessor(&RandomVariableStream::SetAntithetic,
+ &RandomVariableStream::IsAntithetic),
+ MakeBooleanChecker())
+ ;
+ return tid;
+}
+
+RandomVariableStream::RandomVariableStream()
+ : m_rng (0)
+{}
+RandomVariableStream::~RandomVariableStream()
+{
+ delete m_rng;
+}
+
+void
+RandomVariableStream::SetAntithetic(bool isAntithetic)
+{
+ m_isAntithetic = isAntithetic;
+}
+bool
+RandomVariableStream::IsAntithetic(void) const
+{
+ return m_isAntithetic;
+}
+void
+RandomVariableStream::SetStream (int64_t stream)
+{
+ NS_LOG_FUNCTION (this << stream);
+ // negative values are not legal.
+ NS_ASSERT (stream >= -1);
+ delete m_rng;
+ if (stream == -1)
+ {
+ // The first 2^63 streams are reserved for automatic stream
+ // number assignment.
+ uint64_t nextStream = RngSeedManager::GetNextStreamIndex ();
+ NS_ASSERT(nextStream <= ((1ULL)<<63));
+ m_rng = new RngStream (RngSeedManager::GetSeed (),
+ nextStream,
+ RngSeedManager::GetRun ());
+ }
+ else
+ {
+ // The last 2^63 streams are reserved for deterministic stream
+ // number assignment.
+ uint64_t base = ((1ULL)<<63);
+ uint64_t target = base + stream;
+ m_rng = new RngStream (RngSeedManager::GetSeed (),
+ target,
+ RngSeedManager::GetRun ());
+ }
+ m_stream = stream;
+}
+int64_t
+RandomVariableStream::GetStream(void) const
+{
+ return m_stream;
+}
+
+RngStream *
+RandomVariableStream::Peek(void) const
+{
+ return m_rng;
+}
+
+NS_OBJECT_ENSURE_REGISTERED(UniformRandomVariable);
+
+TypeId
+UniformRandomVariable::GetTypeId (void)
+{
+ static TypeId tid = TypeId ("ns3::UniformRandomVariable")
+ .SetParent<RandomVariableStream>()
+ .AddConstructor<UniformRandomVariable> ()
+ .AddAttribute("Min", "The lower bound on the values returned by this RNG stream.",
+ DoubleValue(0),
+ MakeDoubleAccessor(&UniformRandomVariable::m_min),
+ MakeDoubleChecker<double>())
+ .AddAttribute("Max", "The upper bound on the values returned by this RNG stream.",
+ DoubleValue(1.0),
+ MakeDoubleAccessor(&UniformRandomVariable::m_max),
+ MakeDoubleChecker<double>())
+ ;
+ return tid;
+}
+UniformRandomVariable::UniformRandomVariable ()
+{
+ // m_min and m_max are initialized after constructor by attributes
+}
+
+double
+UniformRandomVariable::GetMin (void) const
+{
+ return m_min;
+}
+double
+UniformRandomVariable::GetMax (void) const
+{
+ return m_max;
+}
+
+double
+UniformRandomVariable::GetValue (double min, double max)
+{
+ double v = min + Peek ()->RandU01 () * (max - min);
+ if (IsAntithetic ())
+ {
+ v = min + (max - v);
+ }
+ return v;
+}
+uint32_t
+UniformRandomVariable::GetInteger (uint32_t min, uint32_t max)
+{
+ NS_ASSERT (min <= max);
+ return static_cast<uint32_t> ( GetValue (min, max + 1) );
+}
+
+double
+UniformRandomVariable::GetValue (void)
+{
+ return GetValue (m_min, m_max);
+}
+uint32_t
+UniformRandomVariable::GetInteger (void)
+{
+ return (uint32_t)GetValue (m_min, m_max + 1);
+}
+
+} // namespace ns3
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/core/model/random-variable-stream.h Tue Jul 10 21:31:47 2012 -0700
@@ -0,0 +1,269 @@
+/* -*- Mode:C++; c-file-style:"gnu"; indent-tabs-mode:nil; -*- */
+/*
+ * Copyright (c) 2006 Georgia Tech Research Corporation
+ * Copyright (c) 2011 Mathieu Lacage
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License version 2 as
+ * published by the Free Software Foundation;
+ *
+ * This program 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 this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ * Authors: Rajib Bhattacharjea<raj.b@gatech.edu>
+ * Hadi Arbabi<marbabi@cs.odu.edu>
+ * Mathieu Lacage <mathieu.lacage@gmail.com>
+ *
+ * Modified by Mitch Watrous <watrous@u.washington.edu>
+ *
+ */
+#ifndef RANDOM_VARIABLE_STREAM_H
+#define RANDOM_VARIABLE_STREAM_H
+
+#include "type-id.h"
+#include "object.h"
+#include "attribute-helper.h"
+#include <stdint.h>
+
+namespace ns3 {
+
+class RngStream;
+
+/**
+ * \ingroup randomvariable
+ * \brief The Random Number Generator (RNG) that allows stream numbers to be set deterministically.
+ *
+ * Note: The underlying random number generation method used
+ * by NS-3 is the RngStream code by Pierre L'Ecuyer at
+ * the University of Montreal.
+ *
+ * NS-3 has a rich set of random number generators that allow stream
+ * numbers to be set deterministically if desired. Class
+ * RandomVariableStream defines the base class functionalty required
+ * for all such random number generators.
+ *
+ * By default, the underlying generator is seeded all the time with
+ * the same seed value and run number coming from the ns3::GlobalValue
+ * \ref GlobalValueRngSeed "RngSeed" and \ref GlobalValueRngRun
+ * "RngRun". Also by default, the stream number value for this RNG
+ * stream is automatically allocated.
+ */
+class RandomVariableStream : public Object
+{
+public:
+ static TypeId GetTypeId (void);
+ RandomVariableStream ();
+ virtual ~RandomVariableStream();
+
+ /**
+ * \brief Specifies the stream number for this RNG stream.
+ * \param stream The stream number for this RNG stream. -1 means "allocate a stream number automatically".
+ */
+ void SetStream (int64_t stream);
+
+ /**
+ * \brief Returns the stream number for this RNG stream.
+ * \return The stream number for this RNG stream. -1 means "allocate a stream number automatically".
+ */
+ int64_t GetStream(void) const;
+
+ /**
+ * \brief Specifies whether antithetic values should be generated.
+ * \param isAntithetic Set equal to true if antithetic values should
+ * be generated by this RNG stream.
+ */
+ void SetAntithetic(bool isAntithetic);
+
+ /**
+ * \brief Returns true if antithetic values should be generated.
+ * \return A bool value that indicates if antithetic values should
+ * be generated by this RNG stream.
+ */
+ bool IsAntithetic(void) const;
+
+ /**
+ * \brief Returns a random double from the underlying distribution
+ * \return A floating point random value.
+ */
+ virtual double GetValue (void) = 0;
+
+ /**
+ * \brief Returns a random integer integer from the underlying distribution
+ * \return Integer cast of RandomVariableStream::GetValue
+ */
+ virtual uint32_t GetInteger (void) = 0;
+
+protected:
+ /**
+ * \brief Returns a pointer to the underlying RNG stream.
+ */
+ RngStream *Peek(void) const;
+
+private:
+ // you can't copy these objects.
+ // Theoretically, it is possible to give them good copy semantics
+ // but not enough time to iron out the details.
+ RandomVariableStream (const RandomVariableStream &o);
+ RandomVariableStream &operator = (const RandomVariableStream &o);
+
+ /// Pointer to the underlying RNG stream.
+ RngStream *m_rng;
+
+ /// Indicates if antithetic values should be generated by this RNG stream.
+ bool m_isAntithetic;
+
+ /// The stream number for this RNG stream.
+ int64_t m_stream;
+};
+
+/**
+ * \ingroup randomvariable
+ * \brief The uniform distribution Random Number Generator (RNG) that allows stream numbers to be set deterministically.
+ *
+ * This class supports the creation of objects that return random numbers
+ * from a fixed uniform distribution. It also supports the generation of
+ * single random numbers from various uniform distributions.
+ *
+ * The low end of the range is always included and the high end
+ * of the range is always excluded.
+ *
+ * Here is an example of how to use this class:
+ * \code
+ * double min = 0.0;
+ * double max = 10.0;
+ *
+ * Ptr<UniformRandomVariable> x = CreateObject<UniformRandomVariable> ();
+ * x->SetAttribute ("Min", DoubleValue (min));
+ * x->SetAttribute ("Max", DoubleValue (max));
+ *
+ * // The values returned by a uniformly distributed random
+ * // variable should always be within the range
+ * //
+ * // [min, max) .
+ * //
+ * double value = x->GetValue ();
+ * \endcode
+ */
+class UniformRandomVariable : public RandomVariableStream
+{
+public:
+ static TypeId GetTypeId (void);
+
+ /**
+ * \brief Creates a uniform distribution RNG with the default range.
+ */
+ UniformRandomVariable ();
+
+ /**
+ * \brief Returns the lower bound on values that can be returned by this RNG stream.
+ * \return The lower bound on values that can be returned by this RNG stream.
+ */
+ double GetMin (void) const;
+
+ /**
+ * \brief Returns the upper bound on values that can be returned by this RNG stream.
+ * \return The upper bound on values that can be returned by this RNG stream.
+ */
+ double GetMax (void) const;
+
+ /**
+ * \brief Returns a random double from the uniform distribution with the specified range.
+ * \param min Low end of the range.
+ * \param max High end of the range.
+ * \return A floating point random value.
+ *
+ * Note that antithetic values are being generated if
+ * m_isAntithetic is equal to true. If \f$x\f$ is a value that
+ * would be returned normally, then \f$(max - x\f$) is the distance
+ * that \f$x\f$ would be from \f$max\f$. The value returned in the
+ * antithetic case, \f$x'\f$, is calculated as
+ *
+ * \f[
+ * x' = min + (max - x) ,
+ * \f]
+ *
+ * which is the lower bound plus the distance \f$x\f$ is from the
+ * upper bound.
+ */
+ double GetValue (double min, double max);
+
+ /**
+ * \brief Returns a random unsigned integer from a uniform distribution over the interval [min,max] including both ends.
+ * \param min Low end of the range.
+ * \param max High end of the range.
+ * \return A random unsigned integer value.
+ *
+ * Note that antithetic values are being generated if
+ * m_isAntithetic is equal to true. If \f$x\f$ is a value that
+ * would be returned normally, then \f$(max - x\f$) is the distance
+ * that \f$x\f$ would be from \f$max\f$. The value returned in the
+ * antithetic case, \f$x'\f$, is calculated as
+ *
+ * \f[
+ * x' = min + (max - x) ,
+ * \f]
+ *
+ * which is the lower bound plus the distance \f$x\f$ is from the
+ * upper bound.
+ */
+ uint32_t GetInteger (uint32_t min, uint32_t max);
+
+ /**
+ * \brief Returns a random double from the uniform distribution with the range [min,max), where min and max are the current lower and upper bounds.
+ * \return A floating point random value.
+ *
+ * Note that antithetic values are being generated if
+ * m_isAntithetic is equal to true. If \f$x\f$ is a value that
+ * would be returned normally, then \f$(max - x\f$) is the distance
+ * that \f$x\f$ would be from \f$max\f$. The value returned in the
+ * antithetic case, \f$x'\f$, is calculated as
+ *
+ * \f[
+ * x' = min + (max - x) ,
+ * \f]
+ *
+ * which is the lower bound plus the distance \f$x\f$ is from the
+ * upper bound.
+ *
+ * Note that we have to re-implement this method here because the method is
+ * overloaded above for the two-argument variant and the c++ name resolution
+ * rules don't work well with overloads split between parent and child
+ * classes.
+ */
+ virtual double GetValue (void);
+
+ /**
+ * \brief Returns a random unsigned integer from a uniform distribution over the interval [min,max] including both ends, where min and max are the current lower and upper bounds.
+ * \return A random unsigned integer value.
+ *
+ * Note that antithetic values are being generated if
+ * m_isAntithetic is equal to true. If \f$x\f$ is a value that
+ * would be returned normally, then \f$(max - x\f$) is the distance
+ * that \f$x\f$ would be from \f$max\f$. The value returned in the
+ * antithetic case, \f$x'\f$, is calculated as
+ *
+ * \f[
+ * x' = min + (max - x) ,
+ * \f]
+ *
+ * which is the lower bound plus the distance \f$x\f$ is from the
+ * upper bound.
+ */
+ virtual uint32_t GetInteger (void);
+private:
+ /// The lower bound on values that can be returned by this RNG stream.
+ double m_min;
+
+ /// The upper bound on values that can be returned by this RNG stream.
+ double m_max;
+};
+
+} // namespace ns3
+
+#endif /* RANDOM_VARIABLE_STREAM_H */
--- a/src/core/model/random-variable.cc Mon Jul 09 06:53:14 2012 -0700
+++ b/src/core/model/random-variable.cc Tue Jul 10 21:31:47 2012 -0700
@@ -33,9 +33,8 @@
#include <vector>
#include "assert.h"
-#include "config.h"
-#include "integer.h"
#include "random-variable.h"
+#include "rng-seed-manager.h"
#include "rng-stream.h"
#include "fatal-error.h"
@@ -44,44 +43,6 @@
namespace ns3 {
// -----------------------------------------------------------------------------
-// Seed Manager
-// -----------------------------------------------------------------------------
-
-uint32_t SeedManager::GetSeed ()
-{
- uint32_t s[6];
- RngStream::GetPackageSeed (s);
- NS_ASSERT (
- s[0] == s[1]
- && s[0] == s[2]
- && s[0] == s[3]
- && s[0] == s[4]
- && s[0] == s[5]
- );
- return s[0];
-}
-
-void SeedManager::SetSeed (uint32_t seed)
-{
- Config::SetGlobal ("RngSeed", IntegerValue (seed));
-}
-
-void SeedManager::SetRun (uint32_t run)
-{
- Config::SetGlobal ("RngRun", IntegerValue (run));
-}
-
-uint32_t SeedManager::GetRun ()
-{
- return RngStream::GetPackageRun ();
-}
-
-bool SeedManager::CheckSeed (uint32_t seed)
-{
- return RngStream::CheckSeed (seed);
-}
-
-// -----------------------------------------------------------------------------
// -----------------------------------------------------------------------------
// RandomVariableBase methods
@@ -95,22 +56,24 @@
virtual double GetValue () = 0;
virtual uint32_t GetInteger ();
virtual RandomVariableBase* Copy (void) const = 0;
-
-protected:
+ RngStream *GetStream(void);
+private:
RngStream* m_generator; // underlying generator being wrapped
};
RandomVariableBase::RandomVariableBase ()
- : m_generator (NULL)
+ : m_generator (0)
{
}
RandomVariableBase::RandomVariableBase (const RandomVariableBase& r)
: m_generator (0)
{
- if (r.m_generator)
+ if (r.m_generator != 0)
{
- m_generator = new RngStream (*r.m_generator);
+ m_generator = new RngStream (RngSeedManager::GetSeed (),
+ RngSeedManager::GetNextStreamIndex (),
+ RngSeedManager::GetRun ());
}
}
@@ -124,6 +87,18 @@
return (uint32_t)GetValue ();
}
+RngStream *
+RandomVariableBase::GetStream (void)
+{
+ if (m_generator == 0)
+ {
+ m_generator = new RngStream (RngSeedManager::GetSeed (),
+ RngSeedManager::GetNextStreamIndex (),
+ RngSeedManager::GetRun ());
+ }
+ return m_generator;
+}
+
// -------------------------------------------------------
RandomVariable::RandomVariable ()
@@ -250,20 +225,14 @@
double UniformVariableImpl::GetValue ()
{
- if (!m_generator)
- {
- m_generator = new RngStream ();
- }
- return m_min + m_generator->RandU01 () * (m_max - m_min);
+ RngStream *generator = GetStream ();
+ return m_min + generator->RandU01 () * (m_max - m_min);
}
double UniformVariableImpl::GetValue (double s, double l)
{
- if (!m_generator)
- {
- m_generator = new RngStream ();
- }
- return s + m_generator->RandU01 () * (l - s);
+ RngStream *generator = GetStream ();
+ return s + generator->RandU01 () * (l - s);
}
RandomVariableBase* UniformVariableImpl::Copy () const
@@ -573,13 +542,10 @@
double ExponentialVariableImpl::GetValue ()
{
- if (!m_generator)
- {
- m_generator = new RngStream ();
- }
+ RngStream *generator = GetStream ();
while (1)
{
- double r = -m_mean*log (m_generator->RandU01 ());
+ double r = -m_mean*log (generator->RandU01 ());
if (m_bound == 0 || r <= m_bound)
{
return r;
@@ -737,13 +703,10 @@
double ParetoVariableImpl::GetValue ()
{
- if (!m_generator)
- {
- m_generator = new RngStream ();
- }
+ RngStream *generator = GetStream ();
while (1)
{
- double r = (m_scale * ( 1.0 / pow (m_generator->RandU01 (), 1.0 / m_shape)));
+ double r = (m_scale * ( 1.0 / pow (generator->RandU01 (), 1.0 / m_shape)));
if (m_bound == 0 || r <= m_bound)
{
return r;
@@ -871,14 +834,11 @@
double WeibullVariableImpl::GetValue ()
{
- if (!m_generator)
- {
- m_generator = new RngStream ();
- }
+ RngStream *generator = GetStream ();
double exponent = 1.0 / m_alpha;
while (1)
{
- double r = m_mean * pow ( -log (m_generator->RandU01 ()), exponent);
+ double r = m_mean * pow ( -log (generator->RandU01 ()), exponent);
if (m_bound == 0 || r <= m_bound)
{
return r;
@@ -981,10 +941,7 @@
double NormalVariableImpl::GetValue ()
{
- if (!m_generator)
- {
- m_generator = new RngStream ();
- }
+ RngStream *generator = GetStream ();
if (m_nextValid)
{ // use previously generated
m_nextValid = false;
@@ -994,8 +951,8 @@
{ // See Simulation Modeling and Analysis p. 466 (Averill Law)
// for algorithm; basically a Box-Muller transform:
// http://en.wikipedia.org/wiki/Box-Muller_transform
- double u1 = m_generator->RandU01 ();
- double u2 = m_generator->RandU01 ();
+ double u1 = generator->RandU01 ();
+ double u2 = generator->RandU01 ();
double v1 = 2 * u1 - 1;
double v2 = 2 * u2 - 1;
double w = v1 * v1 + v2 * v2;
@@ -1138,10 +1095,7 @@
double EmpiricalVariableImpl::GetValue ()
{ // Return a value from the empirical distribution
// This code based (loosely) on code by Bruce Mah (Thanks Bruce!)
- if (!m_generator)
- {
- m_generator = new RngStream ();
- }
+ RngStream *generator = GetStream ();
if (emp.size () == 0)
{
return 0.0; // HuH? No empirical data
@@ -1150,7 +1104,7 @@
{
Validate (); // Insure in non-decreasing
}
- double r = m_generator->RandU01 ();
+ double r = generator->RandU01 ();
if (r <= emp.front ().cdf)
{
return emp.front ().value; // Less than first
@@ -1406,18 +1360,15 @@
double
LogNormalVariableImpl::GetValue ()
{
- if (!m_generator)
- {
- m_generator = new RngStream ();
- }
+ RngStream *generator = GetStream ();
double u, v, r2, normal, z;
do
{
/* choose x,y in uniform square (-1,-1) to (+1,+1) */
- u = -1 + 2 * m_generator->RandU01 ();
- v = -1 + 2 * m_generator->RandU01 ();
+ u = -1 + 2 * generator->RandU01 ();
+ v = -1 + 2 * generator->RandU01 ();
/* see if it is in the unit circle */
r2 = u * u + v * v;
@@ -1504,14 +1455,11 @@
double
GammaVariableImpl::GetValue (double alpha, double beta)
{
- if (!m_generator)
- {
- m_generator = new RngStream ();
- }
+ RngStream *generator = GetStream ();
if (alpha < 1)
{
- double u = m_generator->RandU01 ();
+ double u = generator->RandU01 ();
return GetValue (1.0 + alpha, beta) * pow (u, 1.0 / alpha);
}
@@ -1529,7 +1477,7 @@
while (v <= 0);
v = v * v * v;
- u = m_generator->RandU01 ();
+ u = generator->RandU01 ();
if (u < 1 - 0.0331 * x * x * x * x)
{
break;
@@ -1627,11 +1575,8 @@
double
ErlangVariableImpl::GetValue (unsigned int k, double lambda)
{
- if (!m_generator)
- {
- m_generator = new RngStream ();
- }
-
+ // XXX: Fixme: do not create a new
+ // RNG stream every time the function is called !
ExponentialVariable exponential (lambda);
double result = 0;
@@ -1723,11 +1668,8 @@
double TriangularVariableImpl::GetValue ()
{
- if (!m_generator)
- {
- m_generator = new RngStream ();
- }
- double u = m_generator->RandU01 ();
+ RngStream *generator = GetStream ();
+ double u = generator->RandU01 ();
if (u <= (m_mode - m_min) / (m_max - m_min) )
{
return m_min + sqrt (u * (m_max - m_min) * (m_mode - m_min) );
@@ -1811,12 +1753,9 @@
double
ZipfVariableImpl::GetValue ()
{
- if (!m_generator)
- {
- m_generator = new RngStream ();
- }
+ RngStream *generator = GetStream ();
- double u = m_generator->RandU01 ();
+ double u = generator->RandU01 ();
double sum_prob = 0,zipf_value = 0;
for (int i = 1; i <= m_n; i++)
{
@@ -1895,10 +1834,7 @@
double
ZetaVariableImpl::GetValue ()
{
- if (!m_generator)
- {
- m_generator = new RngStream ();
- }
+ RngStream *generator = GetStream ();
double u, v;
double X, T;
@@ -1906,8 +1842,8 @@
do
{
- u = m_generator->RandU01 ();
- v = m_generator->RandU01 ();
+ u = generator->RandU01 ();
+ v = generator->RandU01 ();
X = floor (pow (u, -1.0 / (m_alpha - 1.0)));
T = pow (1.0 + 1.0 / X, m_alpha - 1.0);
test = v * X * (T - 1.0) / (m_b - 1.0);
--- a/src/core/model/random-variable.h Mon Jul 09 06:53:14 2012 -0700
+++ b/src/core/model/random-variable.h Tue Jul 10 21:31:47 2012 -0700
@@ -29,6 +29,7 @@
#include <ostream>
#include "attribute.h"
#include "attribute-helper.h"
+#include "rng-seed-manager.h"
/**
* \ingroup core
@@ -40,69 +41,6 @@
class RandomVariableBase;
-class SeedManager
-{
-public:
- /**
- * \brief set the seed
- * it will duplicate the seed value 6 times
- * \code
- * SeedManger::SetSeed(15);
- * UniformVariable x(2,3); //these will give the same output everytime
- * ExponentialVariable y(120); //as long as the seed stays the same
- * \endcode
- * \param seed
- *
- * Note, while the underlying RNG takes six integer values as a seed;
- * it is sufficient to set these all to the same integer, so we provide
- * a simpler interface here that just takes one integer.
- */
- static void SetSeed (uint32_t seed);
-
- /**
- * \brief Get the seed value
- * \return the seed value
- *
- * Note: returns the first of the six seed values used in the underlying RNG
- */
- static uint32_t GetSeed ();
-
- /**
- * \brief Set the run number of simulation
- *
- * \code
- * SeedManager::SetSeed(12);
- * int N = atol(argv[1]); //read in run number from command line
- * SeedManager::SetRun(N);
- * UniformVariable x(0,10);
- * ExponentialVariable y(2902);
- * \endcode
- * In this example, N could successivly be equal to 1,2,3, etc. and the user
- * would continue to get independent runs out of the single simulation. For
- * this simple example, the following might work:
- * \code
- * ./simulation 0
- * ...Results for run 0:...
- *
- * ./simulation 1
- * ...Results for run 1:...
- * \endcode
- */
- static void SetRun (uint32_t run);
- /**
- * \returns the current run number
- * @sa SetRun
- */
- static uint32_t GetRun (void);
-
- /**
- * \brief Check if seed value is valid if wanted to be used as seed
- * \return true if valid and false if invalid
- */
- static bool CheckSeed (uint32_t seed);
-};
-
-
/**
* \brief The basic RNG for NS-3.
* \ingroup randomvariable
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/core/model/rng-seed-manager.cc Tue Jul 10 21:31:47 2012 -0700
@@ -0,0 +1,52 @@
+#include "rng-seed-manager.h"
+#include "global-value.h"
+#include "attribute-helper.h"
+#include "integer.h"
+#include "config.h"
+
+namespace ns3 {
+
+static uint64_t g_nextStreamIndex = 0;
+static ns3::GlobalValue g_rngSeed ("RngSeed",
+ "The global seed of all rng streams",
+ ns3::IntegerValue(1),
+ ns3::MakeIntegerChecker<uint32_t> ());
+static ns3::GlobalValue g_rngRun ("RngRun",
+ "The run number used to modify the global seed",
+ ns3::IntegerValue (1),
+ ns3::MakeIntegerChecker<int64_t> ());
+
+
+uint32_t RngSeedManager::GetSeed (void)
+{
+ IntegerValue seedValue;
+ g_rngSeed.GetValue (seedValue);
+ return seedValue.Get ();
+}
+void
+RngSeedManager::SetSeed (uint32_t seed)
+{
+ Config::SetGlobal ("RngSeed", IntegerValue(seed));
+}
+
+void RngSeedManager::SetRun (uint64_t run)
+{
+ Config::SetGlobal ("RngRun", IntegerValue (run));
+}
+
+uint64_t RngSeedManager::GetRun ()
+{
+ IntegerValue value;
+ g_rngRun.GetValue (value);
+ int run = value.Get();
+ return run;
+}
+
+uint64_t RngSeedManager::GetNextStreamIndex (void)
+{
+ uint64_t next = g_nextStreamIndex;
+ g_nextStreamIndex++;
+ return next;
+}
+
+} // namespace ns3
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/core/model/rng-seed-manager.h Tue Jul 10 21:31:47 2012 -0700
@@ -0,0 +1,73 @@
+#ifndef RNG_SEED_MANAGER_H
+#define RNG_SEED_MANAGER_H
+
+#include <stdint.h>
+
+namespace ns3 {
+
+class RngSeedManager
+{
+public:
+ /**
+ * \brief set the seed
+ * it will duplicate the seed value 6 times
+ * \code
+ * SeedManger::SetSeed(15);
+ * UniformVariable x(2,3); //these will give the same output everytime
+ * ExponentialVariable y(120); //as long as the seed stays the same
+ * \endcode
+ * \param seed
+ *
+ * Note, while the underlying RNG takes six integer values as a seed;
+ * it is sufficient to set these all to the same integer, so we provide
+ * a simpler interface here that just takes one integer.
+ */
+ static void SetSeed (uint32_t seed);
+
+ /**
+ * \brief Get the seed value
+ * \return the seed value
+ *
+ * Note: returns the first of the six seed values used in the underlying RNG
+ */
+ static uint32_t GetSeed (void);
+
+ /**
+ * \brief Set the run number of simulation
+ *
+ * \code
+ * SeedManager::SetSeed(12);
+ * int N = atol(argv[1]); //read in run number from command line
+ * SeedManager::SetRun(N);
+ * UniformVariable x(0,10);
+ * ExponentialVariable y(2902);
+ * \endcode
+ * In this example, N could successivly be equal to 1,2,3, etc. and the user
+ * would continue to get independent runs out of the single simulation. For
+ * this simple example, the following might work:
+ * \code
+ * ./simulation 0
+ * ...Results for run 0:...
+ *
+ * ./simulation 1
+ * ...Results for run 1:...
+ * \endcode
+ */
+ static void SetRun (uint64_t run);
+ /**
+ * \returns the current run number
+ * @sa SetRun
+ */
+ static uint64_t GetRun (void);
+
+ static uint64_t GetNextStreamIndex(void);
+
+};
+
+// for compatibility
+typedef RngSeedManager SeedManager;
+
+} // namespace ns3
+
+
+#endif /* RNG_SEED_MANAGER_H */
--- a/src/core/model/rng-stream.cc Mon Jul 09 06:53:14 2012 -0700
+++ b/src/core/model/rng-stream.cc Tue Jul 10 21:31:47 2012 -0700
@@ -15,18 +15,20 @@
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
-// Modified for ns-3 by: Rajib Bhattacharjea<raj.b@gatech.edu>
+// Modified for ns-3 by:
+// - Rajib Bhattacharjea<raj.b@gatech.edu>
+// - Mathieu Lacage <mathieu.lacage@gmail.com>
//
#include <cstdlib>
#include <iostream>
#include "rng-stream.h"
-#include "global-value.h"
-#include "integer.h"
-using namespace std;
+#include "fatal-error.h"
namespace
{
+typedef double Matrix[3][3];
+
const double m1 = 4294967087.0;
const double m2 = 4294944443.0;
const double norm = 1.0 / (m1 + 1.0);
@@ -38,58 +40,30 @@
const double two53 = 9007199254740992.0;
const double fact = 5.9604644775390625e-8; /* 1 / 2^24 */
-// The following are the transition matrices of the two MRG components
-// (in matrix form), raised to the powers -1, 1, 2^76, and 2^127, resp.
-
-const double InvA1[3][3] = { // Inverse of A1p0
+const Matrix InvA1 = { // Inverse of A1p0
{ 184888585.0, 0.0, 1945170933.0 },
{ 1.0, 0.0, 0.0 },
{ 0.0, 1.0, 0.0 }
};
-const double InvA2[3][3] = { // Inverse of A2p0
+const Matrix InvA2 = { // Inverse of A2p0
{ 0.0, 360363334.0, 4225571728.0 },
{ 1.0, 0.0, 0.0 },
{ 0.0, 1.0, 0.0 }
};
-const double A1p0[3][3] = {
+const Matrix A1p0 = {
{ 0.0, 1.0, 0.0 },
{ 0.0, 0.0, 1.0 },
{ -810728.0, 1403580.0, 0.0 }
};
-const double A2p0[3][3] = {
+const Matrix A2p0 = {
{ 0.0, 1.0, 0.0 },
{ 0.0, 0.0, 1.0 },
{ -1370589.0, 0.0, 527612.0 }
};
-const double A1p76[3][3] = {
- { 82758667.0, 1871391091.0, 4127413238.0 },
- { 3672831523.0, 69195019.0, 1871391091.0 },
- { 3672091415.0, 3528743235.0, 69195019.0 }
-};
-
-const double A2p76[3][3] = {
- { 1511326704.0, 3759209742.0, 1610795712.0 },
- { 4292754251.0, 1511326704.0, 3889917532.0 },
- { 3859662829.0, 4292754251.0, 3708466080.0 }
-};
-
-const double A1p127[3][3] = {
- { 2427906178.0, 3580155704.0, 949770784.0 },
- { 226153695.0, 1230515664.0, 3580155704.0 },
- { 1988835001.0, 986791581.0, 1230515664.0 }
-};
-
-const double A2p127[3][3] = {
- { 1464411153.0, 277697599.0, 1610723613.0 },
- { 32183930.0, 1464411153.0, 1022607788.0 },
- { 2824425944.0, 32183930.0, 2093834863.0 }
-};
-
-
//-------------------------------------------------------------------------
// Return (a*s + c) MOD m; a, s, c and m must be < 2^35
@@ -101,16 +75,26 @@
v = a * s + c;
- if (v >= two53 || v <= -two53) {
- a1 = static_cast<int32_t> (a / two17); a -= a1 * two17;
+ if (v >= two53 || v <= -two53)
+ {
+ a1 = static_cast<int32_t> (a / two17);
+ a -= a1 * two17;
v = a1 * s;
- a1 = static_cast<int32_t> (v / m); v -= a1 * m;
+ a1 = static_cast<int32_t> (v / m);
+ v -= a1 * m;
v = v * two17 + a * s + c;
}
a1 = static_cast<int32_t> (v / m);
/* in case v < 0)*/
- if ((v -= a1 * m) < 0.0) return v += m; else return v;
+ if ((v -= a1 * m) < 0.0)
+ {
+ return v += m;
+ }
+ else
+ {
+ return v;
+ }
}
@@ -118,19 +102,22 @@
// Compute the vector v = A*s MOD m. Assume that -m < s[i] < m.
// Works also when v = s.
//
-void MatVecModM (const double A[3][3], const double s[3], double v[3],
+void MatVecModM (const Matrix A, const double s[3], double v[3],
double m)
{
int i;
double x[3]; // Necessary if v = s
- for (i = 0; i < 3; ++i) {
+ for (i = 0; i < 3; ++i)
+ {
x[i] = MultModM (A[i][0], s[0], 0.0, m);
x[i] = MultModM (A[i][1], s[1], x[i], m);
x[i] = MultModM (A[i][2], s[2], x[i], m);
}
for (i = 0; i < 3; ++i)
- v[i] = x[i];
+ {
+ v[i] = x[i];
+ }
}
@@ -138,41 +125,55 @@
// Compute the matrix C = A*B MOD m. Assume that -m < s[i] < m.
// Note: works also if A = C or B = C or A = B = C.
//
-void MatMatModM (const double A[3][3], const double B[3][3],
- double C[3][3], double m)
+void MatMatModM (const Matrix A, const Matrix B,
+ Matrix C, double m)
{
int i, j;
- double V[3], W[3][3];
+ double V[3];
+ Matrix W;
- for (i = 0; i < 3; ++i) {
+ for (i = 0; i < 3; ++i)
+ {
for (j = 0; j < 3; ++j)
- V[j] = B[j][i];
+ {
+ V[j] = B[j][i];
+ }
MatVecModM (A, V, V, m);
for (j = 0; j < 3; ++j)
- W[j][i] = V[j];
+ {
+ W[j][i] = V[j];
+ }
}
for (i = 0; i < 3; ++i)
- for (j = 0; j < 3; ++j)
- C[i][j] = W[i][j];
+ {
+ for (j = 0; j < 3; ++j)
+ {
+ C[i][j] = W[i][j];
+ }
+ }
}
//-------------------------------------------------------------------------
// Compute the matrix B = (A^(2^e) Mod m); works also if A = B.
//
-void MatTwoPowModM (const double A[3][3], double B[3][3], double m, int32_t e)
+void MatTwoPowModM (const Matrix src, Matrix dst, double m, int32_t e)
{
int i, j;
- /* initialize: B = A */
- if (A != B) {
- for (i = 0; i < 3; ++i)
- for (j = 0; j < 3; ++j)
- B[i][j] = A[i][j];
+ /* initialize: dst = src */
+ for (i = 0; i < 3; ++i)
+ {
+ for (j = 0; j < 3; ++j)
+ {
+ dst[i][j] = src[i][j];
+ }
}
- /* Compute B = A^(2^e) mod m */
+ /* Compute dst = src^(2^e) mod m */
for (i = 0; i < e; i++)
- MatMatModM (B, B, B, m);
+ {
+ MatMatModM (dst, dst, dst, m);
+ }
}
@@ -186,30 +187,60 @@
/* initialize: W = A; B = I */
for (i = 0; i < 3; ++i)
- for (j = 0; j < 3; ++j) {
- W[i][j] = A[i][j];
- B[i][j] = 0.0;
- }
+ {
+ for (j = 0; j < 3; ++j)
+ {
+ W[i][j] = A[i][j];
+ B[i][j] = 0.0;
+ }
+ }
for (j = 0; j < 3; ++j)
- B[j][j] = 1.0;
+ {
+ B[j][j] = 1.0;
+ }
/* Compute B = A^n mod m using the binary decomposition of n */
- while (n > 0) {
- if (n % 2) MatMatModM (W, B, B, m);
+ while (n > 0)
+ {
+ if (n % 2)
+ {
+ MatMatModM (W, B, B, m);
+ }
MatMatModM (W, W, W, m);
n /= 2;
}
}
-
-static ns3::GlobalValue g_rngSeed ("RngSeed",
- "The global seed of all rng streams",
- ns3::IntegerValue (1),
- ns3::MakeIntegerChecker<uint32_t> ());
-static ns3::GlobalValue g_rngRun ("RngRun",
- "The run number used to modify the global seed",
- ns3::IntegerValue (1),
- ns3::MakeIntegerChecker<uint32_t> ());
+// The following are the transition matrices of the two MRG components
+// (in matrix form), raised to all powers of 2 from 1 to 191
+struct Precalculated
+{
+ Matrix a1[190];
+ Matrix a2[190];
+};
+struct Precalculated PowerOfTwoConstants (void)
+{
+ struct Precalculated precalculated;
+ for (int i = 0; i < 190; i++)
+ {
+ int power = i + 1;
+ MatTwoPowModM (A1p0, precalculated.a1[i], m1, power);
+ MatTwoPowModM (A2p0, precalculated.a2[i], m2, power);
+ }
+ return precalculated;
+}
+void PowerOfTwoMatrix (int n, Matrix a1p, Matrix a2p)
+{
+ static struct Precalculated constants = PowerOfTwoConstants ();
+ for (int i = 0; i < 3; i ++)
+ {
+ for (int j = 0; j < 3; j++)
+ {
+ a1p[i][j] = constants.a1[n-1][i][j];
+ a2p[i][j] = constants.a2[n-1][i][j];
+ }
+ }
+}
} // end of anonymous namespace
@@ -218,348 +249,74 @@
//-------------------------------------------------------------------------
// Generate the next random number.
//
-double RngStream::U01 ()
+double RngStream::RandU01 ()
{
int32_t k;
double p1, p2, u;
/* Component 1 */
- p1 = a12 * Cg[1] - a13n * Cg[0];
+ p1 = a12 * m_currentState[1] - a13n * m_currentState[0];
k = static_cast<int32_t> (p1 / m1);
p1 -= k * m1;
- if (p1 < 0.0) p1 += m1;
- Cg[0] = Cg[1]; Cg[1] = Cg[2]; Cg[2] = p1;
+ if (p1 < 0.0)
+ {
+ p1 += m1;
+ }
+ m_currentState[0] = m_currentState[1]; m_currentState[1] = m_currentState[2]; m_currentState[2] = p1;
/* Component 2 */
- p2 = a21 * Cg[5] - a23n * Cg[3];
+ p2 = a21 * m_currentState[5] - a23n * m_currentState[3];
k = static_cast<int32_t> (p2 / m2);
p2 -= k * m2;
- if (p2 < 0.0) p2 += m2;
- Cg[3] = Cg[4]; Cg[4] = Cg[5]; Cg[5] = p2;
+ if (p2 < 0.0)
+ {
+ p2 += m2;
+ }
+ m_currentState[3] = m_currentState[4]; m_currentState[4] = m_currentState[5]; m_currentState[5] = p2;
/* Combination */
u = ((p1 > p2) ? (p1 - p2) * norm : (p1 - p2 + m1) * norm);
- return (anti == false) ? u : (1 - u);
-}
-
-
-//-------------------------------------------------------------------------
-// Generate the next random number with extended (53 bits) precision.
-//
-double RngStream::U01d ()
-{
- double u;
- u = U01 ();
- if (anti) {
- // Don't forget that U01() returns 1 - u in the antithetic case
- u += (U01 () - 1.0) * fact;
- return (u < 0.0) ? u + 1.0 : u;
- } else {
- u += U01 () * fact;
- return (u < 1.0) ? u : (u - 1.0);
- }
+ return u;
}
-//-------------------------------------------------------------------------
-// Check that the seeds are legitimate values. Returns true if legal seeds,
-// false otherwise.
-//
-bool RngStream::CheckSeed (const uint32_t seed[6])
-{
- int i;
-
- for (i = 0; i < 3; ++i) {
- if (seed[i] >= m1) {
- cerr << "****************************************\n\n"
- << "ERROR: Seed[" << i << "] >= 4294967087, Seed is not set."
- << "\n\n****************************************\n\n";
- return (false);
- }
- }
- for (i = 3; i < 6; ++i) {
- if (seed[i] >= m2) {
- cerr << "Seed[" << i << "] = " << seed[i] << endl;
- cerr << "*****************************************\n\n"
- << "ERROR: Seed[" << i << "] >= 4294944443, Seed is not set."
- << "\n\n*****************************************\n\n";
- return (false);
- }
- }
- if (seed[0] == 0 && seed[1] == 0 && seed[2] == 0) {
- cerr << "****************************\n\n"
- << "ERROR: First 3 seeds = 0.\n\n"
- << "****************************\n\n";
- return (false);
- }
- if (seed[3] == 0 && seed[4] == 0 && seed[5] == 0) {
- cerr << "****************************\n\n"
- << "ERROR: Last 3 seeds = 0.\n\n"
- << "****************************\n\n";
- return (false);
- }
- return true;
-}
-
-uint32_t
-RngStream::EnsureGlobalInitialized (void)
+RngStream::RngStream (uint32_t seedNumber, uint64_t stream, uint64_t substream)
{
- static bool initialized = false;
- static uint32_t run = 0;
- if (!initialized)
+ if (seedNumber >= m1 || seedNumber >= m2 || seedNumber == 0)
{
- initialized = true;
- uint32_t seed;
- IntegerValue value;
- g_rngSeed.GetValue (value);
- seed = value.Get ();
- g_rngRun.GetValue (value);
- run = value.Get ();
- SetPackageSeed (seed);
+ NS_FATAL_ERROR ("invalid Seed " << seedNumber);
}
- return run;
-}
-
-//*************************************************************************
-// Public members of the class start here
-
-
-//-------------------------------------------------------------------------
-// The default seed of the package; at run time, this will be overwritten
-// by the g_rngSeed value the first time RngStream::RngStream is visited;
-// so the value 12345 is for debugging
-//
-double RngStream::nextSeed[6] =
-{
- 12345.0, 12345.0, 12345.0, 12345.0, 12345.0, 12345.0
-};
-
-//-------------------------------------------------------------------------
-// constructor
-//
-RngStream::RngStream ()
-{
- uint32_t run = EnsureGlobalInitialized ();
-
- anti = false;
- incPrec = false;
- // Stream initialization moved to separate method.
- InitializeStream ();
- //move the state of this stream up
- ResetNthSubstream (run);
+ for (int i = 0; i < 6; ++i)
+ {
+ m_currentState[i] = seedNumber;
+ }
+ AdvanceNthBy (stream, 127, m_currentState);
+ AdvanceNthBy (substream, 76, m_currentState);
}
RngStream::RngStream(const RngStream& r)
{
- anti = r.anti;
- incPrec = r.incPrec;
- for (int i = 0; i < 6; ++i) {
- Cg[i] = r.Cg[i];
- Bg[i] = r.Bg[i];
- Ig[i] = r.Ig[i];
+ for (int i = 0; i < 6; ++i)
+ {
+ m_currentState[i] = r.m_currentState[i];
}
}
-
-void RngStream::InitializeStream ()
-{ // Moved from the RngStream constructor above to allow seeding
- // AFTER the global package seed has been set in the Random
- // object constructor.
- /* Information on a stream. The arrays {Cg, Bg, Ig} contain the current
- state of the stream, the starting state of the current SubStream, and the
- starting state of the stream. This stream generates antithetic variates
- if anti = true. It also generates numbers with extended precision (53
- bits if machine follows IEEE 754 standard) if incPrec = true. nextSeed
- will be the seed of the next declared RngStream. */
-
- for (int i = 0; i < 6; ++i) {
- Bg[i] = Cg[i] = Ig[i] = nextSeed[i];
- }
-
- MatVecModM (A1p127, nextSeed, nextSeed, m1);
- MatVecModM (A2p127, &nextSeed[3], &nextSeed[3], m2);
-}
-
-//-------------------------------------------------------------------------
-// Reset Stream to beginning of Stream.
-//
-void RngStream::ResetStartStream ()
-{
- for (int i = 0; i < 6; ++i)
- Cg[i] = Bg[i] = Ig[i];
-}
-
-
-//-------------------------------------------------------------------------
-// Reset Stream to beginning of SubStream.
-//
-void RngStream::ResetStartSubstream ()
-{
- for (int i = 0; i < 6; ++i)
- Cg[i] = Bg[i];
-}
-
-
-//-------------------------------------------------------------------------
-// Reset Stream to NextSubStream.
-//
-void RngStream::ResetNextSubstream ()
+void
+RngStream::AdvanceNthBy (uint64_t nth, int by, double state[6])
{
- MatVecModM (A1p76, Bg, Bg, m1);
- MatVecModM (A2p76, &Bg[3], &Bg[3], m2);
- for (int i = 0; i < 6; ++i)
- Cg[i] = Bg[i];
-}
-
-//-------------------------------------------------------------------------
-// Reset Stream to Nth SubStream.
-//
-void RngStream::ResetNthSubstream (uint32_t N)
-{
- if(N==0) return;
- for(uint32_t i=0; i<N; ++i) {
- MatVecModM (A1p76, Bg, Bg, m1);
- MatVecModM (A2p76, &Bg[3], &Bg[3], m2);
- }
- for (int i = 0; i < 6; ++i)
- Cg[i] = Bg[i];
-}
-
-//-------------------------------------------------------------------------
-bool RngStream::SetPackageSeed (const uint32_t seed[6])
-{
- if (!CheckSeed (seed))
+ Matrix matrix1,matrix2;
+ for (int i = 0; i < 64; i++)
{
- return false;
- }
- for (int i = 0; i < 6; ++i)
- nextSeed[i] = seed[i];
- return true;
-}
-bool
-RngStream::SetPackageSeed (uint32_t seed)
-{
- uint32_t seeds[6] = { seed, seed, seed, seed, seed, seed};
- return SetPackageSeed (seeds);
-}
-void
-RngStream::GetPackageSeed (uint32_t seed[6])
-{
- IntegerValue value;
- g_rngSeed.GetValue (value);
- uint32_t theSeed = value.Get ();
- for (int i = 0; i < 6; i++)
- {
- seed[i] = static_cast<uint32_t> (theSeed);
+ int nbit = 63 - i;
+ int bit = (nth >> nbit) & 0x1;
+ if (bit)
+ {
+ PowerOfTwoMatrix(by + nbit, matrix1, matrix2);
+ MatVecModM (matrix1, state, state, m1);
+ MatVecModM (matrix2, &state[3], &state[3], m2);
+ }
}
}
-void
-RngStream::SetPackageRun (uint32_t run)
-{
- g_rngRun.SetValue (IntegerValue (run));
-}
-uint32_t
-RngStream::GetPackageRun (void)
-{
- IntegerValue run;
- g_rngRun.GetValue (run);
- return run.Get ();
-}
-bool
-RngStream::CheckSeed (uint32_t seed)
-{
- uint32_t seeds[6] = { seed, seed, seed, seed, seed, seed};
- return CheckSeed (seeds);
-}
-
-
-
-//-------------------------------------------------------------------------
-bool RngStream::SetSeeds (const uint32_t seed[6])
-{
- if (!CheckSeed (seed)) return false;
- for (int i = 0; i < 6; ++i)
- Cg[i] = Bg[i] = Ig[i] = seed[i];
- return true;
-}
-
-
-//-------------------------------------------------------------------------
-// if e > 0, let n = 2^e + c;
-// if e < 0, let n = -2^(-e) + c;
-// if e = 0, let n = c.
-// Jump n steps forward if n > 0, backwards if n < 0.
-//
-void RngStream::AdvanceState (int32_t e, int32_t c)
-{
- double B1[3][3], C1[3][3], B2[3][3], C2[3][3];
-
- if (e > 0) {
- MatTwoPowModM (A1p0, B1, m1, e);
- MatTwoPowModM (A2p0, B2, m2, e);
- } else if (e < 0) {
- MatTwoPowModM (InvA1, B1, m1, -e);
- MatTwoPowModM (InvA2, B2, m2, -e);
- }
-
- if (c >= 0) {
- MatPowModM (A1p0, C1, m1, c);
- MatPowModM (A2p0, C2, m2, c);
- } else {
- MatPowModM (InvA1, C1, m1, -c);
- MatPowModM (InvA2, C2, m2, -c);
- }
-
- if (e) {
- MatMatModM (B1, C1, C1, m1);
- MatMatModM (B2, C2, C2, m2);
- }
-
- MatVecModM (C1, Cg, Cg, m1);
- MatVecModM (C2, &Cg[3], &Cg[3], m2);
-}
-
-
-//-------------------------------------------------------------------------
-void RngStream::GetState (uint32_t seed[6]) const
-{
- for (int i = 0; i < 6; ++i)
- seed[i] = static_cast<uint32_t> (Cg[i]);
-}
-
-
-//-------------------------------------------------------------------------
-void RngStream::IncreasedPrecis (bool incp)
-{
- incPrec = incp;
-}
-
-
-//-------------------------------------------------------------------------
-void RngStream::SetAntithetic (bool a)
-{
- anti = a;
-}
-
-
-//-------------------------------------------------------------------------
-// Generate the next random number.
-//
-double RngStream::RandU01 ()
-{
- if (incPrec)
- return U01d ();
- else
- return U01 ();
-}
-
-
-//-------------------------------------------------------------------------
-// Generate the next random integer.
-//
-int32_t RngStream::RandInt (int32_t low, int32_t high)
-{
- return low + static_cast<int32_t> ((high - low + 1) * RandU01 ());
-};
} // namespace ns3
--- a/src/core/model/rng-stream.h Mon Jul 09 06:53:14 2012 -0700
+++ b/src/core/model/rng-stream.h Tue Jul 10 21:31:47 2012 -0700
@@ -36,38 +36,21 @@
* class are explained in:
* http://www.iro.umontreal.ca/~lecuyer/myftp/papers/streams00.pdf
*/
-class RngStream {
-public: //public api
- RngStream ();
+class RngStream
+{
+public:
+ RngStream (uint32_t seed, uint64_t stream, uint64_t substream);
RngStream (const RngStream&);
- void InitializeStream (); // Separate initialization
- void ResetStartStream ();
- void ResetStartSubstream ();
- void ResetNextSubstream ();
- void ResetNthSubstream (uint32_t N);
- void SetAntithetic (bool a);
- void IncreasedPrecis (bool incp);
- bool SetSeeds (const uint32_t seed[6]);
- void AdvanceState (int32_t e, int32_t c);
- void GetState (uint32_t seed[6]) const;
- double RandU01 ();
- int32_t RandInt (int32_t i, int32_t j);
-public: //public static api
- static bool SetPackageSeed (uint32_t seed);
- static bool SetPackageSeed (const uint32_t seed[6]);
- static void GetPackageSeed (uint32_t seed[6]);
- static void SetPackageRun (uint32_t run);
- static uint32_t GetPackageRun (void);
- static bool CheckSeed (const uint32_t seed[6]);
- static bool CheckSeed (uint32_t seed);
-private: //members
- double Cg[6], Bg[6], Ig[6];
- bool anti, incPrec;
- double U01 ();
- double U01d ();
- static uint32_t EnsureGlobalInitialized (void);
-private: //static data
- static double nextSeed[6];
+ /**
+ * Generate the next random number for this stream.
+ * Uniformly distributed between 0 and 1.
+ */
+ double RandU01 (void);
+
+private:
+ void AdvanceNthBy (uint64_t nth, int by, double state[6]);
+
+ double m_currentState[6];
};
} // namespace ns3
--- a/src/core/wscript Mon Jul 09 06:53:14 2012 -0700
+++ b/src/core/wscript Tue Jul 10 21:31:47 2012 -0700
@@ -124,6 +124,9 @@
'model/object.cc',
'model/test.cc',
'model/random-variable.cc',
+ 'model/random-variable-stream.cc',
+ 'model/random-variable-stream-helper.cc',
+ 'model/rng-seed-manager.cc',
'model/rng-stream.cc',
'model/command-line.cc',
'model/type-name.cc',
@@ -160,6 +163,7 @@
'test/object-test-suite.cc',
'test/ptr-test-suite.cc',
'test/random-variable-test-suite.cc',
+ 'test/random-variable-stream-test-suite.cc',
'test/sample-test-suite.cc',
'test/simulator-test-suite.cc',
'test/time-test-suite.cc',
@@ -206,6 +210,9 @@
'model/fatal-error.h',
'model/test.h',
'model/random-variable.h',
+ 'model/random-variable-stream.h',
+ 'model/random-variable-stream-helper.h',
+ 'model/rng-seed-manager.h',
'model/rng-stream.h',
'model/command-line.h',
'model/type-name.h',