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