src/core/rng-stream.cc
author Mathieu Lacage <mathieu.lacage@sophia.inria.fr>
Sat, 04 Jul 2009 08:15:48 +0200
changeset 4654 2eaebe77d66b
parent 4235 c070d2fca46d
permissions -rw-r--r--
Added tag ns-3.5 for changeset c975274c9707
tomh@346
     1
/* -*- Mode:C++; c-file-style:"gnu"; indent-tabs-mode:nil; -*- */
tomh@346
     2
//
tomh@346
     3
//  Copyright (C) 2001  Pierre L'Ecuyer (lecuyer@iro.umontreal.ca)
tomh@346
     4
//
tomh@346
     5
// This program is free software; you can redistribute it and/or modify
tomh@346
     6
// it under the terms of the GNU General Public License version 2 as
tomh@346
     7
// published by the Free Software Foundation;
tomh@346
     8
//
tomh@346
     9
// This program is distributed in the hope that it will be useful,
tomh@346
    10
// but WITHOUT ANY WARRANTY; without even the implied warranty of
tomh@346
    11
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
tomh@346
    12
// GNU General Public License for more details.
tomh@346
    13
//
tomh@346
    14
// You should have received a copy of the GNU General Public License
tomh@346
    15
// along with this program; if not, write to the Free Software
tomh@346
    16
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
tomh@346
    17
//
tomh@346
    18
// Modified for ns-3 by: Rajib Bhattacharjea<raj.b@gatech.edu>
tomh@346
    19
//
tomh@346
    20
tomh@346
    21
#include <cstdlib>
tomh@346
    22
#include <iostream>
tomh@346
    23
#include "rng-stream.h"
raj@4221
    24
#include "global-value.h"
raj@4221
    25
#include "integer.h"
tomh@346
    26
using namespace std;
tomh@346
    27
tomh@346
    28
namespace
tomh@346
    29
{
tomh@346
    30
const double m1   =       4294967087.0;
tomh@346
    31
const double m2   =       4294944443.0;
tomh@346
    32
const double norm =       1.0 / (m1 + 1.0);
tomh@346
    33
const double a12  =       1403580.0;
tomh@346
    34
const double a13n =       810728.0;
tomh@346
    35
const double a21  =       527612.0;
tomh@346
    36
const double a23n =       1370589.0;
tomh@346
    37
const double two17 =      131072.0;
tomh@346
    38
const double two53 =      9007199254740992.0;
tomh@346
    39
const double fact =       5.9604644775390625e-8;     /* 1 / 2^24  */
tomh@346
    40
tomh@346
    41
// The following are the transition matrices of the two MRG components
tomh@346
    42
// (in matrix form), raised to the powers -1, 1, 2^76, and 2^127, resp.
tomh@346
    43
tomh@346
    44
const double InvA1[3][3] = {          // Inverse of A1p0
tomh@346
    45
       { 184888585.0,   0.0,  1945170933.0 },
tomh@346
    46
       {         1.0,   0.0,           0.0 },
tomh@346
    47
       {         0.0,   1.0,           0.0 }
tomh@346
    48
       };
tomh@346
    49
tomh@346
    50
const double InvA2[3][3] = {          // Inverse of A2p0
tomh@346
    51
       {      0.0,  360363334.0,  4225571728.0 },
tomh@346
    52
       {      1.0,          0.0,           0.0 },
tomh@346
    53
       {      0.0,          1.0,           0.0 }
tomh@346
    54
       };
tomh@346
    55
tomh@346
    56
const double A1p0[3][3] = {
tomh@346
    57
       {       0.0,        1.0,       0.0 },
tomh@346
    58
       {       0.0,        0.0,       1.0 },
tomh@346
    59
       { -810728.0,  1403580.0,       0.0 }
tomh@346
    60
       };
tomh@346
    61
tomh@346
    62
const double A2p0[3][3] = {
tomh@346
    63
       {        0.0,        1.0,       0.0 },
tomh@346
    64
       {        0.0,        0.0,       1.0 },
tomh@346
    65
       { -1370589.0,        0.0,  527612.0 }
tomh@346
    66
       };
tomh@346
    67
tomh@346
    68
const double A1p76[3][3] = {
tomh@346
    69
       {      82758667.0, 1871391091.0, 4127413238.0 },
tomh@346
    70
       {    3672831523.0,   69195019.0, 1871391091.0 },
tomh@346
    71
       {    3672091415.0, 3528743235.0,   69195019.0 }
tomh@346
    72
       };
tomh@346
    73
tomh@346
    74
const double A2p76[3][3] = {
tomh@346
    75
       {    1511326704.0, 3759209742.0, 1610795712.0 },
tomh@346
    76
       {    4292754251.0, 1511326704.0, 3889917532.0 },
tomh@346
    77
       {    3859662829.0, 4292754251.0, 3708466080.0 }
tomh@346
    78
       };
tomh@346
    79
tomh@346
    80
const double A1p127[3][3] = {
tomh@346
    81
       {    2427906178.0, 3580155704.0,  949770784.0 },
tomh@346
    82
       {     226153695.0, 1230515664.0, 3580155704.0 },
tomh@346
    83
       {    1988835001.0,  986791581.0, 1230515664.0 }
tomh@346
    84
       };
tomh@346
    85
tomh@346
    86
const double A2p127[3][3] = {
tomh@346
    87
       {    1464411153.0,  277697599.0, 1610723613.0 },
tomh@346
    88
       {      32183930.0, 1464411153.0, 1022607788.0 },
tomh@346
    89
       {    2824425944.0,   32183930.0, 2093834863.0 }
tomh@346
    90
       };
tomh@346
    91
tomh@346
    92
tomh@346
    93
tomh@346
    94
//-------------------------------------------------------------------------
tomh@346
    95
// Return (a*s + c) MOD m; a, s, c and m must be < 2^35
tomh@346
    96
//
tomh@346
    97
double MultModM (double a, double s, double c, double m)
tomh@346
    98
{
tomh@346
    99
    double v;
tomh@346
   100
    int32_t a1;
tomh@346
   101
tomh@346
   102
    v = a * s + c;
tomh@346
   103
tomh@346
   104
    if (v >= two53 || v <= -two53) {
tomh@346
   105
        a1 = static_cast<int32_t> (a / two17);    a -= a1 * two17;
tomh@346
   106
        v  = a1 * s;
tomh@346
   107
        a1 = static_cast<int32_t> (v / m);     v -= a1 * m;
tomh@346
   108
        v = v * two17 + a * s + c;
tomh@346
   109
    }
tomh@346
   110
tomh@346
   111
    a1 = static_cast<int32_t> (v / m);
tomh@346
   112
    /* in case v < 0)*/
tomh@346
   113
    if ((v -= a1 * m) < 0.0) return v += m;   else return v;
tomh@346
   114
}
tomh@346
   115
tomh@346
   116
tomh@346
   117
//-------------------------------------------------------------------------
tomh@346
   118
// Compute the vector v = A*s MOD m. Assume that -m < s[i] < m.
tomh@346
   119
// Works also when v = s.
tomh@346
   120
//
tomh@346
   121
void MatVecModM (const double A[3][3], const double s[3], double v[3],
tomh@346
   122
                 double m)
tomh@346
   123
{
tomh@346
   124
    int i;
tomh@346
   125
    double x[3];               // Necessary if v = s
tomh@346
   126
tomh@346
   127
    for (i = 0; i < 3; ++i) {
tomh@346
   128
        x[i] = MultModM (A[i][0], s[0], 0.0, m);
tomh@346
   129
        x[i] = MultModM (A[i][1], s[1], x[i], m);
tomh@346
   130
        x[i] = MultModM (A[i][2], s[2], x[i], m);
tomh@346
   131
    }
tomh@346
   132
    for (i = 0; i < 3; ++i)
tomh@346
   133
        v[i] = x[i];
tomh@346
   134
}
tomh@346
   135
tomh@346
   136
tomh@346
   137
//-------------------------------------------------------------------------
tomh@346
   138
// Compute the matrix C = A*B MOD m. Assume that -m < s[i] < m.
tomh@346
   139
// Note: works also if A = C or B = C or A = B = C.
tomh@346
   140
//
tomh@346
   141
void MatMatModM (const double A[3][3], const double B[3][3],
tomh@346
   142
                 double C[3][3], double m)
tomh@346
   143
{
tomh@346
   144
    int i, j;
tomh@346
   145
    double V[3], W[3][3];
tomh@346
   146
tomh@346
   147
    for (i = 0; i < 3; ++i) {
tomh@346
   148
        for (j = 0; j < 3; ++j)
tomh@346
   149
            V[j] = B[j][i];
tomh@346
   150
        MatVecModM (A, V, V, m);
tomh@346
   151
        for (j = 0; j < 3; ++j)
tomh@346
   152
            W[j][i] = V[j];
tomh@346
   153
    }
tomh@346
   154
    for (i = 0; i < 3; ++i)
tomh@346
   155
        for (j = 0; j < 3; ++j)
tomh@346
   156
            C[i][j] = W[i][j];
tomh@346
   157
}
tomh@346
   158
tomh@346
   159
tomh@346
   160
//-------------------------------------------------------------------------
tomh@346
   161
// Compute the matrix B = (A^(2^e) Mod m);  works also if A = B. 
tomh@346
   162
//
tomh@346
   163
void MatTwoPowModM (const double A[3][3], double B[3][3], double m, int32_t e)
tomh@346
   164
{
tomh@346
   165
   int i, j;
tomh@346
   166
tomh@346
   167
   /* initialize: B = A */
tomh@346
   168
   if (A != B) {
tomh@346
   169
      for (i = 0; i < 3; ++i)
tomh@346
   170
         for (j = 0; j < 3; ++j)
tomh@346
   171
            B[i][j] = A[i][j];
tomh@346
   172
   }
tomh@346
   173
   /* Compute B = A^(2^e) mod m */
tomh@346
   174
   for (i = 0; i < e; i++)
tomh@346
   175
      MatMatModM (B, B, B, m);
tomh@346
   176
}
tomh@346
   177
tomh@346
   178
tomh@346
   179
//-------------------------------------------------------------------------
tomh@346
   180
// Compute the matrix B = (A^n Mod m);  works even if A = B.
tomh@346
   181
//
tomh@346
   182
void MatPowModM (const double A[3][3], double B[3][3], double m, int32_t n)
tomh@346
   183
{
tomh@346
   184
    int i, j;
tomh@346
   185
    double W[3][3];
tomh@346
   186
tomh@346
   187
    /* initialize: W = A; B = I */
tomh@346
   188
    for (i = 0; i < 3; ++i)
tomh@346
   189
        for (j = 0; j < 3; ++j) {
tomh@346
   190
            W[i][j] = A[i][j];
tomh@346
   191
            B[i][j] = 0.0;
tomh@346
   192
        }
tomh@346
   193
    for (j = 0; j < 3; ++j)
tomh@346
   194
        B[j][j] = 1.0;
tomh@346
   195
tomh@346
   196
    /* Compute B = A^n mod m using the binary decomposition of n */
tomh@346
   197
    while (n > 0) {
tomh@346
   198
        if (n % 2) MatMatModM (W, B, B, m);
tomh@346
   199
        MatMatModM (W, W, W, m);
tomh@346
   200
        n /= 2;
tomh@346
   201
    }
tomh@346
   202
}
tomh@346
   203
raj@4234
   204
raj@4235
   205
static ns3::GlobalValue g_rngSeed ("RngSeed", 
raj@4235
   206
                                   "The global seed of all rng streams",
mathieu@4244
   207
                                   ns3::IntegerValue (1),
raj@4235
   208
                                   ns3::MakeIntegerChecker<uint32_t> ());
raj@4235
   209
static ns3::GlobalValue g_rngRun ("RngRun", 
raj@4235
   210
                                  "The run number used to modify the global seed",
mathieu@4244
   211
                                  ns3::IntegerValue (1),
raj@4235
   212
                                  ns3::MakeIntegerChecker<uint32_t> ());
raj@4235
   213
tomh@346
   214
} // end of anonymous namespace
tomh@346
   215
tomh@346
   216
tomh@346
   217
namespace ns3{
tomh@346
   218
//-------------------------------------------------------------------------
tomh@346
   219
// Generate the next random number.
tomh@346
   220
//
tomh@346
   221
double RngStream::U01 ()
tomh@346
   222
{
tomh@346
   223
    int32_t k;
tomh@346
   224
    double p1, p2, u;
tomh@346
   225
tomh@346
   226
    /* Component 1 */
tomh@346
   227
    p1 = a12 * Cg[1] - a13n * Cg[0];
tomh@346
   228
    k = static_cast<int32_t> (p1 / m1);
tomh@346
   229
    p1 -= k * m1;
tomh@346
   230
    if (p1 < 0.0) p1 += m1;
tomh@346
   231
    Cg[0] = Cg[1]; Cg[1] = Cg[2]; Cg[2] = p1;
tomh@346
   232
tomh@346
   233
    /* Component 2 */
tomh@346
   234
    p2 = a21 * Cg[5] - a23n * Cg[3];
tomh@346
   235
    k = static_cast<int32_t> (p2 / m2);
tomh@346
   236
    p2 -= k * m2;
tomh@346
   237
    if (p2 < 0.0) p2 += m2;
tomh@346
   238
    Cg[3] = Cg[4]; Cg[4] = Cg[5]; Cg[5] = p2;
tomh@346
   239
tomh@346
   240
    /* Combination */
tomh@346
   241
    u = ((p1 > p2) ? (p1 - p2) * norm : (p1 - p2 + m1) * norm);
tomh@346
   242
tomh@346
   243
    return (anti == false) ? u : (1 - u);
tomh@346
   244
}
tomh@346
   245
tomh@346
   246
tomh@346
   247
//-------------------------------------------------------------------------
tomh@346
   248
// Generate the next random number with extended (53 bits) precision.
tomh@346
   249
//
tomh@346
   250
double RngStream::U01d ()
tomh@346
   251
{
tomh@346
   252
    double u;
tomh@346
   253
    u = U01();
tomh@346
   254
    if (anti) {
tomh@346
   255
        // Don't forget that U01() returns 1 - u in the antithetic case
tomh@346
   256
        u += (U01() - 1.0) * fact;
tomh@346
   257
        return (u < 0.0) ? u + 1.0 : u;
tomh@346
   258
    } else {
tomh@346
   259
        u += U01() * fact;
tomh@346
   260
        return (u < 1.0) ? u : (u - 1.0);
tomh@346
   261
    }
tomh@346
   262
}
tomh@346
   263
tomh@346
   264
//-------------------------------------------------------------------------
tomh@346
   265
// Check that the seeds are legitimate values. Returns true if legal seeds,
tomh@346
   266
// false otherwise.
tomh@346
   267
//
tomh@346
   268
bool RngStream::CheckSeed (const uint32_t seed[6])
tomh@346
   269
{
tomh@346
   270
    int i;
tomh@346
   271
tomh@346
   272
    for (i = 0; i < 3; ++i) {
tomh@346
   273
        if (seed[i] >= m1) {
tomh@346
   274
            cerr << "****************************************\n\n"
tomh@346
   275
                 << "ERROR: Seed[" << i << "] >= 4294967087, Seed is not set."
tomh@346
   276
                 << "\n\n****************************************\n\n";
tomh@346
   277
            return (false);
tomh@346
   278
        }
tomh@346
   279
    }
tomh@346
   280
    for (i = 3; i < 6; ++i) {
tomh@346
   281
        if (seed[i] >= m2) {
gjc@1020
   282
	  cerr << "Seed[" << i << "] = " << seed[i] << endl; 
gjc@1020
   283
          cerr << "*****************************************\n\n"
gjc@1020
   284
               << "ERROR: Seed[" << i << "] >= 4294944443, Seed is not set."
gjc@1020
   285
               << "\n\n*****************************************\n\n";
gjc@1020
   286
          return (false);
tomh@346
   287
        }
tomh@346
   288
    }
tomh@346
   289
    if (seed[0] == 0 && seed[1] == 0 && seed[2] == 0) {
tomh@346
   290
         cerr << "****************************\n\n"
tomh@346
   291
              << "ERROR: First 3 seeds = 0.\n\n"
tomh@346
   292
              << "****************************\n\n";
tomh@346
   293
         return (false);
tomh@346
   294
    }
tomh@346
   295
    if (seed[3] == 0 && seed[4] == 0 && seed[5] == 0) {
tomh@346
   296
         cerr << "****************************\n\n"
tomh@346
   297
              << "ERROR: Last 3 seeds = 0.\n\n"
tomh@346
   298
              << "****************************\n\n";
tomh@346
   299
         return (false);
tomh@346
   300
    }
tomh@346
   301
    return true;
tomh@346
   302
}
tomh@346
   303
raj@4224
   304
uint32_t
raj@4221
   305
RngStream::EnsureGlobalInitialized (void)
raj@4221
   306
{
raj@4221
   307
  static bool initialized = false;
raj@4224
   308
  static uint32_t run = 0;
raj@4221
   309
  if (!initialized)
raj@4221
   310
    {
raj@4221
   311
      initialized = true;
raj@4221
   312
      uint32_t seed;
mathieu@4244
   313
      IntegerValue value;
mathieu@4244
   314
      g_rngSeed.GetValue (value);
mathieu@4244
   315
      seed = value.Get ();
mathieu@4244
   316
      g_rngRun.GetValue (value);
mathieu@4244
   317
      run = value.Get ();
mathieu@4244
   318
      SetPackageSeed (seed);
raj@4221
   319
    }
raj@4224
   320
  return run;
raj@4221
   321
}
tomh@346
   322
tomh@346
   323
//*************************************************************************
tomh@346
   324
// Public members of the class start here
tomh@346
   325
tomh@346
   326
tomh@346
   327
//-------------------------------------------------------------------------
raj@4234
   328
// The default seed of the package; at run time, this will be overwritten
raj@4234
   329
// by the g_rngSeed value the first time RngStream::RngStream is visited;
raj@4234
   330
// so the value 12345 is for debugging
tomh@346
   331
//
tomh@346
   332
double RngStream::nextSeed[6] =
tomh@346
   333
{
raj@4224
   334
  12345.0, 12345.0, 12345.0, 12345.0, 12345.0, 12345.0
tomh@346
   335
};
tomh@346
   336
tomh@346
   337
//-------------------------------------------------------------------------
tomh@346
   338
// constructor
tomh@346
   339
//
tomh@346
   340
RngStream::RngStream ()
tomh@346
   341
{
mathieu@4244
   342
  uint32_t run = EnsureGlobalInitialized ();
raj@4234
   343
  
raj@4224
   344
  anti = false;
raj@4224
   345
  incPrec = false;
raj@4224
   346
  // Stream initialization moved to separate method.
raj@4224
   347
  InitializeStream ();
raj@4234
   348
  //move the state of this stream up
raj@4224
   349
  ResetNthSubstream (run);
tomh@346
   350
}
tomh@346
   351
raj@360
   352
RngStream::RngStream(const RngStream& r)
raj@360
   353
{
raj@4224
   354
  anti = r.anti;
raj@4224
   355
  incPrec = r.incPrec;
raj@4224
   356
  for (int i = 0; i < 6; ++i) {
raj@4224
   357
    Cg[i] = r.Cg[i];
raj@4224
   358
    Bg[i] = r.Bg[i];
raj@4224
   359
    Ig[i] = r.Ig[i];
raj@4224
   360
  }
raj@360
   361
}
raj@360
   362
      
raj@360
   363
tomh@346
   364
void RngStream::InitializeStream()
tomh@346
   365
{ // Moved from the RngStream constructor above to allow seeding
tomh@346
   366
  // AFTER the global package seed has been set in the Random
tomh@346
   367
  // object constructor.
raj@4224
   368
  /* Information on a stream. The arrays {Cg, Bg, Ig} contain the current
raj@4224
   369
     state of the stream, the starting state of the current SubStream, and the
raj@4224
   370
     starting state of the stream. This stream generates antithetic variates
raj@4224
   371
     if anti = true. It also generates numbers with extended precision (53
raj@4224
   372
     bits if machine follows IEEE 754 standard) if incPrec = true. nextSeed
raj@4224
   373
     will be the seed of the next declared RngStream. */
tomh@346
   374
raj@4224
   375
  for (int i = 0; i < 6; ++i) {
raj@4224
   376
    Bg[i] = Cg[i] = Ig[i] = nextSeed[i];
raj@4224
   377
  }
tomh@346
   378
raj@4224
   379
  MatVecModM (A1p127, nextSeed, nextSeed, m1);
raj@4224
   380
  MatVecModM (A2p127, &nextSeed[3], &nextSeed[3], m2);
tomh@346
   381
}
tomh@346
   382
tomh@346
   383
//-------------------------------------------------------------------------
tomh@346
   384
// Reset Stream to beginning of Stream.
tomh@346
   385
//
tomh@346
   386
void RngStream::ResetStartStream ()
tomh@346
   387
{
raj@4224
   388
  for (int i = 0; i < 6; ++i)
raj@4224
   389
    Cg[i] = Bg[i] = Ig[i];
tomh@346
   390
}
tomh@346
   391
tomh@346
   392
tomh@346
   393
//-------------------------------------------------------------------------
tomh@346
   394
// Reset Stream to beginning of SubStream.
tomh@346
   395
//
tomh@346
   396
void RngStream::ResetStartSubstream ()
tomh@346
   397
{
raj@4224
   398
  for (int i = 0; i < 6; ++i)
raj@4224
   399
    Cg[i] = Bg[i];
tomh@346
   400
}
tomh@346
   401
tomh@346
   402
tomh@346
   403
//-------------------------------------------------------------------------
tomh@346
   404
// Reset Stream to NextSubStream.
tomh@346
   405
//
tomh@346
   406
void RngStream::ResetNextSubstream ()
tomh@346
   407
{
raj@4224
   408
  MatVecModM(A1p76, Bg, Bg, m1);
raj@4224
   409
  MatVecModM(A2p76, &Bg[3], &Bg[3], m2);
raj@4224
   410
  for (int i = 0; i < 6; ++i)
raj@4224
   411
    Cg[i] = Bg[i];
tomh@346
   412
}
tomh@346
   413
raj@364
   414
//-------------------------------------------------------------------------
raj@364
   415
// Reset Stream to Nth SubStream.
raj@364
   416
//
raj@364
   417
void RngStream::ResetNthSubstream (uint32_t N)
raj@364
   418
{
raj@4224
   419
  if(N==0) return;
raj@4224
   420
  for(uint32_t i=0;i<N;++i) {
raj@4224
   421
    MatVecModM(A1p76, Bg, Bg, m1);
raj@4224
   422
    MatVecModM(A2p76, &Bg[3], &Bg[3], m2);
raj@4224
   423
  }
raj@4224
   424
  for (int i = 0; i < 6; ++i)
raj@4224
   425
    Cg[i] = Bg[i];
raj@364
   426
}
tomh@346
   427
tomh@346
   428
//-------------------------------------------------------------------------
tomh@346
   429
bool RngStream::SetPackageSeed (const uint32_t seed[6])
tomh@346
   430
{
raj@4224
   431
  if (!CheckSeed (seed))
raj@4224
   432
    {
raj@4224
   433
      return false;
raj@4224
   434
    }
tomh@346
   435
  for (int i = 0; i < 6; ++i)
raj@4234
   436
    nextSeed[i] = seed[i];
tomh@346
   437
  return true;
tomh@346
   438
}
raj@4224
   439
bool 
raj@4224
   440
RngStream::SetPackageSeed (uint32_t seed)
raj@4224
   441
{
raj@4224
   442
  uint32_t seeds[6] = {seed, seed, seed, seed, seed, seed};
raj@4224
   443
  return SetPackageSeed (seeds);
raj@4224
   444
}
raj@4224
   445
void 
raj@4224
   446
RngStream::GetPackageSeed (uint32_t seed[6])
raj@4224
   447
{
raj@4234
   448
  IntegerValue value;
raj@4234
   449
  g_rngSeed.GetValue (value);
raj@4234
   450
  uint32_t theSeed = value.Get ();
raj@4224
   451
  for (int i = 0; i < 6; i++)
raj@4224
   452
    {
raj@4234
   453
      seed[i] = static_cast<uint32_t> (theSeed);
raj@4224
   454
    }
raj@4224
   455
}
raj@4224
   456
void 
raj@4224
   457
RngStream::SetPackageRun (uint32_t run)
raj@4224
   458
{
raj@4224
   459
  g_rngRun.SetValue (IntegerValue (run));
raj@4224
   460
}
raj@4224
   461
uint32_t 
raj@4224
   462
RngStream::GetPackageRun (void)
raj@4224
   463
{
raj@4224
   464
  IntegerValue run;
raj@4224
   465
  g_rngRun.GetValue (run);
raj@4224
   466
  return run.Get ();
raj@4224
   467
}
raj@4224
   468
bool 
raj@4224
   469
RngStream::CheckSeed(uint32_t seed)
raj@4224
   470
{
raj@4224
   471
  uint32_t seeds[6] = {seed, seed, seed, seed, seed, seed};
raj@4224
   472
  return CheckSeed (seeds);
raj@4224
   473
}
raj@4224
   474
tomh@346
   475
tomh@346
   476
tomh@346
   477
//-------------------------------------------------------------------------
tomh@346
   478
bool RngStream::SetSeeds (const uint32_t seed[6])
tomh@346
   479
{
tomh@346
   480
  if (!CheckSeed (seed)) return false;
raj@4224
   481
  for (int i = 0; i < 6; ++i)
raj@4224
   482
    Cg[i] = Bg[i] = Ig[i] = seed[i];
raj@4224
   483
  return true;
tomh@346
   484
}
tomh@346
   485
tomh@346
   486
tomh@346
   487
//-------------------------------------------------------------------------
tomh@346
   488
// if e > 0, let n = 2^e + c;
tomh@346
   489
// if e < 0, let n = -2^(-e) + c;
tomh@346
   490
// if e = 0, let n = c.
tomh@346
   491
// Jump n steps forward if n > 0, backwards if n < 0.
tomh@346
   492
//
tomh@346
   493
void RngStream::AdvanceState (int32_t e, int32_t c)
tomh@346
   494
{
tomh@346
   495
    double B1[3][3], C1[3][3], B2[3][3], C2[3][3];
tomh@346
   496
tomh@346
   497
    if (e > 0) {
tomh@346
   498
        MatTwoPowModM (A1p0, B1, m1, e);
tomh@346
   499
        MatTwoPowModM (A2p0, B2, m2, e);
tomh@346
   500
    } else if (e < 0) {
tomh@346
   501
        MatTwoPowModM (InvA1, B1, m1, -e);
tomh@346
   502
        MatTwoPowModM (InvA2, B2, m2, -e);
tomh@346
   503
    }
tomh@346
   504
tomh@346
   505
    if (c >= 0) {
tomh@346
   506
        MatPowModM (A1p0, C1, m1, c);
tomh@346
   507
        MatPowModM (A2p0, C2, m2, c);
tomh@346
   508
    } else {
tomh@346
   509
        MatPowModM (InvA1, C1, m1, -c);
tomh@346
   510
        MatPowModM (InvA2, C2, m2, -c);
tomh@346
   511
    }
tomh@346
   512
tomh@346
   513
    if (e) {
tomh@346
   514
        MatMatModM (B1, C1, C1, m1);
tomh@346
   515
        MatMatModM (B2, C2, C2, m2);
tomh@346
   516
    }
tomh@346
   517
tomh@346
   518
    MatVecModM (C1, Cg, Cg, m1);
tomh@346
   519
    MatVecModM (C2, &Cg[3], &Cg[3], m2);
tomh@346
   520
}
tomh@346
   521
tomh@346
   522
tomh@346
   523
//-------------------------------------------------------------------------
tomh@346
   524
void RngStream::GetState (uint32_t seed[6]) const
tomh@346
   525
{
tomh@346
   526
   for (int i = 0; i < 6; ++i)
tomh@346
   527
      seed[i] = static_cast<uint32_t> (Cg[i]);
tomh@346
   528
}
tomh@346
   529
tomh@346
   530
tomh@346
   531
//-------------------------------------------------------------------------
tomh@346
   532
void RngStream::IncreasedPrecis (bool incp)
tomh@346
   533
{
tomh@346
   534
   incPrec = incp;
tomh@346
   535
}
tomh@346
   536
tomh@346
   537
tomh@346
   538
//-------------------------------------------------------------------------
tomh@346
   539
void RngStream::SetAntithetic (bool a)
tomh@346
   540
{
tomh@346
   541
   anti = a;
tomh@346
   542
}
tomh@346
   543
tomh@346
   544
tomh@346
   545
//-------------------------------------------------------------------------
tomh@346
   546
// Generate the next random number.
tomh@346
   547
//
tomh@346
   548
double RngStream::RandU01 ()
tomh@346
   549
{
tomh@346
   550
   if (incPrec)
tomh@346
   551
      return U01d();
tomh@346
   552
   else
tomh@346
   553
      return U01();
tomh@346
   554
}
tomh@346
   555
tomh@346
   556
tomh@346
   557
//-------------------------------------------------------------------------
tomh@346
   558
// Generate the next random integer.
tomh@346
   559
//
tomh@346
   560
int32_t RngStream::RandInt (int32_t low, int32_t high)
tomh@346
   561
{
tomh@346
   562
    return low + static_cast<int32_t> ((high - low + 1) * RandU01 ());
tomh@346
   563
};
tomh@346
   564
tomh@346
   565
} //namespace ns3