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

源码开发语言/平台
当前位置：首页 > 源码/资料 > 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.cc 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.,
  50.  *  //make -D OLD_RNG).
  51.  *  // - Michele Weigle, University of North Carolina
  52.  *  // (mcweigle@cs.unc.edu)  October 10, 2001
  53.  */
  55. #ifndef lint
  56. static const char rcsid[] =
  57. "@(#) $Header: /cvsroot/nsnam/ns-2/tools/rng.cc,v 1.31 2006/02/21 15:20:20 mahrenho Exp $ (LBL)";
  58. #endif
  60. /* new random number generator */
  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
  74. /*
  75.  * RNGImplementation
  76.  */
  78. /*
  79.  * Generate a periodic sequence of pseudo-random numbers with
  80.  * a period of 2^31 - 2.  The generator is the "minimal standard"
  81.  * multiplicative linear congruential generator of Park, S.K. and
  82.  * Miller, K.W., "Random Number Generators: Good Ones are Hard to Find,"
  83.  * CACM 31:10, Oct. 88, pp. 1192-1201.
  84.  *
  85.  * The algorithm implemented is:  Sn = (a*s) mod m.
  86.  *   The modulus m can be approximately factored as: m = a*q + r,
  87.  *   where q = m div a and r = m mod a.
  88.  *
  89.  * Then Sn = g(s) + m*d(s)
  90.  *   where g(s) = a(s mod q) - r(s div q)
  91.  *   and d(s) = (s div q) - ((a*s) div m)
  92.  *
  93.  * Observations:
  94.  *   - d(s) is either 0 or 1.
  95.  *   - both terms of g(s) are in 0, 1, 2, . . ., m - 1.
  96.  *   - |g(s)| <= m - 1.
  97.  *   - if g(s) > 0, d(s) = 0, else d(s) = 1.
  98.  *   - s mod q = s - k*q, where k = s div q.
  99.  *
  100.  * Thus Sn = a(s - k*q) - r*k,
  101.  *   if (Sn <= 0), then Sn += m.
  102.  *
  103.  * To test an implementation for A = 16807, M = 2^31-1, you should
  104.  *   get the following sequences for the given starting seeds:
  105.  *
  106.  *   s0, s1,    s2,        s3,          . . . , s10000,     . . . , s551246 
  107.  *    1, 16807, 282475249, 1622650073,  . . . , 1043618065, . . . , 1003 
  108.  *    1973272912, 1207871363, 531082850, 967423018
  109.  *
  110.  * It is important to check for s10000 and s551246 with s0=1, to guard 
  111.  * against overflow.
  112. */
  114. /*
  115.  * The sparc assembly code [no longer here] is based on Carta, D.G., "Two Fast
  116.  * Implementations of the 'Minimal Standard' Random Number
  117.  * Generator," CACM 33:1, Jan. 90, pp. 87-88.
  118.  *
  119.  * ASSUME that "the product of two [signed 32-bit] integers (a, sn)
  120.  *        will occupy two [32-bit] registers (p, q)."
  121.  * Thus: a*s = (2^31)p + q
  122.  *
  123.  * From the observation that: x = y mod z is but
  124.  *   x = z * the fraction part of (y/z),
  125.  * Let: sn = m * Frac(as/m)
  126.  *
  127.  * For m = 2^31 - 1,
  128.  *   sn = (2^31 - 1) * Frac[as/(2^31 -1)]
  129.  *      = (2^31 - 1) * Frac[as(2^-31 + 2^-2(31) + 2^-3(31) + . . .)]
  130.  *      = (2^31 - 1) * Frac{[(2^31)p + q] [2^-31 + 2^-2(31) + 2^-3(31) + . . 
  131. .]}
  132.  *      = (2^31 - 1) * Frac[p+(p+q)2^-31+(p+q)2^-2(31)+(p+q)3^(-31)+ . . .]
  133.  *
  134.  * if p+q < 2^31:
  135.  *   sn = (2^31 - 1) * Frac[p + a fraction + a fraction + a fraction + . . .]
  136.  *      = (2^31 - 1) * [(p+q)2^-31 + (p+q)2^-2(31) + (p+q)3^(-31) + . . .]
  137.  *      = p + q
  138.  *
  139.  * otherwise:
  140.  *   sn = (2^31 - 1) * Frac[p + 1.frac . . .]
  141.  *      = (2^31 - 1) * (-1 + 1.frac . . .)
  142.  *      = (2^31 - 1) * [-1 + (p+q)2^-31 + (p+q)2^-2(31) + (p+q)3^(-31) + . . .]
  143.  *      = p + q - 2^31 + 1
  144.  */
  145.  
  147. #define A 16807L /* multiplier, 7**5 */
  148. #define M 2147483647L /* modulus, 2**31 - 1; both used in random */
  149. #define INVERSE_M ((double)4.656612875e-10) /* (1.0/(double)M) */
  151. long
  152. RNGImplementation::next()
  153. {
  154. long L, H;
  155. L = A * (seed_ & 0xffff);
  156. H = A * (seed_ >> 16);
  158. seed_ = ((H & 0x7fff) << 16) + L;
  159. seed_ -= 0x7fffffff;
  160. seed_ += H >> 15;
  162. if (seed_ <= 0) {
  163. seed_ += 0x7fffffff;
  164. }
  165. return(seed_);
  166. }
  168. double
  169. RNGImplementation::next_double()
  170. {
  171. long i = next();
  172. return i * INVERSE_M;
  173. }
  174. #endif /* OLD_RNG */
  176. /*
  177.  * RNG implements a nice front-end around RNGImplementation
  178.  */
  179. #ifndef stand_alone
  180. static class RNGClass : public TclClass {
  181. public:
  182. RNGClass() : TclClass("RNG") {}
  183. TclObject* create(int, const char*const*) {
  184. return(new RNG());
  185. }
  186. } class_rng;
  187. #endif /* stand_alone */
  189. /* default RNG */
  191. RNG* RNG::default_ = NULL;
  193. double
  194. RNG::normal(double avg, double std)
  195. {
  196. static int parity = 0;
  197. static double nextresult;
  198. double sam1, sam2, rad;
  199.    
  200. if (std == 0) return avg;
  201. if (parity == 0) {
  202. sam1 = 2*uniform() - 1;
  203. sam2 = 2*uniform() - 1;
  204. while ((rad = sam1*sam1 + sam2*sam2) >= 1) {
  205. sam1 = 2*uniform() - 1;
  206. sam2 = 2*uniform() - 1;
  207. }
  208. rad = sqrt((-2*log(rad))/rad);
  209. nextresult = sam2 * rad;
  210. parity = 1;
  211. return (sam1 * rad * std + avg);
  212. }
  213. else {
  214. parity = 0;
  215. return (nextresult * std + avg);
  216. }
  217. }
  219. #ifndef stand_alone
  220. int
  221. RNG::command(int argc, const char*const* argv)
  222. {
  223. Tcl& tcl = Tcl::instance();
  225. if (argc == 3) {
  226. if (strcmp(argv[1], "testint") == 0) {
  227. int n = atoi(argv[2]);
  228. tcl.resultf("%d", uniform(n));
  229. return (TCL_OK);
  230. }
  231. if (strcmp(argv[1], "testdouble") == 0) {
  232. double d = atof(argv[2]);
  233. tcl.resultf("%6e", uniform(d));
  234. return (TCL_OK);
  235. }
  236. if (strcmp(argv[1], "seed") == 0) {
  237. int s = atoi(argv[2]);
  238. // NEEDSWORK: should be a way to set seed to PRDEF_SEED_SOURCE
  239. if (s) {
  240. if (s <= 0 || (unsigned int)s >= MAXINT) {
  241. tcl.resultf("Setting random number seed to known bad value.");
  242. return TCL_ERROR;
  243. };
  244. set_seed(RAW_SEED_SOURCE, s);
  245. } else set_seed(HEURISTIC_SEED_SOURCE, 0);
  246. return(TCL_OK);
  247. }
  248. } else if (argc == 2) {
  249. if (strcmp(argv[1], "next-random") == 0) {
  250. tcl.resultf("%u", uniform_positive_int());
  251. return(TCL_OK);
  252. }
  253. if (strcmp(argv[1], "seed") == 0) {
  254. #ifdef OLD_RNG
  255. tcl.resultf("%u", stream_.seed());
  256. #else
  257. tcl.resultf("%u", seed());
  258. #endif /* OLD_RNG */
  259. return(TCL_OK);
  260. }
  261. #ifndef OLD_RNG
  262. if (strcmp (argv[1], "next-substream") == 0) {
  263. reset_next_substream();
  264. return (TCL_OK);
  265. }
  266. if (strcmp (argv[1], "all-seeds") == 0) {
  267. unsigned long seeds[6];
  268. get_state (seeds);
  269. tcl.resultf ("%lu %lu %lu %lu %lu %lu",
  270.      seeds[0], seeds[1], seeds[2], 
  271.      seeds[3], seeds[4], seeds[5]);
  272. return (TCL_OK);
  273. }
  274. if (strcmp (argv[1], "reset-start-substream") == 0) {
  275. reset_start_substream();
  276. return (TCL_OK);
  277. }
  278. #endif /* !OLD_RNG */
  279. if (strcmp(argv[1], "default") == 0) {
  280. default_ = this;
  281. return(TCL_OK);
  282. }
  283. //#if 0
  284. if (strcmp(argv[1], "test") == 0) {
  285. RNGTest test; test.verbose_mil();
  286. //if (test())
  287. // tcl.resultf("RNG test failed");
  288. //else
  289. // tcl.resultf("RNG test passed");
  290. return(TCL_OK);
  291. }
  292. //#endif
  293. } else if (argc == 4) {
  294. if (strcmp(argv[1], "seed") == 0) {
  295. int s = atoi(argv[3]);
  296. if (strcmp(argv[2], "raw") == 0) {
  297. set_seed(RNG::RAW_SEED_SOURCE, s);
  298. } else if (strcmp(argv[2], "predef") == 0) {
  299. set_seed(RNG::PREDEF_SEED_SOURCE, s);
  300. // s is the index in predefined seed array
  301. // 0 <= s < 64
  302. } else if (strcmp(argv[2], "heuristic") == 0) {
  303. set_seed(RNG::HEURISTIC_SEED_SOURCE, 0);
  304. }
  305. return(TCL_OK);
  306. }
  307. if (strcmp(argv[1], "normal") == 0) {
  308. double avg = strtod(argv[2], NULL);
  309. double std = strtod(argv[3], NULL);
  310. tcl.resultf("%g", normal(avg, std));
  311. return (TCL_OK);
  312. }
  313. if (strcmp(argv[1], "lognormal") == 0) {
  314. double avg = strtod(argv[2], NULL);
  315. double std = strtod(argv[3], NULL);
  316. tcl.resultf("%g", lognormal(avg, std));
  317. return (TCL_OK);
  318. }
  319. }
  320. return(TclObject::command(argc, argv));
  321. }
  322. #endif /* stand_alone */
  324. void
  325. RNG::set_seed(RNGSources source, int seed)
  326. {
  327. /* The following predefined seeds are evenly spaced around
  328.  * the 2^31 cycle.  Each is approximately 33,000,000 elements
  329.  * apart.
  330.  */
  331. #define N_SEEDS_ 64
  332. static long predef_seeds[N_SEEDS_] = {
  333. 1973272912L,  188312339L, 1072664641L,  694388766L,
  334. 2009044369L,  934100682L, 1972392646L, 1936856304L,
  335. 1598189534L, 1822174485L, 1871883252L,  558746720L,
  336. 605846893L, 1384311643L, 2081634991L, 1644999263L,
  337. 773370613L,  358485174L, 1996632795L, 1000004583L,
  338. 1769370802L, 1895218768L,  186872697L, 1859168769L,
  339. 349544396L, 1996610406L,  222735214L, 1334983095L,
  340. 144443207L,  720236707L,  762772169L,  437720306L,
  341. 939612284L,  425414105L, 1998078925L,  981631283L,
  342. 1024155645L,  822780843L,  701857417L,  960703545L,
  343. 2101442385L, 2125204119L, 2041095833L,   89865291L,
  344. 898723423L, 1859531344L,  764283187L, 1349341884L,
  345. 678622600L,  778794064L, 1319566104L, 1277478588L,
  346. 538474442L,  683102175L,  999157082L,  985046914L,
  347. 722594620L, 1695858027L, 1700738670L, 1995749838L,
  348. 1147024708L,  346983590L,  565528207L,  513791680L
  349. };
  350. static long heuristic_sequence = 0;
  352. switch (source) {
  353. case RAW_SEED_SOURCE:
  354. if (seed <= 0 || (unsigned int)seed >= MAXINT)    // Wei Ye
  355. abort();
  356. // use it as it is
  357. break;
  358. case PREDEF_SEED_SOURCE:
  359. if (seed < 0 || seed >= N_SEEDS_)
  360. abort();
  361. seed = predef_seeds[seed];
  362. break;
  363. case HEURISTIC_SEED_SOURCE:
  364. timeval tv;
  365. gettimeofday(&tv, 0);
  366. heuristic_sequence++;   // Always make sure we're different than last time.
  367. seed = (tv.tv_sec ^ tv.tv_usec ^ (heuristic_sequence << 8)) & 0x7fffffff;
  368. break;
  369. };
  370. // set it
  371. // NEEDSWORK: should we throw out known bad seeds?
  372. // (are there any?)
  373. if (seed < 0)
  374. seed = -seed;
  375. #ifdef OLD_RNG
  376. stream_.set_seed(seed);
  377. #else
  378. set_seed(seed);
  379. #endif /* OLD_RNG */
  381. // Toss away the first few values of heuristic seed.
  382. // In practice this makes sequential heuristic seeds
  383. // generate different first values.
  384. //
  385. // How many values to throw away should be the subject
  386. // of careful analysis.  Until then, I just throw away
  387. // ``a bunch''.  --johnh
  388. if (source == HEURISTIC_SEED_SOURCE) {
  389. int i;
  390. for (i = 0; i < 128; i++) {
  391. #ifdef OLD_RNG
  392. stream_.next();
  393. #else
  394. next();
  395. #endif /* OLD_RNG */
  396. };
  397. };
  398. }
  401. /*
  402.  * RNGTest:
  403.  * Make sure the RNG makes known values.
  404.  * Optionally, print out some stuff.
  405.  *
  406.  * gcc rng.cc -Drng_stand_alone -o rng_test -lm
  407.  *
  408.  * Simple test program:
  409.  */
  410. #ifdef rng_stand_alone
  411. int main() { RNGTest test; test.verbose(); }
  412. #endif /* rng_stand_alone */
  414. #ifdef rng_test
  415. RNGTest::RNGTest()
  416. {
  417.   RNG rng(RNG::RAW_SEED_SOURCE, 1L);
  419.   int i; 
  420.   long r;
  421.   for (i = 0; i < 10000; i++) 
  422.   r = rng.uniform_positive_int();
  424. #ifdef OLD_RNG
  425.   if (r != 1043618065L) {
  426.   fprintf (stderr, "r (%lu) != 1043618065Ln", r);
  427.   exit(-1);
  428.   }
  429. #else
  430.   if (r != 1179983147L) {
  431.   fprintf (stderr, "r (%lu) != 1179983147Ln", r);
  432.   exit(-1);
  433.   }
  434. #endif /* OLD_RNG */
  436.   for (i = 10000; i < 551246; i++)
  437.   r = rng.uniform_positive_int();
  439. #ifdef OLD_RNG
  440.   if (r != 1003L) {
  441.   fprintf (stderr, "r (%lu) != 1003Ln", r);
  442.   exit(-1);
  443.   }
  444. #else
  445.   if (r != 817829295L) {
  446.   fprintf (stderr, "r (%lu) != 817829295Ln", r);
  447.   exit(-1);
  448.   }
  449. #endif /* OLD_RNG */
  451. }
  453. void
  454. RNGTest::first_n(RNG::RNGSources source, long seed, int n)
  455. {
  456. RNG rng(source, seed);
  458. for (int i = 0; i < n; i++) {
  459. long r = rng.uniform_positive_int();
  460. printf("%10lu ", r);
  461. };
  462. printf("n");
  463. }
  465. void
  466. RNGTest::first_n_mil(RNG::RNGSources source, long seed, int n, FILE *outfile)
  467. {
  468. RNG rng(source, seed);
  469. // print the first 1000 no, then every millionth upto n millions
  470. long m = n * 1000000;
  471. for (int i = 0; i < m; i++) {
  472. long r = rng.uniform_positive_int();
  473. if (i < 100 || (i % 1000000 == 0))
  474. fprintf(outfile, "%10lu ", r);
  475. };
  476. fprintf(outfile, "n");
  477. }
  479. void 
  480. RNGTest::verbose_mil()
  481. {
  482. FILE *outfile = NULL;
  483. outfile = fopen("temp.rands", "w");
  485. if (outfile == NULL)
  486. fprintf(stderr, "Cannot create temp.randsn");
  488. fprintf (outfile, "default: ");
  489. first_n_mil(RNG::RAW_SEED_SOURCE, 188312339, 5, outfile);
  490. int i = 1;
  491. fprintf (outfile, "predef source %2u: ", i);
  492. first_n_mil(RNG::PREDEF_SEED_SOURCE, i, 5, outfile);
  494. }
  496. void
  497. RNGTest::verbose()
  498. {
  499. printf ("default: ");
  500. first_n(RNG::RAW_SEED_SOURCE, 1L, 5);
  502. int i;
  503. for (i = 0; i < 64; i++) {
  504. printf ("predef source %2u: ", i);
  505. first_n(RNG::PREDEF_SEED_SOURCE, i, 5);
  506. };
  508. printf("heuristic seeds should be different from each other and on each run.n");
  509. for (i = 0; i < 64; i++) {
  510. printf ("heuristic source %2u: ", i);
  511. first_n(RNG::HEURISTIC_SEED_SOURCE, i, 5);
  512. };
  513. }
  515. #endif /* rng_test */
  517. #ifndef OLD_RNG
  518. using namespace std; 
  519. namespace 
  520. { 
  521. const double m1 = 4294967087.0; 
  522. const double m2 = 4294944443.0; 
  523. const double norm = 1.0 / (m1 + 1.0); 
  524. const double a12 = 1403580.0; 
  525. const double a13n = 810728.0; 
  526. const double a21 = 527612.0; 
  527. const double a23n = 1370589.0; 
  528. const double two17 = 131072.0; 
  529. const double two53 = 9007199254740992.0; 
  530. const double fact = 5.9604644775390625e-8; /* 1 / 2^24 */ 
  532. // The following are the transition matrices of the two MRG 
  533. // components (in matrix form), raised to the powers -1, 1, 
  534. // 2^76, and 2^127, resp. 
  536. const double InvA1[3][3] = { // Inverse of A1p0 
  537. { 184888585.0, 0.0, 1945170933.0 }, 
  538. { 1.0, 0.0, 0.0 }, 
  539. { 0.0, 1.0, 0.0 } 
  540. }; 
  542. const double InvA2[3][3] = { // Inverse of A2p0 
  543. { 0.0, 360363334.0, 4225571728.0 }, 
  544. { 1.0, 0.0, 0.0 }, 
  545. { 0.0, 1.0, 0.0 } 
  546. }; 
  548. const double A1p0[3][3] = { 
  549. { 0.0, 1.0, 0.0 }, 
  550. { 0.0, 0.0, 1.0 }, 
  551. { -810728.0, 1403580.0, 0.0 } 
  552. }; 
  554. const double A2p0[3][3] = { 
  555. { 0.0, 1.0, 0.0 }, 
  556. { 0.0, 0.0, 1.0 }, 
  557. { -1370589.0, 0.0, 527612.0 } 
  558. }; 
  560. const double A1p76[3][3] = { 
  561. { 82758667.0, 1871391091.0, 4127413238.0 }, 
  562. { 3672831523.0, 69195019.0, 1871391091.0 }, 
  563. { 3672091415.0, 3528743235.0, 69195019.0 } 
  564. }; 
  566. const double A2p76[3][3] = { 
  567. { 1511326704.0, 3759209742.0, 1610795712.0 }, 
  568. { 4292754251.0, 1511326704.0, 3889917532.0 }, 
  569. { 3859662829.0, 4292754251.0, 3708466080.0 } 
  570. }; 
  572. const double A1p127[3][3] = { 
  573. { 2427906178.0, 3580155704.0, 949770784.0 }, 
  574. { 226153695.0, 1230515664.0, 3580155704.0 }, 
  575. { 1988835001.0, 986791581.0, 1230515664.0 } 
  576. }; 
  578. const double A2p127[3][3] = { 
  579. { 1464411153.0, 277697599.0, 1610723613.0 }, 
  580. { 32183930.0, 1464411153.0, 1022607788.0 }, 
  581. { 2824425944.0, 32183930.0, 2093834863.0 } 
  582. }; 
  584. //------------------------------------------------------------------- 
  585. // Return (a*s + c) MOD m; a, s, c and m must be < 2^35 
  586. // 
  588. double MultModM (double a, double s, double c, double m) 
  589. { 
  590. double v; 
  591. long a1; 
  592. v=a*s+c; 
  594. if (v >= two53 || v <= -two53) { 
  595. a1 = static_cast<long> (a / two17); a -= a1 * two17; 
  596. v =a1*s; 
  597. a1 = static_cast<long> (v / m); v -= a1 * m; 
  598. v = v * two17 + a * s + c; 
  599. } 
  600. a1 = static_cast<long> (v / m); 
  601. /* in case v < 0)*/ 
  602. if ((v -= a1 * m) < 0.0) return v += m; else return v; 
  603. } 
  605. //------------------------------------------------------------------- 
  606. // Compute the vector v = A*s MOD m. Assume that -m < s[i] < m. 
  607. // Works also when v = s. 
  608. // 
  609. void MatVecModM (const double A[3][3], const double s[3], double v[3], 
  610.  double m) 
  611. { 
  612. int i; 
  613. double x[3]; // Necessary if v = s 
  614. for (i = 0; i < 3; ++i) { 
  615. x[i] = MultModM (A[i][0], s[0], 0.0, m); 
  616. x[i] = MultModM (A[i][1], s[1], x[i], m); 
  617. x[i] = MultModM (A[i][2], s[2], x[i], m); 
  618. } 
  619. for (i = 0; i < 3; ++i) 
  620. v[i] = x[i]; 
  621. } 
  623. //------------------------------------------------------------------- 
  624. // Compute the matrix C = A*B MOD m. Assume that -m < s[i] < m. 
  625. // Note: works also if A = C or B = C or A = B = C. 
  626. // 
  627. void MatMatModM (const double A[3][3], const double B[3][3], 
  628.  double C[3][3], double m) 
  629. { 
  630. int i, j; 
  631. double V[3], W[3][3]; 
  632. for (i = 0; i < 3; ++i) { 
  633. for (j = 0; j < 3; ++j) 
  634. V[j] = B[j][i]; 
  635. MatVecModM (A, V, V, m); 
  636. for (j = 0; j < 3; ++j) 
  638. W[j][i] = V[j]; 
  639. } 
  640. for (i = 0; i < 3; ++i) 
  641. for (j = 0; j < 3; ++j) 
  642. C[i][j] = W[i][j]; 
  643. } 
  645. //------------------------------------------------------------------- 
  646. // Compute the matrix B = (A^(2^e) Mod m); works also if A = B. 
  647. // 
  648. void MatTwoPowModM (const double A[3][3], double B[3][3], double m, 
  649.     long e) 
  650. { 
  651. int i, j; 
  652. /* initialize: B = A */ 
  653. if (A != B) { 
  654. for (i = 0; i < 3; ++i) 
  655. for (j = 0; j < 3; ++j) 
  656. B[i][j] = A[i][j]; 
  657. } 
  658. /* Compute B = A^(2^e) mod m */ 
  659. for (i = 0; i < e; i++) 
  660. MatMatModM (B, B, B, m); 
  661. } 
  663. //------------------------------------------------------------------- 
  664. // Compute the matrix B = (A^n Mod m); works even if A = B. 
  665. // 
  666. void MatPowModM (const double A[3][3], double B[3][3], double m, 
  667.  long n) 
  668. { 
  669. int i, j; 
  670. double W[3][3]; 
  671. /* initialize: W = A; B = I */ 
  672. for (i = 0; i < 3; ++i) 
  673. for (j = 0; j < 3; ++j) { 
  674. W[i][j] = A[i][j]; 
  675. B[i][j] = 0.0; 
  676. } 
  677. for (j = 0; j < 3; ++j) 
  678. B[j][j] = 1.0; 
  679. /* Compute B = A^n mod m using the binary decomposition of n */
  680. while (n > 0) { 
  681. if (n % 2) MatMatModM (W, B, B, m); 
  682. MatMatModM (W, W, W, m); 
  684. n/=2; 
  685. } 
  686. } 
  688. //-------------------------------------------------------------------- 
  689. // Check that the seeds are legitimate values. Returns 0 if legal 
  690. // seeds, -1 otherwise. 
  691. // 
  692. int CheckSeed (const unsigned long seed[6]) 
  693. { 
  694. int i; 
  695. for (i = 0; i < 3; ++i) { 
  696. if (seed[i] >= m1) { 
  697. fprintf (stderr, "****************************************nn");
  698. fprintf (stderr, "ERROR: Seed[%d] >= 4294967087, Seed is not set.", i);
  699. fprintf (stderr, "nn****************************************nn");
  700. return (-1); 
  701. } 
  702. } 
  703. for (i = 3; i < 6; ++i) { 
  704. if (seed[i] >= m2) { 
  705. fprintf (stderr, "****************************************nn");
  706. fprintf (stderr, "ERROR: Seed[%d] >= 429444443, Seed is not set.", i);
  707. fprintf (stderr, "nn****************************************nn");
  708. return (-1); 
  709. } 
  710. } 
  711. if (seed[0] == 0 && seed[1] == 0 && seed[2] == 0) { 
  712. fprintf (stderr, "****************************************nn");
  713. fprintf (stderr, "ERROR: First 3 seeds = 0.nn");
  714. fprintf (stderr, "****************************************nn");
  715. return (-1); 
  716. } 
  717. if (seed[3] == 0 && seed[4] == 0 && seed[5] == 0) { 
  718. fprintf (stderr, "****************************************nn");
  719. fprintf (stderr, "ERROR: Last 3 seeds = 0.nn");
  720. fprintf (stderr, "****************************************nn");
  721. return (-1); 
  722. } 
  723. return 0; 
  724. } 
  725. } // end of anonymous namespace 
  727. //------------------------------------------------------------------------- 
  728. // Generate the next random number. 
  729. // 
  730. double RNG::U01 () 
  731. { 
  732. long k; 
  733. double p1, p2, u; 
  734. /* Component 1 */ 
  735. p1 = a12 * Cg_[1] - a13n * Cg_[0]; 
  736. k = static_cast<long> (p1 / m1); 
  737. p1 -= k * m1; 
  738. if (p1 < 0.0) p1 += m1; 
  739. Cg_[0] = Cg_[1]; Cg_[1] = Cg_[2]; Cg_[2] = p1; 
  740. /* Component 2 */ 
  741. p2 = a21 * Cg_[5] - a23n * Cg_[3]; 
  742. k = static_cast<long> (p2 / m2); 
  743. p2 -= k * m2; 
  744. if (p2 < 0.0) p2 += m2; 
  745. Cg_[3] = Cg_[4]; Cg_[4] = Cg_[5]; Cg_[5] = p2; 
  746. /* Combination */ 
  747. u = ((p1 > p2) ? (p1 - p2) * norm : (p1 - p2 + m1) * norm); 
  748. return (anti_ == false) ? u : (1 - u); 
  749. } 
  751. //------------------------------------------------------------------------- 
  752. // Generate the next random number with extended (53 bits) precision. 
  753. // 
  754. double RNG::U01d () 
  755. { 
  756. double u; 
  757. u = U01(); 
  758. if (anti_) { 
  759. // Don't forget that U01() returns 1 - u in 
  760. // the antithetic case 
  761. u += (U01() - 1.0) * fact; 
  762. return (u < 0.0) ? u + 1.0 : u; 
  763. } else { 
  764. u += U01() * fact; 
  765. return (u < 1.0) ? u : (u - 1.0); 
  766. } 
  767. } 
  769. //************************************************************************* 
  770. // Public members of the class start here 
  771. //------------------------------------------------------------------------- 
  773. /*
  774.  * Functions for backward compatibility:
  775.  *
  776.  *   RNG (long seed)
  777.  *   void set_seed (long seed)
  778.  *   long seed()
  779.  *   long next()
  780.  *   double next_double()
  781.  */
  782. RNG::RNG (long seed) 
  783. {
  784. set_seed (seed);
  785. init();
  786. }
  788. void RNG::init()
  789. {
  790. anti_ = false; 
  791. inc_prec_ = false; 
  793. /* Information on a stream. The arrays {Cg_, Bg_, Ig_} contain the
  794.    current state of the stream, the starting state of the current
  795.    SubStream, and the starting state of the stream. This stream
  796.    generates antithetic variates if anti_ = true. It also generates
  797.    numbers with extended precision (53 bits if machine follows IEEE
  798.    754 standard) if inc_prec_ = true. next_seed_ will be the seed of
  799.    the next declared RngStream. */
  801. for (int i = 0; i < 6; ++i) { 
  802. Bg_[i] = Cg_[i] = Ig_[i] = next_seed_[i]; 
  803. } 
  804. MatVecModM (A1p127, next_seed_, next_seed_, m1); 
  805. MatVecModM (A2p127, &next_seed_[3], &next_seed_[3], m2); 
  806. }
  808. void RNG::set_seed (long seed) 
  809. {
  810. if (seed == 0) {
  811. set_seed (HEURISTIC_SEED_SOURCE, seed);
  812. }
  813. else {
  814. unsigned long seed_vec[6] = {0, 0, 0, 0, 0, 0};
  815. for (int i=0; i<6; i++) {
  816. seed_vec[i] = (unsigned long) seed;
  817. }
  818. set_package_seed (seed_vec);
  819. init();
  820. }
  821. }
  823. long RNG::seed() 
  824. {
  825. unsigned long seed[6];
  826. get_state(seed);
  827. return ((long) seed[0]);
  828. }
  830. long RNG::next()
  831. {
  832. return (rand_int(0, MAXINT));
  833. }
  835. double RNG::next_double()
  836. {
  837. return (rand_u01());
  838. }
  839. /* End of backward compatibility functions */
  841. // The default seed of the package; will be the seed of the first 
  842. // declared RNG, unless set_package_seed is called. 
  843. // 
  844. double RNG::next_seed_[6] = 
  845. { 
  846. 12345.0, 12345.0, 12345.0, 12345.0, 12345.0, 12345.0 
  847. }; 
  849. //------------------------------------------------------------------------- 
  850. // constructor 
  851. // 
  852. RNG::RNG (const char *s) 
  853. { 
  854. if (strlen (s) > 99) {
  855. strncpy (name_, s, 99);
  856. name_[100] = 0;
  857. }
  858. else 
  859. strcpy (name_, s);
  861. init();
  862. }
  865. //------------------------------------------------------------------------- 
  866. // Reset Stream to beginning of Stream. 
  867. // 
  868. void RNG::reset_start_stream () 
  869. { 
  870. for (int i = 0; i < 6; ++i) 
  871. Cg_[i] = Bg_[i] = Ig_[i]; 
  872. } 
  874. //------------------------------------------------------------------------- 
  875. // Reset Stream to beginning of SubStream. 
  876. // 
  877. void RNG::reset_start_substream () 
  878. { 
  879. for (int i = 0; i < 6; ++i) 
  880. Cg_[i] = Bg_[i]; 
  881. } 
  883. //------------------------------------------------------------------------- 
  884. // Reset Stream to NextSubStream. 
  885. // 
  886. void RNG::reset_next_substream () 
  887. { 
  888. MatVecModM(A1p76, Bg_, Bg_, m1); 
  889. MatVecModM(A2p76, &Bg_[3], &Bg_[3], m2); 
  890. for (int i = 0; i < 6; ++i) 
  891. Cg_[i] = Bg_[i]; 
  892. } 
  894. //------------------------------------------------------------------------- 
  895. void RNG::set_package_seed (const unsigned long seed[6]) 
  896. { 
  897. if (CheckSeed (seed)) 
  898. abort();
  899. for (int i = 0; i < 6; ++i) 
  900. next_seed_[i] = seed[i]; 
  901. } 
  903. //------------------------------------------------------------------------- 
  904. void RNG::set_seed (const unsigned long seed[6]) 
  905. { 
  906. if (CheckSeed (seed)) 
  907. abort();
  908. for (int i = 0; i < 6; ++i) 
  909. Cg_[i] = Bg_[i] = Ig_[i] = seed[i]; 
  910. } 
  912. //------------------------------------------------------------------------- 
  913. // if e > 0, let n = 2^e + c; 
  914. // if e < 0, let n = -2^(-e) + c; 
  915. // if e = 0, let n = c. 
  916. // Jump n steps forward if n > 0, backwards if n < 0. 
  917. // 
  918. void RNG::advance_state (long e, long c) 
  919. { 
  920. double B1[3][3], C1[3][3], B2[3][3], C2[3][3]; 
  921. if (e > 0) { 
  922. MatTwoPowModM (A1p0, B1, m1, e); 
  923. MatTwoPowModM (A2p0, B2, m2, e); 
  924. } else if (e < 0) { 
  925. MatTwoPowModM (InvA1, B1, m1, -e); 
  926. MatTwoPowModM (InvA2, B2, m2, -e); 
  927. } 
  928. if (c >= 0) { 
  929. MatPowModM (A1p0, C1, m1, c); 
  930. MatPowModM (A2p0, C2, m2, c); 
  931. } else { 
  932. MatPowModM (InvA1, C1, m1, -c); 
  933. MatPowModM (InvA2, C2, m2, -c); 
  934. } 
  935. if (e) { 
  936. MatMatModM (B1, C1, C1, m1); 
  937. MatMatModM (B2, C2, C2, m2); 
  938. } 
  939. MatVecModM (C1, Cg_, Cg_, m1); 
  940. MatVecModM (C2, &Cg_[3], &Cg_[3], m2); 
  941. } 
  943. //------------------------------------------------------------------------- 
  944. void RNG::get_state (unsigned long seed[6]) const 
  945. { 
  946. for (int i = 0; i < 6; ++i) 
  947. seed[i] = static_cast<unsigned long> (Cg_[i]); 
  948. } 
  950. //------------------------------------------------------------------------- 
  951. void RNG::write_state () const 
  952. { 
  953. printf ("The current state of the Rngstream %s:n", name_);
  954. printf (" Cg_ = { ");
  955. for(int i=0;i<5;i++) { 
  956. printf ("%lu, ", (unsigned long) Cg_[i]);
  957. } 
  958. printf ("%lu }nn", (unsigned long) Cg_[5]);
  959. } 
  961. //------------------------------------------------------------------------- 
  962. void RNG::write_state_full () const 
  963. { 
  964. int i; 
  965. printf ("The RNG %s:n", name_);
  966. printf (" anti_ = %s", (anti_ ? "true" : "false")); 
  967. printf (" inc_prec_ = %sn", (inc_prec_ ? "true" : "false")); 
  969. printf (" Ig_ = { ");
  970. for (i = 0; i < 5; i++) { 
  971. printf ("%lu, ", (unsigned long) Ig_[i]);
  972. } 
  973. printf ("%lu }n", (unsigned long) Ig_[5]);
  975. printf (" Bg_ = { ");
  976. for (i = 0; i < 5; i++) { 
  977. printf ("%lu, ", (unsigned long) Bg_[i]);
  978. } 
  979. printf ("%lu }n", (unsigned long) Bg_[5]);
  981. printf (" Cg_ = { ");
  982. for (i = 0; i < 5; i++) { 
  983. printf ("%lu, ", (unsigned long) Cg_[i]);
  984. } 
  985. printf ("%lu }nn", (unsigned long) Cg_[5]);
  986. } 
  988. //------------------------------------------------------------------------- 
  989. void RNG::increased_precis (bool incp) 
  990. { 
  991. inc_prec_ = incp; 
  992. } 
  994. //------------------------------------------------------------------------- 
  995. void RNG::set_antithetic (bool a) 
  996. { 
  997. anti_ = a; 
  998. } 
  1000. //------------------------------------------------------------------------- 
  1001. // Generate the next random number. 
  1002. // 
  1003. double RNG::rand_u01 () 
  1004. { 
  1005. if (inc_prec_) 
  1006. return U01d(); 
  1007. else 
  1008. return U01(); 
  1009. } 
  1011. //------------------------------------------------------------------------- 
  1012. // Generate the next random integer. 
  1013. // 
  1014. long RNG::rand_int (long low, long high) 
  1015. { 
  1016. // return (long) low + (long) (((double) (high-low) * drn) + 0.5);
  1017. return ((long) (low + (unsigned long) (((unsigned long) 
  1018. (high-low+1)) * rand_u01())));
  1019. } 
  1020. #endif /* !OLD_RNG */