开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > Windows编程 > 通讯编程 > 查看源码
rng.cc
资源名称：ns-allinone-2.33.tar.gz　[点击查看]
上传用户：rrhhcc
上传日期：2015-12-11
资源大小：54129k
文件大小：27k
源码类别：

通讯编程

开发平台：

Visual C++

rng.cc：源码内容
  1. /* -*- Mode:C++; c-basic-offset:8; tab-width:8; indent-tabs-mode:t -*- */
  2. /*
  3.  *  File:  RngStream.cc for multiple streams of Random Numbers
  4.  *  Copyright (C) 2001  Pierre L'Ecuyer (lecuyer@iro.umontreal.ca)
  5.  *
  6.  *  This program is free software; you can redistribute it and/or modify
  7.  *  it under the terms of the GNU General Public License as published by
  8.  *  the Free Software Foundation; either version 2 of the License, or
  9.  *  (at your option) any later version.
  10.  *
  11.  *  This program is distributed in the hope that it will be useful,
  12.  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13.  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14.  *  GNU General Public License for more details.
  15.  *
  16.  *  You should have received a copy of the GNU General Public License
  17.  *  along with this program; if not, write to the Free Software
  18.  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
  19.  *  02110-1301 USA
  20.  *
  21.  *  Linking this random number generator statically or dynamically with
  22.  *  other modules is making a combined work based on this random number
  23.  *  generator.  Thus, the terms and conditions of the GNU General Public
  24.  *  License cover the whole combination.
  25.  *
  26.  *  In addition, as a special exception, the copyright holders of this
  27.  *  random number generator give you permission to combine this random
  28.  *  number generator program with free software programs or libraries
  29.  *  that are released under the GNU LGPL and with code included in the
  30.  *  standard release of ns-2 under the Apache 2.0 license or under
  31.  *  otherwise-compatible licenses with advertising requirements (or
  32.  *  modified versions of such code, with unchanged license).  You may
  33.  *  copy and distribute such a system following the terms of the GNU GPL
  34.  *  for this random number generator and the licenses of the other code
  35.  *  concerned, provided  that you include the source code of that other
  36.  *  code when and as the GNU GPL requires distribution of source code.
  37.  *
  38.  *  Note that people who make modified versions of this random number
  39.  *  generator are not obligated to grant this special exception for
  40.  *  their modified versions; it is their choice whether to do so.
  41.  *  The GNU General Public License gives permission to release a
  42.  *  modified  version without this exception; this exception also makes
  43.  *  it possible to release a modified version which carries forward
  44.  *  this exception.
  45.  *
  46.  * Incorporated into rng.h and modified to maintain backward 
  47.  * compatibility with ns-2.1b8.  Users can use their current scripts 
  48.  * unmodified with the new RNG.  To get the same results as with the 
  49.  * previous RNG, define OLD_RNG when compiling (e.g., make -D OLD_RNG).
  50.  * - Michele Weigle, University of North Carolina (mcweigle@cs.unc.edu)
  51.  * October 10, 2001
  52.  *
  53.  */
  55. #ifndef lint
  56. static const char rcsid[] =
  57.     "@(#) $Header: /cvsroot/nsnam/nam-1/rng.cc,v 1.8 2007/02/12 07:18:09 tom_henderson Exp $ (LBL)";
  58. #endif
  60. /* new random number generator */
  61.   
  62. #ifndef WIN32
  63. #include <sys/time.h> // for gettimeofday
  64. #include <unistd.h>
  65. #endif
  67. #include <stdio.h>
  68. #ifndef OLD_RNG
  69. #include <string.h>
  70. #endif /* !OLD_RNG */
  71. #include "rng.h"
  73. #ifdef OLD_RNG
  75. /*
  76.  * RNGImplementation
  77.  */
  79. /*
  80.  * Generate a periodic sequence of pseudo-random numbers with
  81.  * a period of 2^31 - 2.  The generator is the "minimal standard"
  82.  * multiplicative linear congruential generator of Park, S.K. and
  83.  * Miller, K.W., "Random Number Generators: Good Ones are Hard to Find,"
  84.  * CACM 31:10, Oct. 88, pp. 1192-1201.
  85.  *
  86.  * The algorithm implemented is:  Sn = (a*s) mod m.
  87.  *   The modulus m can be approximately factored as: m = a*q + r,
  88.  *   where q = m div a and r = m mod a.
  89.  *
  90.  * Then Sn = g(s) + m*d(s)
  91.  *   where g(s) = a(s mod q) - r(s div q)
  92.  *   and d(s) = (s div q) - ((a*s) div m)
  93.  *
  94.  * Observations:
  95.  *   - d(s) is either 0 or 1.
  96.  *   - both terms of g(s) are in 0, 1, 2, . . ., m - 1.
  97.  *   - |g(s)| <= m - 1.
  98.  *   - if g(s) > 0, d(s) = 0, else d(s) = 1.
  99.  *   - s mod q = s - k*q, where k = s div q.
  100.  *
  101.  * Thus Sn = a(s - k*q) - r*k,
  102.  *   if (Sn <= 0), then Sn += m.
  103.  *
  104.  * To test an implementation for A = 16807, M = 2^31-1, you should
  105.  *   get the following sequences for the given starting seeds:
  106.  *
  107.  *   s0, s1,    s2,        s3,          . . . , s10000,     . . . , s551246 
  108.  *    1, 16807, 282475249, 1622650073,  . . . , 1043618065, . . . , 1003 
  109.  *    1973272912, 1207871363, 531082850, 967423018
  110.  *
  111.  * It is important to check for s10000 and s551246 with s0=1, to guard 
  112.  * against overflow.
  113. */
  115. /*
  116.  * The sparc assembly code [no longer here] is based on Carta, D.G., "Two Fast
  117.  * Implementations of the 'Minimal Standard' Random Number
  118.  * Generator," CACM 33:1, Jan. 90, pp. 87-88.
  119.  *
  120.  * ASSUME that "the product of two [signed 32-bit] integers (a, sn)
  121.  *        will occupy two [32-bit] registers (p, q)."
  122.  * Thus: a*s = (2^31)p + q
  123.  *
  124.  * From the observation that: x = y mod z is but
  125.  *   x = z * the fraction part of (y/z),
  126.  * Let: sn = m * Frac(as/m)
  127.  *
  128.  * For m = 2^31 - 1,
  129.  *   sn = (2^31 - 1) * Frac[as/(2^31 -1)]
  130.  *      = (2^31 - 1) * Frac[as(2^-31 + 2^-2(31) + 2^-3(31) + . . .)]
  131.  *      = (2^31 - 1) * Frac{[(2^31)p + q] [2^-31 + 2^-2(31) + 2^-3(31) + . . 
  132. .]}
  133.  *      = (2^31 - 1) * Frac[p+(p+q)2^-31+(p+q)2^-2(31)+(p+q)3^(-31)+ . . .]
  134.  *
  135.  * if p+q < 2^31:
  136.  *   sn = (2^31 - 1) * Frac[p + a fraction + a fraction + a fraction + . . .]
  137.  *      = (2^31 - 1) * [(p+q)2^-31 + (p+q)2^-2(31) + (p+q)3^(-31) + . . .]
  138.  *      = p + q
  139.  *
  140.  * otherwise:
  141.  *   sn = (2^31 - 1) * Frac[p + 1.frac . . .]
  142.  *      = (2^31 - 1) * (-1 + 1.frac . . .)
  143.  *      = (2^31 - 1) * [-1 + (p+q)2^-31 + (p+q)2^-2(31) + (p+q)3^(-31) + . . .]
  144.  *      = p + q - 2^31 + 1
  145.  */
  146.  
  148. #define A 16807L /* multiplier, 7**5 */
  149. #define M 2147483647L /* modulus, 2**31 - 1; both used in random */
  150. #define INVERSE_M ((double)4.656612875e-10) /* (1.0/(double)M) */
  152. long
  153. RNGImplementation::next()
  154. {
  155. long L, H;
  156. L = A * (seed_ & 0xffff);
  157. H = A * (seed_ >> 16);
  159. seed_ = ((H & 0x7fff) << 16) + L;
  160. seed_ -= 0x7fffffff;
  161. seed_ += H >> 15;
  163. if (seed_ <= 0) {
  164. seed_ += 0x7fffffff;
  165. }
  166. return(seed_);
  167. }
  169. double
  170. RNGImplementation::next_double()
  171. {
  172. long i = next();
  173. return i * INVERSE_M;
  174. }
  175. #endif /* OLD_RNG */
  177. /*
  178.  * RNG implements a nice front-end around RNGImplementation
  179.  */
  180. #ifndef stand_alone
  181. static class RNGClass : public TclClass {
  182. public:
  183. RNGClass() : TclClass("RNG") {}
  184. TclObject* create(int, const char*const*) {
  185. return(new RNG());
  186. }
  187. } class_rng;
  188. #endif /* stand_alone */
  190. /* default RNG */
  192. RNG* RNG::default_ = NULL;
  194. double
  195. RNG::normal(double avg, double std)
  196. {
  197.         static int parity = 0;
  198.         static double nextresult;
  199.         double sam1, sam2, rad;
  200.    
  201.         if (std == 0) return avg;
  202.         if (parity == 0) {
  203.                 sam1 = 2*uniform() - 1;
  204.                 sam2 = 2*uniform() - 1;
  205.                 while ((rad = sam1*sam1 + sam2*sam2) >= 1) {
  206.                         sam1 = 2*uniform() - 1;
  207.                         sam2 = 2*uniform() - 1;
  208.                 }
  209.                 rad = sqrt((-2*log(rad))/rad);
  210.                 nextresult = sam2 * rad;
  211.                 parity = 1;
  212.                 return (sam1 * rad * std + avg);
  213.         }
  214.         else {
  215.                 parity = 0;
  216.                 return (nextresult * std + avg);
  217.         }
  218. }
  220. #ifndef stand_alone
  221. int
  222. RNG::command(int argc, const char*const* argv)
  223. {
  224. Tcl& tcl = Tcl::instance();
  226. if (argc == 3) {
  227. if (strcmp(argv[1], "testint") == 0) {
  228. int n = atoi(argv[2]);
  229. tcl.resultf("%d", uniform(n));
  230. return (TCL_OK);
  231. }
  232. if (strcmp(argv[1], "testdouble") == 0) {
  233. double d = atof(argv[2]);
  234. tcl.resultf("%6e", uniform(d));
  235. return (TCL_OK);
  236. }
  237. if (strcmp(argv[1], "seed") == 0) {
  238. int s = atoi(argv[2]);
  239. // NEEDSWORK: should be a way to set seed to PRDEF_SEED_SORUCE
  240. if (s) {
  241. if (s <= 0 || (unsigned int)s >= MAXINT) {
  242. // MAXINT on freebsd is a ``black hole'' (it says at MAXINT).
  243. tcl.resultf("Setting random number seed to known bad value.");
  244. return TCL_ERROR;
  245. };
  246. set_seed(RAW_SEED_SOURCE, s);
  247. } else set_seed(HEURISTIC_SEED_SOURCE, 0);
  248. return(TCL_OK);
  249. }
  250. }else if (argc == 2) {
  251. if (strcmp(argv[1], "next-random") == 0) {
  252. tcl.resultf("%u", uniform_positive_int());
  253. return(TCL_OK);
  254. }
  255. if (strcmp(argv[1], "seed") == 0) {
  256. #ifdef OLD_RNG
  257. tcl.resultf("%u", stream_.seed());
  258. #else
  259. tcl.resultf("%u", seed());
  260. #endif /* OLD_RNG */
  261. return(TCL_OK);
  262. }
  263. #ifndef OLD_RNG
  264. if (strcmp (argv[1], "next-substream") == 0) {
  265. reset_next_substream();
  266. return (TCL_OK);
  267. }
  268. if (strcmp (argv[1], "all-seeds") == 0) {
  269. unsigned long seeds[6];
  270. get_state (seeds);
  271. tcl.resultf ("%lu %lu %lu %lu %lu %lu",
  272.      seeds[0], seeds[1], seeds[2], 
  273.      seeds[3], seeds[4], seeds[5]);
  274. return (TCL_OK);
  275. }
  276. if (strcmp (argv[1], "reset-start-substream") == 0) {
  277. reset_start_substream();
  278. return (TCL_OK);
  279. }
  280. #endif /* !OLD_RNG */
  281. if (strcmp(argv[1], "default") == 0) {
  282. default_ = this;
  283. return(TCL_OK);
  284. }
  285. #if 0
  286. if (strcmp(argv[1], "test") == 0) {
  287. if (test())
  288. tcl.resultf("RNG test failed");
  289. else
  290. tcl.resultf("RNG test passed");
  291. return(TCL_OK);
  292. }
  293. #endif
  294. } else if (argc == 4) {
  295.                 if (strcmp(argv[1], "seed") == 0) {
  296.                         int s = atoi(argv[3]);
  297.                         if (strcmp(argv[2], "raw") == 0) {
  298.                                 set_seed(RNG::RAW_SEED_SOURCE, s);
  299.                         } else if (strcmp(argv[2], "predef") == 0) {
  300.                                 set_seed(RNG::PREDEF_SEED_SOURCE, s);
  301.                                 // s is the index in predefined seed array
  302.                                 // 0 <= s < 64
  303.                         } else if (strcmp(argv[2], "heuristic") == 0) {
  304.                                 set_seed(RNG::HEURISTIC_SEED_SOURCE, 0);
  305.                         }
  306.                         return(TCL_OK);
  307.                 }
  308.                 if (strcmp(argv[1], "normal") == 0) {
  309.                         double avg = strtod(argv[2], NULL);
  310.                         double std = strtod(argv[3], NULL);
  311.                         tcl.resultf("%g", normal(avg, std));
  312.                         return (TCL_OK);
  313.                 }
  314.                 if (strcmp(argv[1], "lognormal") == 0) {
  315.                         double avg = strtod(argv[2], NULL);
  316.                         double std = strtod(argv[3], NULL);
  317.                         tcl.resultf("%g", lognormal(avg, std));
  318.                         return (TCL_OK);
  319.                 }
  320. }
  321. return(TclObject::command(argc, argv));
  322. }
  323. #endif /* stand_alone */
  325. void
  326. RNG::set_seed(RNGSources source, int seed)
  327. {
  328. /* The following predefined seeds are evenly spaced around
  329.  * the 2^31 cycle.  Each is approximately 33,000,000 elements
  330.  * apart.
  331.  */
  332. #define N_SEEDS_ 64
  333. static long predef_seeds[N_SEEDS_] = {
  334. 1973272912L,  188312339L, 1072664641L,  694388766L,
  335. 2009044369L,  934100682L, 1972392646L, 1936856304L,
  336. 1598189534L, 1822174485L, 1871883252L,  558746720L,
  337. 605846893L, 1384311643L, 2081634991L, 1644999263L,
  338. 773370613L,  358485174L, 1996632795L, 1000004583L,
  339. 1769370802L, 1895218768L,  186872697L, 1859168769L,
  340. 349544396L, 1996610406L,  222735214L, 1334983095L,
  341. 144443207L,  720236707L,  762772169L,  437720306L,
  342. 939612284L,  425414105L, 1998078925L,  981631283L,
  343. 1024155645L,  822780843L,  701857417L,  960703545L,
  344. 2101442385L, 2125204119L, 2041095833L,   89865291L,
  345. 898723423L, 1859531344L,  764283187L, 1349341884L,
  346. 678622600L,  778794064L, 1319566104L, 1277478588L,
  347. 538474442L,  683102175L,  999157082L,  985046914L,
  348. 722594620L, 1695858027L, 1700738670L, 1995749838L,
  349. 1147024708L,  346983590L,  565528207L,  513791680L
  350. };
  351. static long heuristic_sequence = 0;
  353. switch (source) {
  354. case RAW_SEED_SOURCE:
  355. if (seed <= 0 || (unsigned int)seed >= MAXINT) // Wei Ye
  356. abort();
  357. // use it as it is
  358. break;
  359. case PREDEF_SEED_SOURCE:
  360. if (seed < 0 || seed >= N_SEEDS_)
  361. abort();
  362. seed = predef_seeds[seed];
  363. break;
  364. case HEURISTIC_SEED_SOURCE:
  365. timeval tv;
  366. gettimeofday(&tv, 0);
  367. heuristic_sequence++;   // Always make sure we're different than last time.
  368. seed = (tv.tv_sec ^ tv.tv_usec ^ (heuristic_sequence << 8)) & 0x7fffffff;
  369. break;
  370. };
  371. // set it
  372. // NEEDSWORK: should we throw out known bad seeds?
  373. // (are there any?)
  374. if (seed < 0)
  375. seed = -seed;
  376. #ifdef OLD_RNG
  377. stream_.set_seed(seed);
  378. #else
  379. set_seed(seed);
  380. #endif /* OLD_RNG */
  382. // Toss away the first few values of heuristic seed.
  383. // In practice this makes sequential heuristic seeds
  384. // generate different first values.
  385. //
  386. // How many values to throw away should be the subject
  387. // of careful analysis.  Until then, I just throw away
  388. // ``a bunch''.  --johnh
  389. if (source == HEURISTIC_SEED_SOURCE) {
  390. int i;
  391. for (i = 0; i < 128; i++) {
  392. #ifdef OLD_RNG
  393. stream_.next();
  394. #else
  395. next();
  396. #endif /* OLD_RNG */
  397. };
  398. };
  399. }
  403. /*
  404.  * RNGTest:
  405.  * Make sure the RNG makes known values.
  406.  * Optionally, print out some stuff.
  407.  *
  408.  * Simple test program:
  409.  *
  410.  * gcc rng.cc -Drng_stand_alone -o rng_test -lm
  411.  *
  412.  */
  413. #ifdef rng_stand_alone
  414. int main() { RNGTest test; test.verbose(); }
  415. #endif
  417. #ifdef rng_test
  418. RNGTest::RNGTest()
  419. {
  420. RNG rng(RNG::RAW_SEED_SOURCE, 1L);
  421. int i;
  422. long r;
  424. for (i = 0; i < 10000; i++) 
  425. r = rng.uniform_positive_int();
  427. #ifdef OLD_RNG
  428. if (r != 1043618065L) {
  429. fprintf (stderr, "r (%lu) != 1043618065Ln", r);
  430. abort();
  431. }
  432. #else
  433.   if (r != 582989707L) {
  434.   fprintf (stderr, "r (%lu) != 582989707n", r);
  435.   abort();
  436.   }
  437. #endif /* OLD_RNG */
  439. for (i = 10000; i < 551246; i++)
  440. r = rng.uniform_positive_int();
  442. #ifdef OLD_RNG
  443. if (r != 1003L) {
  444. fprintf (stderr, "r (%lu) != 1003Ln", r);
  445. abort();
  446. }
  447. #else
  448.   if (r != 1625334359L) {
  449.   fprintf (stderr, "r (%lu) != 1625334359Ln", r);
  450.   abort();
  451.   }
  452. #endif /* OLD_RNG */
  453. }
  455. void
  456. RNGTest::first_n(RNG::RNGSources source, long seed, int n)
  457. {
  458. RNG rng(source, seed);
  460. for (int i = 0; i < n; i++) {
  461. long r = rng.uniform_positive_int();
  462. printf("%10lu ", r);
  463. };
  464. printf("n");
  465. }
  467. void
  468. RNGTest::verbose()
  469. {
  470. printf ("default: ");
  471. first_n(RNG::RAW_SEED_SOURCE, 1L, 5);
  473. int i;
  474. for (i = 0; i < 64; i++) {
  475. printf ("predef source %2u: ", i);
  476. first_n(RNG::PREDEF_SEED_SOURCE, i, 5);
  477. };
  479. printf("heuristic seeds should be different from each other and on each run.n");
  480. for (i = 0; i < 64; i++) {
  481. printf ("heuristic source %2u: ", i);
  482. first_n(RNG::HEURISTIC_SEED_SOURCE, i, 5);
  483. };
  484. }
  485. #endif /* rng_test */
  487. #ifndef OLD_RNG
  488. namespace 
  489. { 
  490. const double m1 = 4294967087.0; 
  491. const double m2 = 4294944443.0; 
  492. const double norm = 1.0 / (m1 + 1.0); 
  493. const double a12 = 1403580.0; 
  494. const double a13n = 810728.0; 
  495. const double a21 = 527612.0; 
  496. const double a23n = 1370589.0; 
  497. const double two17 = 131072.0; 
  498. const double two53 = 9007199254740992.0; 
  499. const double fact = 5.9604644775390625e-8; /* 1 / 2^24 */ 
  501. // The following are the transition matrices of the two MRG 
  502. // components (in matrix form), raised to the powers -1, 1, 
  503. // 2^76, and 2^127, resp. 
  505. const double InvA1[3][3] = { // Inverse of A1p0 
  506. { 184888585.0, 0.0, 1945170933.0 }, 
  507. { 1.0, 0.0, 0.0 }, 
  508. { 0.0, 1.0, 0.0 } 
  509. }; 
  511. const double InvA2[3][3] = { // Inverse of A2p0 
  512. { 0.0, 360363334.0, 4225571728.0 }, 
  513. { 1.0, 0.0, 0.0 }, 
  514. { 0.0, 1.0, 0.0 } 
  515. }; 
  517. const double A1p0[3][3] = { 
  518. { 0.0, 1.0, 0.0 }, 
  519. { 0.0, 0.0, 1.0 }, 
  520. { -810728.0, 1403580.0, 0.0 } 
  521. }; 
  523. const double A2p0[3][3] = { 
  524. { 0.0, 1.0, 0.0 }, 
  525. { 0.0, 0.0, 1.0 }, 
  526. { -1370589.0, 0.0, 527612.0 } 
  527. }; 
  529. const double A1p76[3][3] = { 
  530. { 82758667.0, 1871391091.0, 4127413238.0 }, 
  531. { 3672831523.0, 69195019.0, 1871391091.0 }, 
  532. { 3672091415.0, 3528743235.0, 69195019.0 } 
  533. }; 
  535. const double A2p76[3][3] = { 
  536. { 1511326704.0, 3759209742.0, 1610795712.0 }, 
  537. { 4292754251.0, 1511326704.0, 3889917532.0 }, 
  538. { 3859662829.0, 4292754251.0, 3708466080.0 } 
  539. }; 
  541. const double A1p127[3][3] = { 
  542. { 2427906178.0, 3580155704.0, 949770784.0 }, 
  543. { 226153695.0, 1230515664.0, 3580155704.0 }, 
  544. { 1988835001.0, 986791581.0, 1230515664.0 } 
  545. }; 
  547. const double A2p127[3][3] = { 
  548. { 1464411153.0, 277697599.0, 1610723613.0 }, 
  549. { 32183930.0, 1464411153.0, 1022607788.0 }, 
  550. { 2824425944.0, 32183930.0, 2093834863.0 } 
  551. }; 
  553. //------------------------------------------------------------------- 
  554. // Return (a*s + c) MOD m; a, s, c and m must be < 2^35 
  555. // 
  557. double MultModM (double a, double s, double c, double m) 
  558. { 
  559. double v; 
  560. long a1; 
  561. v=a*s+c; 
  563. if (v >= two53 || v <= -two53) { 
  564. a1 = static_cast<long> (a / two17); a -= a1 * two17; 
  565. v =a1*s; 
  566. a1 = static_cast<long> (v / m); v -= a1 * m; 
  567. v = v * two17 + a * s + c; 
  568. } 
  569. a1 = static_cast<long> (v / m); 
  570. /* in case v < 0)*/ 
  571. if ((v -= a1 * m) < 0.0) return v += m; else return v; 
  572. } 
  574. //------------------------------------------------------------------- 
  575. // Compute the vector v = A*s MOD m. Assume that -m < s[i] < m. 
  576. // Works also when v = s. 
  577. // 
  578. void MatVecModM (const double A[3][3], const double s[3], double v[3], 
  579.  double m) 
  580. { 
  581. int i; 
  582. double x[3]; // Necessary if v = s 
  583. for (i = 0; i < 3; ++i) { 
  584. x[i] = MultModM (A[i][0], s[0], 0.0, m); 
  585. x[i] = MultModM (A[i][1], s[1], x[i], m); 
  586. x[i] = MultModM (A[i][2], s[2], x[i], m); 
  587. } 
  588. for (i = 0; i < 3; ++i) 
  589. v[i] = x[i]; 
  590. } 
  592. //------------------------------------------------------------------- 
  593. // Compute the matrix C = A*B MOD m. Assume that -m < s[i] < m. 
  594. // Note: works also if A = C or B = C or A = B = C. 
  595. // 
  596. void MatMatModM (const double A[3][3], const double B[3][3], 
  597.  double C[3][3], double m) 
  598. { 
  599. int i, j; 
  600. double V[3], W[3][3]; 
  601. for (i = 0; i < 3; ++i) { 
  602. for (j = 0; j < 3; ++j) 
  603. V[j] = B[j][i]; 
  604. MatVecModM (A, V, V, m); 
  605. for (j = 0; j < 3; ++j) 
  607. W[j][i] = V[j]; 
  608. } 
  609. for (i = 0; i < 3; ++i) 
  610. for (j = 0; j < 3; ++j) 
  611. C[i][j] = W[i][j]; 
  612. } 
  614. //------------------------------------------------------------------- 
  615. // Compute the matrix B = (A^(2^e) Mod m); works also if A = B. 
  616. // 
  617. void MatTwoPowModM (const double A[3][3], double B[3][3], double m, 
  618.     long e) 
  619. { 
  620. int i, j; 
  621. /* initialize: B = A */ 
  622. if (A != B) { 
  623. for (i = 0; i < 3; ++i) 
  624. for (j = 0; j < 3; ++j) 
  625. B[i][j] = A[i][j]; 
  626. } 
  627. /* Compute B = A^(2^e) mod m */ 
  628. for (i = 0; i < e; i++) 
  629. MatMatModM (B, B, B, m); 
  630. } 
  632. //------------------------------------------------------------------- 
  633. // Compute the matrix B = (A^n Mod m); works even if A = B. 
  634. // 
  635. void MatPowModM (const double A[3][3], double B[3][3], double m, 
  636.  long n) 
  637. { 
  638. int i, j; 
  639. double W[3][3]; 
  640. /* initialize: W = A; B = I */ 
  641. for (i = 0; i < 3; ++i) 
  642. for (j = 0; j < 3; ++j) { 
  643. W[i][j] = A[i][j]; 
  644. B[i][j] = 0.0; 
  645. } 
  646. for (j = 0; j < 3; ++j) 
  647. B[j][j] = 1.0; 
  648. /* Compute B = A^n mod m using the binary decomposition of n */
  649. while (n > 0) { 
  650. if (n % 2) MatMatModM (W, B, B, m); 
  651. MatMatModM (W, W, W, m); 
  653. n/=2; 
  654. } 
  655. } 
  657. //-------------------------------------------------------------------- 
  658. // Check that the seeds are legitimate values. Returns 0 if legal 
  659. // seeds, -1 otherwise. 
  660. // 
  661. int CheckSeed (const unsigned long seed[6]) 
  662. { 
  663. int i; 
  664. for (i = 0; i < 3; ++i) { 
  665. if (seed[i] >= m1) { 
  666. fprintf (stderr, "****************************************nn");
  667. fprintf (stderr, "ERROR: Seed[%d] >= 4294967087, Seed is not set.", i);
  668. fprintf (stderr, "nn****************************************nn");
  669. return (-1); 
  670. } 
  671. } 
  672. for (i = 3; i < 6; ++i) { 
  673. if (seed[i] >= m2) { 
  674. fprintf (stderr, "****************************************nn");
  675. fprintf (stderr, "ERROR: Seed[%d] >= 429444443, Seed is not set.", i);
  676. fprintf (stderr, "nn****************************************nn");
  677. return (-1); 
  678. } 
  679. } 
  680. if (seed[0] == 0 && seed[1] == 0 && seed[2] == 0) { 
  681. fprintf (stderr, "****************************************nn");
  682. fprintf (stderr, "ERROR: First 3 seeds = 0.nn");
  683. fprintf (stderr, "****************************************nn");
  684. return (-1); 
  685. } 
  686. if (seed[3] == 0 && seed[4] == 0 && seed[5] == 0) { 
  687. fprintf (stderr, "****************************************nn");
  688. fprintf (stderr, "ERROR: Last 3 seeds = 0.nn");
  689. fprintf (stderr, "****************************************nn");
  690. return (-1); 
  691. } 
  692. return 0; 
  693. } 
  694. } // end of anonymous namespace 
  696. //------------------------------------------------------------------------- 
  697. // Generate the next random number. 
  698. // 
  699. double RNG::U01 () 
  700. { 
  701. long k; 
  702. double p1, p2, u; 
  703. /* Component 1 */ 
  704. p1 = a12 * Cg_[1] - a13n * Cg_[0]; 
  705. k = static_cast<long> (p1 / m1); 
  706. p1 -= k * m1; 
  707. if (p1 < 0.0) p1 += m1; 
  708. Cg_[0] = Cg_[1]; Cg_[1] = Cg_[2]; Cg_[2] = p1; 
  709. /* Component 2 */ 
  710. p2 = a21 * Cg_[5] - a23n * Cg_[3]; 
  711. k = static_cast<long> (p2 / m2); 
  712. p2 -= k * m2; 
  713. if (p2 < 0.0) p2 += m2; 
  714. Cg_[3] = Cg_[4]; Cg_[4] = Cg_[5]; Cg_[5] = p2; 
  715. /* Combination */ 
  716. u = ((p1 > p2) ? (p1 - p2) * norm : (p1 - p2 + m1) * norm); 
  717. return (anti_ == false) ? u : (1 - u); 
  718. } 
  720. //------------------------------------------------------------------------- 
  721. // Generate the next random number with extended (53 bits) precision. 
  722. // 
  723. double RNG::U01d () 
  724. { 
  725. double u; 
  726. u = U01(); 
  727. if (anti_) { 
  728. // Don't forget that U01() returns 1 - u in 
  729. // the antithetic case 
  730. u += (U01() - 1.0) * fact; 
  731. return (u < 0.0) ? u + 1.0 : u; 
  732. } else { 
  733. u += U01() * fact; 
  734. return (u < 1.0) ? u : (u - 1.0); 
  735. } 
  736. } 
  738. //************************************************************************* 
  739. // Public members of the class start here 
  740. //------------------------------------------------------------------------- 
  742. /*
  743.  * Functions for backward compatibility:
  744.  *
  745.  *   RNG (long seed)
  746.  *   void set_seed (long seed)
  747.  *   long seed()
  748.  *   long next()
  749.  *   double next_double()
  750.  */
  751. RNG::RNG (long seed) 
  752. {
  753. set_seed (seed);
  754. init();
  755. }
  757. void RNG::init()
  758. {
  759. anti_ = false; 
  760. inc_prec_ = false; 
  762. /* Information on a stream. The arrays {Cg_, Bg_, Ig_} contain the
  763.    current state of the stream, the starting state of the current
  764.    SubStream, and the starting state of the stream. This stream
  765.    generates antithetic variates if anti_ = true. It also generates
  766.    numbers with extended precision (53 bits if machine follows IEEE
  767.    754 standard) if inc_prec_ = true. next_seed_ will be the seed of
  768.    the next declared RngStream. */
  770. for (int i = 0; i < 6; ++i) { 
  771. Bg_[i] = Cg_[i] = Ig_[i] = next_seed_[i]; 
  772. } 
  773. MatVecModM (A1p127, next_seed_, next_seed_, m1); 
  774. MatVecModM (A2p127, &next_seed_[3], &next_seed_[3], m2); 
  775. }
  777. void RNG::set_seed (long seed) 
  778. {
  779. unsigned long seed_vec[6] = {0, 0, 0, 0, 0, 0};
  780. for (int i=0; i<6; i++) {
  781. seed_vec[i] = (unsigned long) seed;
  782. }
  783. set_package_seed (seed_vec);
  784. init();
  785. }
  787. long RNG::seed() 
  788. {
  789. unsigned long seed[6];
  790. get_state(seed);
  791. return ((long) seed[0]);
  792. }
  794. long RNG::next()
  795. {
  796. return (rand_int(0, MAXINT));
  797. }
  799. double RNG::next_double()
  800. {
  801. return (rand_u01());
  802. }
  803. /* End of backward compatibility functions */
  805. // The default seed of the package; will be the seed of the first 
  806. // declared RNG, unless set_package_seed is called. 
  807. // 
  808. double RNG::next_seed_[6] = 
  809. { 
  810. 12345.0, 12345.0, 12345.0, 12345.0, 12345.0, 12345.0 
  811. }; 
  813. //------------------------------------------------------------------------- 
  814. // constructor 
  815. // 
  816. RNG::RNG (const char *s) 
  817. { 
  818. if (strlen (s) > 99) {
  819. strncpy (name_, s, 99);
  820. name_[100] = 0;
  821. }
  822. else 
  823. strcpy (name_, s);
  825. init();
  826. }
  829. //------------------------------------------------------------------------- 
  830. // Reset Stream to beginning of Stream. 
  831. // 
  832. void RNG::reset_start_stream () 
  833. { 
  834. for (int i = 0; i < 6; ++i) 
  835. Cg_[i] = Bg_[i] = Ig_[i]; 
  836. } 
  838. //------------------------------------------------------------------------- 
  839. // Reset Stream to beginning of SubStream. 
  840. // 
  841. void RNG::reset_start_substream () 
  842. { 
  843. for (int i = 0; i < 6; ++i) 
  844. Cg_[i] = Bg_[i]; 
  845. } 
  847. //------------------------------------------------------------------------- 
  848. // Reset Stream to NextSubStream. 
  849. // 
  850. void RNG::reset_next_substream () 
  851. { 
  852. MatVecModM(A1p76, Bg_, Bg_, m1); 
  853. MatVecModM(A2p76, &Bg_[3], &Bg_[3], m2); 
  854. for (int i = 0; i < 6; ++i) 
  855. Cg_[i] = Bg_[i]; 
  856. } 
  858. //------------------------------------------------------------------------- 
  859. void RNG::set_package_seed (const unsigned long seed[6]) 
  860. { 
  861. if (CheckSeed (seed)) 
  862. abort();
  863. for (int i = 0; i < 6; ++i) 
  864. next_seed_[i] = seed[i]; 
  865. } 
  867. //------------------------------------------------------------------------- 
  868. void RNG::set_seed (const unsigned long seed[6]) 
  869. { 
  870. if (CheckSeed (seed)) 
  871. abort();
  872. for (int i = 0; i < 6; ++i) 
  873. Cg_[i] = Bg_[i] = Ig_[i] = seed[i]; 
  874. } 
  876. //------------------------------------------------------------------------- 
  877. // if e > 0, let n = 2^e + c; 
  878. // if e < 0, let n = -2^(-e) + c; 
  879. // if e = 0, let n = c. 
  880. // Jump n steps forward if n > 0, backwards if n < 0. 
  881. // 
  882. void RNG::advance_state (long e, long c) 
  883. { 
  884. double B1[3][3], C1[3][3], B2[3][3], C2[3][3]; 
  885. if (e > 0) { 
  886. MatTwoPowModM (A1p0, B1, m1, e); 
  887. MatTwoPowModM (A2p0, B2, m2, e); 
  888. } else if (e < 0) { 
  889. MatTwoPowModM (InvA1, B1, m1, -e); 
  890. MatTwoPowModM (InvA2, B2, m2, -e); 
  891. } 
  892. if (c >= 0) { 
  893. MatPowModM (A1p0, C1, m1, c); 
  894. MatPowModM (A2p0, C2, m2, c); 
  895. } else { 
  896. MatPowModM (InvA1, C1, m1, -c); 
  897. MatPowModM (InvA2, C2, m2, -c); 
  898. } 
  899. if (e) { 
  900. MatMatModM (B1, C1, C1, m1); 
  901. MatMatModM (B2, C2, C2, m2); 
  902. } 
  903. MatVecModM (C1, Cg_, Cg_, m1); 
  904. MatVecModM (C2, &Cg_[3], &Cg_[3], m2); 
  905. } 
  907. //------------------------------------------------------------------------- 
  908. void RNG::get_state (unsigned long seed[6]) const 
  909. { 
  910. for (int i = 0; i < 6; ++i) 
  911. seed[i] = static_cast<unsigned long> (Cg_[i]); 
  912. } 
  914. //------------------------------------------------------------------------- 
  915. void RNG::write_state () const 
  916. { 
  917. printf ("The current state of the Rngstream %s:n", name_);
  918. printf (" Cg_ = { ");
  919. for(int i=0;i<5;i++) { 
  920. printf ("%lu, ", (unsigned long) Cg_[i]);
  921. } 
  922. printf ("%lu }nn", (unsigned long) Cg_[5]);
  923. } 
  925. //------------------------------------------------------------------------- 
  926. void RNG::write_state_full () const 
  927. { 
  928. int i; 
  929. printf ("The RNG %s:n", name_);
  930. printf (" anti_ = %s", (anti_ ? "true" : "false")); 
  931. printf (" inc_prec_ = %sn", (inc_prec_ ? "true" : "false")); 
  933. printf (" Ig_ = { ");
  934. for (i = 0; i < 5; i++) { 
  935. printf ("%lu, ", (unsigned long) Ig_[i]);
  936. } 
  937. printf ("%lu }n", (unsigned long) Ig_[5]);
  939. printf (" Bg_ = { ");
  940. for (i = 0; i < 5; i++) { 
  941. printf ("%lu, ", (unsigned long) Bg_[i]);
  942. } 
  943. printf ("%lu }n", (unsigned long) Bg_[5]);
  945. printf (" Cg_ = { ");
  946. for (i = 0; i < 5; i++) { 
  947. printf ("%lu, ", (unsigned long) Cg_[i]);
  948. } 
  949. printf ("%lu }nn", (unsigned long) Cg_[5]);
  950. } 
  952. //------------------------------------------------------------------------- 
  953. void RNG::increased_precis (bool incp) 
  954. { 
  955. inc_prec_ = incp; 
  956. } 
  958. //------------------------------------------------------------------------- 
  959. void RNG::set_antithetic (bool a) 
  960. { 
  961. anti_ = a; 
  962. } 
  964. //------------------------------------------------------------------------- 
  965. // Generate the next random number. 
  966. // 
  967. double RNG::rand_u01 () 
  968. { 
  969. if (inc_prec_) 
  970. return U01d(); 
  971. else 
  972. return U01(); 
  973. } 
  975. //------------------------------------------------------------------------- 
  976. // Generate the next random integer. 
  977. // 
  978. long RNG::rand_int (long low, long high) 
  979. { 
  980. // return (long) low + (long) (((double) (high-low) * drn) + 0.5);
  981. return ((long) (low + (unsigned long) (((unsigned long) 
  982. (high-low+1)) * rand_u01())));
  983. }; 
  984. #endif /* !OLD_RNG */