rng-stream.cc

Go to the documentation of this file.

//
//  Copyright (C) 2001  Pierre L'Ecuyer (lecuyer@iro.umontreal.ca)
//
// 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
//
// Modified for ns-3 by:
//   - Rajib Bhattacharjea<raj.b@gatech.edu>
//   - Mathieu Lacage <mathieu.lacage@gmail.com>
//
 
#include "rng-stream.h"
 
#include "fatal-error.h"
#include "log.h"
 
#include <cstdlib>
#include <iostream>
 
/**
 * \file
 * \ingroup rngimpl
 * ns3::RngStream and MRG32k3a implementations.
 */
 
namespace ns3
{
 
// Note:  Logging in this file is largely avoided due to the
// number of calls that are made to these functions and the possibility
// of causing recursions leading to stack overflow
NS_LOG_COMPONENT_DEFINE("RngStream");
 
} // namespace ns3
 
/**
 * \ingroup rngimpl
 * @{
 */
/** Namespace for MRG32k3a implementation details. */
namespace MRG32k3a
{
 
// clang-format off
 
/** Type for 3x3 matrix of doubles. */
typedef double Matrix[3][3];
 
/** First component modulus, 2<sup>32</sup> - 209. */
const double m1   =       4294967087.0;
 
/** Second component modulus, 2<sup>32</sup> - 22853. */
const double m2   =       4294944443.0;
 
/** Normalization to obtain randoms on [0,1). */
const double norm =       1.0 / (m1 + 1.0);
 
/** First component multiplier of <i>n</i> - 2 value. */
const double a12  =       1403580.0;
 
/** First component multiplier of <i>n</i> - 3 value. */
const double a13n =       810728.0;
 
/** Second component multiplier of <i>n</i> - 1 value. */
const double a21  =       527612.0;
 
/** Second component multiplier of <i>n</i> - 3 value. */
const double a23n =       1370589.0;
 
/** Decomposition factor for computing a*s in less than 53 bits, 2<sup>17</sup> */
const double two17 =      131072.0;
 
/** IEEE-754 floating point precision, 2<sup>53</sup> */
const double two53 =      9007199254740992.0;
 
/** First component transition matrix. */
const Matrix A1p0 = {
  {       0.0,        1.0,       0.0 },
  {       0.0,        0.0,       1.0 },
  { -810728.0,  1403580.0,       0.0 }
};
 
/** Second component transition matrix. */
const Matrix A2p0 = {
  {        0.0,        1.0,       0.0 },
  {        0.0,        0.0,       1.0 },
  { -1370589.0,        0.0,  527612.0 }
};
 
 
//-------------------------------------------------------------------------
/**
 * Return (a*s + c) MOD m; a, s, c and m must be < 2^35
 *
 * This computes the result exactly, without exceeding the 53 bit
 * precision of doubles.
 *
 * \param [in] a First multiplicative argument.
 * \param [in] s Second multiplicative argument.
 * \param [in] c Additive argument.
 * \param [in] m Modulus.
 * \returns <tt>(a*s +c) MOD m</tt>
 */
double MultModM (double a, double s, double c, double m)
{
  double v;
  int32_t a1;
 
  v = a * s + c;
 
  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;
      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;
    }
}
 
 
//-------------------------------------------------------------------------
/**
 * Compute the vector v = A*s MOD m. Assume that -m < s[i] < m.
 * Works also when v = s.
 *
 * \param [in] A Matrix argument, 3x3.
 * \param [in] s Three component input vector.
 * \param [out] v Three component output vector.
 * \param [in] m Modulus.
 */
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)
    {
      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];
    }
}
 
 
//-------------------------------------------------------------------------
/**
 * 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.
 *
 * \param [in] A First matrix argument.
 * \param [in] B Second matrix argument.
 * \param [out] C Result matrix.
 * \param [in] m Modulus.
 */
void MatMatModM (const Matrix A, const Matrix B,
                 Matrix C, double m)
{
  int i;
  int j;
  double V[3];
  Matrix W;
 
  for (i = 0; i < 3; ++i)
    {
      for (j = 0; j < 3; ++j)
        {
          V[j] = B[j][i];
        }
      MatVecModM (A, V, V, m);
      for (j = 0; j < 3; ++j)
        {
          W[j][i] = V[j];
        }
    }
  for (i = 0; i < 3; ++i)
    {
      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.
 *
 * \param [in] src Matrix input argument \c A.
 * \param [out] dst Matrix output \c B.
 * \param [in] m Modulus.
 * \param [in] e The exponent.
 */
void MatTwoPowModM (const Matrix src, Matrix dst, double m, int32_t e)
{
  int i;
  int j;
 
  /* initialize: dst = src */
  for (i = 0; i < 3; ++i)
    {
      for (j = 0; j < 3; ++j)
        {
          dst[i][j] = src[i][j];
        }
    }
  /* Compute dst = src^(2^e) mod m */
  for (i = 0; i < e; i++)
    {
      MatMatModM (dst, dst, dst, m);
    }
}
 
 
//-------------------------------------------------------------------------
/**
 * Compute the matrix B = (A^n Mod m);  works even if A = B.
 *
 * \param [in] A Matrix input argument.
 * \param [out] B Matrix output.
 * \param [in] m Modulus.
 * \param [in] n Exponent.
 */
void MatPowModM (const double A[3][3], double B[3][3], double m, int32_t n)
{
  int i;
  int j;
  double W[3][3];
 
  // 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)
    {
      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);
        }
      MatMatModM (W, W, W, m);
      n /= 2;
    }
}
 
/**
 * 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];  //!< First component transition matrix powers.
  Matrix a2[190];  //!< Second component transition matrix powers.
};
 
/**
 * Compute the transition matrices of the two MRG components
 * raised to all powers of 2 from 1 to 191.
 *
 * \returns The precalculated powers of the transition matrices.
 */
Precalculated PowerOfTwoConstants ()
{
  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;
}
/**
 * Get the transition matrices raised to a power of 2.
 *
 * \param [in] n The power of 2.
 * \param [out] a1p The first transition matrix power.
 * \param [out] a2p The second transition matrix power.
 */
void PowerOfTwoMatrix (int n, Matrix a1p, Matrix a2p)
{
  static  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];
        }
    }
}
 
} // namespace MRG32k3a
 
// clang-format on
 
namespace ns3
{
 
using namespace MRG32k3a;
 
double
RngStream::RandU01()
{
    int32_t k;
    double p1;
    double p2;
    double u;
 
    /* Component 1 */
    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;
    }
    m_currentState[0] = m_currentState[1];
    m_currentState[1] = m_currentState[2];
    m_currentState[2] = p1;
 
    /* Component 2 */
    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;
    }
    m_currentState[3] = m_currentState[4];
    m_currentState[4] = m_currentState[5];
    m_currentState[5] = p2;
 
    /* Combination */
    u = ((p1 > p2) ? (p1 - p2) * MRG32k3a::norm : (p1 - p2 + m1) * MRG32k3a::norm);
 
    return u;
}
 
RngStream::RngStream(uint32_t seedNumber, uint64_t stream, uint64_t substream)
{
    if (seedNumber >= m1 || seedNumber >= m2 || seedNumber == 0)
    {
        NS_FATAL_ERROR("invalid Seed " << seedNumber);
    }
    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)
{
    for (int i = 0; i < 6; ++i)
    {
        m_currentState[i] = r.m_currentState[i];
    }
}
 
void
RngStream::AdvanceNthBy(uint64_t nth, int by, double state[6])
{
    Matrix matrix1;
    Matrix matrix2;
    for (int i = 0; i < 64; i++)
    {
        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);
        }
    }
}
 
} // namespace ns3
 
/**@}*/ // \ingroup rngimpl