1 /* Copyright (c) 2012, 2014. The SimGrid Team.
2 * All rights reserved. */
4 /* This program is free software; you can redistribute it and/or modify it
5 * under the terms of the license (GNU LGPL) which comes with this package. */
7 /***********************************************************************\
9 * File: RngStream.c for multiple streams of Random Numbers
11 * Copyright: Pierre L'Ecuyer, University of Montreal
12 * Date: 14 August 2001
13 * License: GPL version 2 or later
15 * Notice: Please contact P. L'Ecuyer at <lecuyer@iro.UMontreal.ca>
16 * for commercial purposes.
18 \***********************************************************************/
21 #include "xbt/RngStream.h"
22 #include "xbt/sysdep.h"
28 /*---------------------------------------------------------------------*/
30 /*---------------------------------------------------------------------*/
33 #define norm 2.328306549295727688e-10
34 #define m1 4294967087.0
35 #define m2 4294944443.0
39 #define a23n 1370589.0
41 #define two17 131072.0
42 #define two53 9007199254740992.0
43 #define fact 5.9604644775390625e-8 /* 1 / 2^24 */
48 /* Default initial seed of the package. Will be updated to become
49 the seed of the next created stream. */
50 static double nextSeed[6] = { 12345, 12345, 12345, 12345, 12345, 12345 };
53 /* The following are the transition matrices of the two MRG components */
54 /* (in matrix form), raised to the powers -1, 1, 2^76, and 2^127, resp.*/
55 static double InvA1[3][3] = { /* Inverse of A1p0 */
56 { 184888585.0, 0.0, 1945170933.0 },
61 static double InvA2[3][3] = { /* Inverse of A2p0 */
62 { 0.0, 360363334.0, 4225571728.0 },
67 static double A1p0[3][3] = {
70 { -810728.0, 1403580.0, 0.0 }
73 static double A2p0[3][3] = {
76 { -1370589.0, 0.0, 527612.0 }
79 static double A1p76[3][3] = {
80 { 82758667.0, 1871391091.0, 4127413238.0 },
81 { 3672831523.0, 69195019.0, 1871391091.0 },
82 { 3672091415.0, 3528743235.0, 69195019.0 }
85 static double A2p76[3][3] = {
86 { 1511326704.0, 3759209742.0, 1610795712.0 },
87 { 4292754251.0, 1511326704.0, 3889917532.0 },
88 { 3859662829.0, 4292754251.0, 3708466080.0 }
91 static double A1p127[3][3] = {
92 { 2427906178.0, 3580155704.0, 949770784.0 },
93 { 226153695.0, 1230515664.0, 3580155704.0 },
94 { 1988835001.0, 986791581.0, 1230515664.0 }
97 static double A2p127[3][3] = {
98 { 1464411153.0, 277697599.0, 1610723613.0 },
99 { 32183930.0, 1464411153.0, 1022607788.0 },
100 { 2824425944.0, 32183930.0, 2093834863.0 }
107 /*-------------------------------------------------------------------------*/
109 static double MultModM (double a, double s, double c, double m)
110 /* Compute (a*s + c) % m. m must be < 2^35. Works also for s, c < 0 */
115 if ((v >= two53) || (v <= -two53)) {
116 a1 = (long) (a / two17);
121 v = v * two17 + a * s + c;
124 if ((v -= a1 * m) < 0.0)
131 /*-------------------------------------------------------------------------*/
133 static void MatVecModM (double A[3][3], double s[3], double v[3], double m)
134 /* Returns v = A*s % m. Assumes that -m < s[i] < m. */
135 /* Works even if v = s. */
139 for (i = 0; i < 3; ++i) {
140 x[i] = MultModM (A[i][0], s[0], 0.0, m);
141 x[i] = MultModM (A[i][1], s[1], x[i], m);
142 x[i] = MultModM (A[i][2], s[2], x[i], m);
144 for (i = 0; i < 3; ++i)
149 /*-------------------------------------------------------------------------*/
151 static void MatMatModM (double A[3][3], double B[3][3], double C[3][3],
153 /* Returns C = A*B % m. Work even if A = C or B = C or A = B = C. */
156 double V[3], W[3][3];
157 for (i = 0; i < 3; ++i) {
158 for (j = 0; j < 3; ++j)
160 MatVecModM (A, V, V, m);
161 for (j = 0; j < 3; ++j)
164 for (i = 0; i < 3; ++i) {
165 for (j = 0; j < 3; ++j)
171 /*-------------------------------------------------------------------------*/
173 static void MatTwoPowModM (double A[3][3], double B[3][3], double m, long e)
174 /* Compute matrix B = (A^(2^e) % m); works even if A = B */
178 /* initialize: B = A */
180 for (i = 0; i < 3; i++) {
181 for (j = 0; j < 3; ++j)
185 /* Compute B = A^{2^e} */
186 for (i = 0; i < e; i++)
187 MatMatModM (B, B, B, m);
191 /*-------------------------------------------------------------------------*/
193 static void MatPowModM (double A[3][3], double B[3][3], double m, long n)
194 /* Compute matrix B = A^n % m ; works even if A = B */
199 /* initialize: W = A; B = I */
200 for (i = 0; i < 3; i++) {
201 for (j = 0; j < 3; ++j) {
206 for (j = 0; j < 3; ++j)
209 /* Compute B = A^n % m using the binary decomposition of n */
212 MatMatModM (W, B, B, m);
213 MatMatModM (W, W, W, m);
219 /*-------------------------------------------------------------------------*/
221 static double U01 (RngStream g)
227 p1 = a12 * g->Cg[1] - a13n * g->Cg[0];
237 p2 = a21 * g->Cg[5] - a23n * g->Cg[3];
247 u = ((p1 > p2) ? (p1 - p2) * norm : (p1 - p2 + m1) * norm);
248 return (g->Anti) ? (1 - u) : u;
252 /*-------------------------------------------------------------------------*/
254 static double U01d (RngStream g)
260 return (u < 1.0) ? u : (u - 1.0);
262 /* Don't forget that U01() returns 1 - u in the antithetic case */
263 u += (U01(g) - 1.0) * fact;
264 return (u < 0.0) ? u + 1.0 : u;
269 /*-------------------------------------------------------------------------*/
271 static int CheckSeed (unsigned long seed[6])
273 /* Check that the seeds are legitimate values. Returns 0 if legal seeds,
277 for (i = 0; i < 3; ++i) {
279 fprintf (stderr, "****************************************\n"
280 "ERROR: Seed[%1d] >= m1, Seed is not set.\n"
281 "****************************************\n\n", i);
285 for (i = 3; i < 6; ++i) {
287 fprintf (stderr, "****************************************\n"
288 "ERROR: Seed[%1d] >= m1, Seed is not set.\n"
289 "****************************************\n\n", i);
293 if (seed[0] == 0 && seed[1] == 0 && seed[2] == 0) {
294 fprintf (stderr, "****************************\n"
295 "ERROR: First 3 seeds = 0.\n"
296 "****************************\n\n");
299 if (seed[3] == 0 && seed[4] == 0 && seed[5] == 0) {
300 fprintf (stderr, "****************************\n"
301 "ERROR: Last 3 seeds = 0.\n"
302 "****************************\n\n");
310 /*---------------------------------------------------------------------*/
312 /*---------------------------------------------------------------------*/
315 RngStream RngStream_CreateStream (const char name[])
321 g = (RngStream) xbt_malloc (sizeof (struct RngStream_InfoState));
323 printf ("RngStream_CreateStream: No more memory\n\n");
328 g->name = (char *) xbt_malloc ((len + 1) * sizeof (char));
329 strncpy (g->name, name, len + 1);
335 for (i = 0; i < 6; ++i) {
336 g->Bg[i] = g->Cg[i] = g->Ig[i] = nextSeed[i];
338 MatVecModM (A1p127, nextSeed, nextSeed, m1);
339 MatVecModM (A2p127, &nextSeed[3], &nextSeed[3], m2);
343 /*-------------------------------------------------------------------------*/
345 void RngStream_DeleteStream (RngStream * p)
354 /*-------------------------------------------------------------------------*/
356 RngStream RngStream_CopyStream (const RngStream src)
361 printf ("RngStream_CopyStream: 'src' not initialized\n\n");
365 g = (RngStream) xbt_malloc (sizeof (struct RngStream_InfoState));
367 printf ("RngStream_CopyStream: No more memory\n\n");
370 memcpy((void*) g, (void*) src, sizeof (struct RngStream_InfoState));
375 /*-------------------------------------------------------------------------*/
377 void RngStream_ResetStartStream (RngStream g)
380 for (i = 0; i < 6; ++i)
381 g->Cg[i] = g->Bg[i] = g->Ig[i];
384 /*-------------------------------------------------------------------------*/
386 void RngStream_ResetNextSubstream (RngStream g)
389 MatVecModM (A1p76, g->Bg, g->Bg, m1);
390 MatVecModM (A2p76, &g->Bg[3], &g->Bg[3], m2);
391 for (i = 0; i < 6; ++i)
395 /*-------------------------------------------------------------------------*/
397 void RngStream_ResetStartSubstream (RngStream g)
400 for (i = 0; i < 6; ++i)
404 /*-------------------------------------------------------------------------*/
406 int RngStream_SetPackageSeed (unsigned long seed[6])
409 if (CheckSeed (seed))
410 return -1; /* FAILURE */
411 for (i = 0; i < 6; ++i)
412 nextSeed[i] = seed[i];
413 return 0; /* SUCCESS */
416 /*-------------------------------------------------------------------------*/
418 int RngStream_SetSeed (RngStream g, unsigned long seed[6])
421 if (CheckSeed (seed))
422 return -1; /* FAILURE */
423 for (i = 0; i < 6; ++i)
424 g->Cg[i] = g->Bg[i] = g->Ig[i] = seed[i];
425 return 0; /* SUCCESS */
428 /*-------------------------------------------------------------------------*/
430 void RngStream_AdvanceState (RngStream g, long e, long c)
432 double B1[3][3], C1[3][3], B2[3][3], C2[3][3];
435 MatTwoPowModM (A1p0, B1, m1, e);
436 MatTwoPowModM (A2p0, B2, m2, e);
438 MatTwoPowModM (InvA1, B1, m1, -e);
439 MatTwoPowModM (InvA2, B2, m2, -e);
443 MatPowModM (A1p0, C1, m1, c);
444 MatPowModM (A2p0, C2, m2, c);
446 MatPowModM (InvA1, C1, m1, -c);
447 MatPowModM (InvA2, C2, m2, -c);
451 MatMatModM (B1, C1, C1, m1);
452 MatMatModM (B2, C2, C2, m2);
455 MatVecModM (C1, g->Cg, g->Cg, m1);
456 MatVecModM (C2, &g->Cg[3], &g->Cg[3], m2);
459 /*-------------------------------------------------------------------------*/
461 void RngStream_GetState (RngStream g, unsigned long seed[6])
464 for (i = 0; i < 6; ++i)
468 /*-------------------------------------------------------------------------*/
470 void RngStream_WriteState (RngStream g)
475 printf ("The current state of the Rngstream");
476 if (g->name && (strlen (g->name) > 0))
477 printf (" %s", g->name);
478 printf (":\n Cg = { ");
480 for (i = 0; i < 5; i++) {
481 printf ("%lu, ", (unsigned long) g->Cg[i]);
483 printf ("%lu }\n\n", (unsigned long) g->Cg[5]);
486 /*-------------------------------------------------------------------------*/
488 void RngStream_WriteStateFull (RngStream g)
493 printf ("The RngStream");
494 if (g->name && (strlen (g->name) > 0))
495 printf (" %s", g->name);
496 printf (":\n Anti = %s\n", (g->Anti ? "true" : "false"));
497 printf (" IncPrec = %s\n", (g->IncPrec ? "true" : "false"));
500 for (i = 0; i < 5; i++) {
501 printf ("%lu, ", (unsigned long) (g->Ig[i]));
503 printf ("%lu }\n", (unsigned long) g->Ig[5]);
506 for (i = 0; i < 5; i++) {
507 printf ("%lu, ", (unsigned long) (g->Bg[i]));
509 printf ("%lu }\n", (unsigned long) g->Bg[5]);
512 for (i = 0; i < 5; i++) {
513 printf ("%lu, ", (unsigned long) (g->Cg[i]));
515 printf ("%lu }\n\n", (unsigned long) g->Cg[5]);
518 /*-------------------------------------------------------------------------*/
520 void RngStream_IncreasedPrecis (RngStream g, int incp)
525 /*-------------------------------------------------------------------------*/
527 void RngStream_SetAntithetic (RngStream g, int a)
532 /*-------------------------------------------------------------------------*/
534 double RngStream_RandU01 (RngStream g)
542 /*-------------------------------------------------------------------------*/
544 int RngStream_RandInt (RngStream g, int i, int j)
546 return i + (int) ((j - i + 1.0) * RngStream_RandU01 (g));