src/xbt/RngStream.c

   1 /***********************************************************************\
   2  *
   3  * File:           RngStream.c for multiple streams of Random Numbers
   4  * Language:       ANSI C
   5  * Copyright:      Pierre L'Ecuyer, University of Montreal
   6  * Date:           14 August 2001
   7  * License:      GPL version 2 or later
   8  *
   9  * Notice:         Please contact P. L'Ecuyer at <lecuyer@iro.UMontreal.ca>
  10  *                 for commercial purposes.
  11  *
  12 \***********************************************************************/
  13
  14
  15 #include "xbt/RngStream.h"
  16 #include <stdio.h>
  17 #include <stdlib.h>
  18 #include <string.h>
  19
  20
  21 /*---------------------------------------------------------------------*/
  22 /* Private part.                                                       */
  23 /*---------------------------------------------------------------------*/
  24
  25
  26 #define norm  2.328306549295727688e-10
  27 #define m1    4294967087.0
  28 #define m2    4294944443.0
  29 #define a12     1403580.0
  30 #define a13n     810728.0
  31 #define a21      527612.0
  32 #define a23n    1370589.0
  33
  34 #define two17   131072.0
  35 #define two53   9007199254740992.0
  36 #define fact  5.9604644775390625e-8    /* 1 / 2^24 */
  37
  38
  39
  40
  41 /* Default initial seed of the package. Will be updated to become
  42    the seed of the next created stream. */
  43 static double nextSeed[6] = { 12345, 12345, 12345, 12345, 12345, 12345 };
  44
  45
  46 /* The following are the transition matrices of the two MRG components */
  47 /* (in matrix form), raised to the powers -1, 1, 2^76, and 2^127, resp.*/
  48 static double InvA1[3][3] = {          /* Inverse of A1p0 */
  49           { 184888585.0,   0.0,  1945170933.0 },
  50           {         1.0,   0.0,           0.0 },
  51           {         0.0,   1.0,           0.0 }
  52           };
  53
  54 static double InvA2[3][3] = {          /* Inverse of A2p0 */
  55           {      0.0,  360363334.0,  4225571728.0 },
  56           {      1.0,          0.0,           0.0 },
  57           {      0.0,          1.0,           0.0 }
  58           };
  59
  60 static double A1p0[3][3] = {
  61           {       0.0,        1.0,       0.0 },
  62           {       0.0,        0.0,       1.0 },
  63           { -810728.0,  1403580.0,       0.0 }
  64           };
  65
  66 static double A2p0[3][3] = {
  67           {        0.0,        1.0,       0.0 },
  68           {        0.0,        0.0,       1.0 },
  69           { -1370589.0,        0.0,  527612.0 }
  70           };
  71
  72 static double A1p76[3][3] = {
  73           {      82758667.0, 1871391091.0, 4127413238.0 },
  74           {    3672831523.0,   69195019.0, 1871391091.0 },
  75           {    3672091415.0, 3528743235.0,   69195019.0 }
  76           };
  77
  78 static double A2p76[3][3] = {
  79           {    1511326704.0, 3759209742.0, 1610795712.0 },
  80           {    4292754251.0, 1511326704.0, 3889917532.0 },
  81           {    3859662829.0, 4292754251.0, 3708466080.0 }
  82           };
  83
  84 static double A1p127[3][3] = {
  85           {    2427906178.0, 3580155704.0,  949770784.0 },
  86           {     226153695.0, 1230515664.0, 3580155704.0 },
  87           {    1988835001.0,  986791581.0, 1230515664.0 }
  88           };
  89
  90 static double A2p127[3][3] = {
  91           {    1464411153.0,  277697599.0, 1610723613.0 },
  92           {      32183930.0, 1464411153.0, 1022607788.0 },
  93           {    2824425944.0,   32183930.0, 2093834863.0 }
  94           };
  95
  96
  97
  98
  99
 100 /*-------------------------------------------------------------------------*/
 101
 102 static double MultModM (double a, double s, double c, double m)
 103    /* Compute (a*s + c) % m. m must be < 2^35.  Works also for s, c < 0 */
 104 {
 105    double v;
 106    long a1;
 107    v = a * s + c;
 108    if ((v >= two53) || (v <= -two53)) {
 109       a1 = (long) (a / two17);
 110       a -= a1 * two17;
 111       v = a1 * s;
 112       a1 = (long) (v / m);
 113       v -= a1 * m;
 114       v = v * two17 + a * s + c;
 115    }
 116    a1 = (long) (v / m);
 117    if ((v -= a1 * m) < 0.0)
 118       return v += m;
 119    else
 120       return v;
 121 }
 122
 123
 124 /*-------------------------------------------------------------------------*/
 125
 126 static void MatVecModM (double A[3][3], double s[3], double v[3], double m)
 127    /* Returns v = A*s % m.  Assumes that -m < s[i] < m. */
 128    /* Works even if v = s. */
 129 {
 130    int i;
 131    double x[3];
 132    for (i = 0; i < 3; ++i) {
 133       x[i] = MultModM (A[i][0], s[0], 0.0, m);
 134       x[i] = MultModM (A[i][1], s[1], x[i], m);
 135       x[i] = MultModM (A[i][2], s[2], x[i], m);
 136    }
 137    for (i = 0; i < 3; ++i)
 138       v[i] = x[i];
 139 }
 140
 141
 142 /*-------------------------------------------------------------------------*/
 143
 144 static void MatMatModM (double A[3][3], double B[3][3], double C[3][3],
 145                         double m)
 146    /* Returns C = A*B % m. Work even if A = C or B = C or A = B = C. */
 147 {
 148    int i, j;
 149    double V[3], W[3][3];
 150    for (i = 0; i < 3; ++i) {
 151       for (j = 0; j < 3; ++j)
 152          V[j] = B[j][i];
 153       MatVecModM (A, V, V, m);
 154       for (j = 0; j < 3; ++j)
 155          W[j][i] = V[j];
 156    }
 157    for (i = 0; i < 3; ++i) {
 158       for (j = 0; j < 3; ++j)
 159          C[i][j] = W[i][j];
 160    }
 161 }
 162
 163
 164 /*-------------------------------------------------------------------------*/
 165
 166 static void MatTwoPowModM (double A[3][3], double B[3][3], double m, long e)
 167   /* Compute matrix B = (A^(2^e) % m);  works even if A = B */
 168 {
 169    int i, j;
 170
 171    /* initialize: B = A */
 172    if (A != B) {
 173       for (i = 0; i < 3; i++) {
 174          for (j = 0; j < 3; ++j)
 175             B[i][j] = A[i][j];
 176       }
 177    }
 178    /* Compute B = A^{2^e} */
 179    for (i = 0; i < e; i++)
 180       MatMatModM (B, B, B, m);
 181 }
 182
 183
 184 /*-------------------------------------------------------------------------*/
 185
 186 static void MatPowModM (double A[3][3], double B[3][3], double m, long n)
 187    /* Compute matrix B = A^n % m ;  works even if A = B */
 188 {
 189    int i, j;
 190    double W[3][3];
 191
 192    /* initialize: W = A; B = I */
 193    for (i = 0; i < 3; i++) {
 194       for (j = 0; j < 3; ++j) {
 195          W[i][j] = A[i][j];
 196          B[i][j] = 0.0;
 197       }
 198    }
 199    for (j = 0; j < 3; ++j)
 200       B[j][j] = 1.0;
 201
 202    /* Compute B = A^n % m using the binary decomposition of n */
 203    while (n > 0) {
 204       if (n % 2)
 205          MatMatModM (W, B, B, m);
 206       MatMatModM (W, W, W, m);
 207       n /= 2;
 208    }
 209 }
 210
 211
 212 /*-------------------------------------------------------------------------*/
 213
 214 static double U01 (RngStream g)
 215 {
 216    long k;
 217    double p1, p2, u;
 218
 219    /* Component 1 */
 220    p1 = a12 * g->Cg[1] - a13n * g->Cg[0];
 221    k = p1 / m1;
 222    p1 -= k * m1;
 223    if (p1 < 0.0)
 224       p1 += m1;
 225    g->Cg[0] = g->Cg[1];
 226    g->Cg[1] = g->Cg[2];
 227    g->Cg[2] = p1;
 228
 229    /* Component 2 */
 230    p2 = a21 * g->Cg[5] - a23n * g->Cg[3];
 231    k = p2 / m2;
 232    p2 -= k * m2;
 233    if (p2 < 0.0)
 234       p2 += m2;
 235    g->Cg[3] = g->Cg[4];
 236    g->Cg[4] = g->Cg[5];
 237    g->Cg[5] = p2;
 238
 239    /* Combination */
 240    u = ((p1 > p2) ? (p1 - p2) * norm : (p1 - p2 + m1) * norm);
 241    return (g->Anti) ? (1 - u) : u;
 242 }
 243
 244
 245 /*-------------------------------------------------------------------------*/
 246
 247 static double U01d (RngStream g)
 248 {
 249    double u;
 250    u = U01(g);
 251    if (g->Anti == 0) {
 252       u += U01(g) * fact;
 253       return (u < 1.0) ? u : (u - 1.0);
 254    } else {
 255       /* Don't forget that U01() returns 1 - u in the antithetic case */
 256       u += (U01(g) - 1.0) * fact;
 257       return (u < 0.0) ? u + 1.0 : u;
 258    }
 259 }
 260
 261
 262 /*-------------------------------------------------------------------------*/
 263
 264 static int CheckSeed (unsigned long seed[6])
 265 {
 266    /* Check that the seeds are legitimate values. Returns 0 if legal seeds,
 267      -1 otherwise */
 268    int i;
 269
 270    for (i = 0; i < 3; ++i) {
 271       if (seed[i] >= m1) {
 272    fprintf (stderr, "****************************************\n"
 273      "ERROR: Seed[%1d] >= m1, Seed is not set.\n"
 274      "****************************************\n\n", i);
 275    return (-1);
 276        }
 277    }
 278    for (i = 3; i < 6; ++i) {
 279       if (seed[i] >= m2) {
 280    fprintf (stderr, "****************************************\n"
 281      "ERROR: Seed[%1d] >= m1, Seed is not set.\n"
 282      "****************************************\n\n", i);
 283    return (-1);
 284        }
 285    }
 286    if (seed[0] == 0 && seed[1] == 0 && seed[2] == 0) {
 287       fprintf (stderr, "****************************\n"
 288         "ERROR: First 3 seeds = 0.\n"
 289         "****************************\n\n");
 290       return (-1);
 291    }
 292    if (seed[3] == 0 && seed[4] == 0 && seed[5] == 0) {
 293       fprintf (stderr, "****************************\n"
 294         "ERROR: Last 3 seeds = 0.\n"
 295         "****************************\n\n");
 296       return (-1);
 297    }
 298
 299    return 0;
 300 }
 301
 302
 303 /*---------------------------------------------------------------------*/
 304 /* Public part.                                                        */
 305 /*---------------------------------------------------------------------*/
 306
 307
 308 RngStream RngStream_CreateStream (const char name[])
 309 {
 310    int i;
 311    RngStream g;
 312    size_t len;
 313
 314    g = (RngStream) malloc (sizeof (struct RngStream_InfoState));
 315    if (g == NULL) {
 316       printf ("RngStream_CreateStream: No more memory\n\n");
 317       exit (EXIT_FAILURE);
 318    }
 319    if (name) {
 320       len = strlen (name);
 321       g->name = (char *) malloc ((len + 1) * sizeof (char));
 322       strncpy (g->name, name, len + 1);
 323    } else
 324       g->name = 0;
 325    g->Anti = 0;
 326    g->IncPrec = 0;
 327
 328    for (i = 0; i < 6; ++i) {
 329       g->Bg[i] = g->Cg[i] = g->Ig[i] = nextSeed[i];
 330    }
 331    MatVecModM (A1p127, nextSeed, nextSeed, m1);
 332    MatVecModM (A2p127, &nextSeed[3], &nextSeed[3], m2);
 333    return g;
 334 }
 335
 336 /*-------------------------------------------------------------------------*/
 337
 338 void RngStream_DeleteStream (RngStream * p)
 339 {
 340    if (*p == NULL)
 341       return;
 342    free((*p)->name);
 343    free (*p);
 344    *p = NULL;
 345 }
 346
 347 /*-------------------------------------------------------------------------*/
 348
 349 RngStream RngStream_CopyStream (const RngStream src)
 350 {
 351    RngStream g;
 352
 353    if(src == NULL) {
 354      printf ("RngStream_CopyStream: 'src' not initialized\n\n");
 355      exit (EXIT_FAILURE);
 356    }
 357
 358    g = (RngStream) malloc (sizeof (struct RngStream_InfoState));
 359    if (g == NULL) {
 360       printf ("RngStream_CopyStream: No more memory\n\n");
 361       exit (EXIT_FAILURE);
 362    }
 363    memcpy((void*) g, (void*) src, sizeof (struct RngStream_InfoState));
 364
 365    return g;
 366 }
 367
 368 /*-------------------------------------------------------------------------*/
 369
 370 void RngStream_ResetStartStream (RngStream g)
 371 {
 372    int i;
 373    for (i = 0; i < 6; ++i)
 374       g->Cg[i] = g->Bg[i] = g->Ig[i];
 375 }
 376
 377 /*-------------------------------------------------------------------------*/
 378
 379 void RngStream_ResetNextSubstream (RngStream g)
 380 {
 381    int i;
 382    MatVecModM (A1p76, g->Bg, g->Bg, m1);
 383    MatVecModM (A2p76, &g->Bg[3], &g->Bg[3], m2);
 384    for (i = 0; i < 6; ++i)
 385       g->Cg[i] = g->Bg[i];
 386 }
 387
 388 /*-------------------------------------------------------------------------*/
 389
 390 void RngStream_ResetStartSubstream (RngStream g)
 391 {
 392    int i;
 393    for (i = 0; i < 6; ++i)
 394       g->Cg[i] = g->Bg[i];
 395 }
 396
 397 /*-------------------------------------------------------------------------*/
 398
 399 int RngStream_SetPackageSeed (unsigned long seed[6])
 400 {
 401    int i;
 402    if (CheckSeed (seed))
 403       return -1;                    /* FAILURE */
 404    for (i = 0; i < 6; ++i)
 405       nextSeed[i] = seed[i];
 406    return 0;                       /* SUCCESS */
 407 }
 408
 409 /*-------------------------------------------------------------------------*/
 410
 411 int RngStream_SetSeed (RngStream g, unsigned long seed[6])
 412 {
 413    int i;
 414    if (CheckSeed (seed))
 415       return -1;                    /* FAILURE */
 416    for (i = 0; i < 6; ++i)
 417       g->Cg[i] = g->Bg[i] = g->Ig[i] = seed[i];
 418    return 0;                       /* SUCCESS */
 419 }
 420
 421 /*-------------------------------------------------------------------------*/
 422
 423 void RngStream_AdvanceState (RngStream g, long e, long c)
 424 {
 425    double B1[3][3], C1[3][3], B2[3][3], C2[3][3];
 426
 427    if (e > 0) {
 428       MatTwoPowModM (A1p0, B1, m1, e);
 429       MatTwoPowModM (A2p0, B2, m2, e);
 430    } else if (e < 0) {
 431       MatTwoPowModM (InvA1, B1, m1, -e);
 432       MatTwoPowModM (InvA2, B2, m2, -e);
 433    }
 434
 435    if (c >= 0) {
 436       MatPowModM (A1p0, C1, m1, c);
 437       MatPowModM (A2p0, C2, m2, c);
 438    } else {
 439       MatPowModM (InvA1, C1, m1, -c);
 440       MatPowModM (InvA2, C2, m2, -c);
 441    }
 442
 443    if (e) {
 444       MatMatModM (B1, C1, C1, m1);
 445       MatMatModM (B2, C2, C2, m2);
 446    }
 447
 448    MatVecModM (C1, g->Cg, g->Cg, m1);
 449    MatVecModM (C2, &g->Cg[3], &g->Cg[3], m2);
 450 }
 451
 452 /*-------------------------------------------------------------------------*/
 453
 454 void RngStream_GetState (RngStream g, unsigned long seed[6])
 455 {
 456    int i;
 457    for (i = 0; i < 6; ++i)
 458       seed[i] = g->Cg[i];
 459 }
 460
 461 /*-------------------------------------------------------------------------*/
 462
 463 void RngStream_WriteState (RngStream g)
 464 {
 465    int i;
 466    if (g == NULL)
 467       return;
 468    printf ("The current state of the Rngstream");
 469    if (g->name && (strlen (g->name) > 0))
 470       printf (" %s", g->name);
 471    printf (":\n   Cg = { ");
 472
 473    for (i = 0; i < 5; i++) {
 474       printf ("%lu, ", (unsigned long) g->Cg[i]);
 475    }
 476    printf ("%lu }\n\n", (unsigned long) g->Cg[5]);
 477 }
 478
 479 /*-------------------------------------------------------------------------*/
 480
 481 void RngStream_WriteStateFull (RngStream g)
 482 {
 483    int i;
 484    if (g == NULL)
 485       return;
 486    printf ("The RngStream");
 487    if (g->name && (strlen (g->name) > 0))
 488       printf (" %s", g->name);
 489    printf (":\n   Anti = %s\n", (g->Anti ? "true" : "false"));
 490    printf ("   IncPrec = %s\n", (g->IncPrec ? "true" : "false"));
 491
 492    printf ("   Ig = { ");
 493    for (i = 0; i < 5; i++) {
 494       printf ("%lu, ", (unsigned long) (g->Ig[i]));
 495    }
 496    printf ("%lu }\n", (unsigned long) g->Ig[5]);
 497
 498    printf ("   Bg = { ");
 499    for (i = 0; i < 5; i++) {
 500       printf ("%lu, ", (unsigned long) (g->Bg[i]));
 501    }
 502    printf ("%lu }\n", (unsigned long) g->Bg[5]);
 503
 504    printf ("   Cg = { ");
 505    for (i = 0; i < 5; i++) {
 506       printf ("%lu, ", (unsigned long) (g->Cg[i]));
 507    }
 508    printf ("%lu }\n\n", (unsigned long) g->Cg[5]);
 509 }
 510
 511 /*-------------------------------------------------------------------------*/
 512
 513 void RngStream_IncreasedPrecis (RngStream g, int incp)
 514 {
 515    g->IncPrec = incp;
 516 }
 517
 518 /*-------------------------------------------------------------------------*/
 519
 520 void RngStream_SetAntithetic (RngStream g, int a)
 521 {
 522    g->Anti = a;
 523 }
 524
 525 /*-------------------------------------------------------------------------*/
 526
 527 double RngStream_RandU01 (RngStream g)
 528 {
 529    if (g->IncPrec)
 530       return U01d (g);
 531    else
 532       return U01 (g);
 533 }
 534
 535 /*-------------------------------------------------------------------------*/
 536
 537 int RngStream_RandInt (RngStream g, int i, int j)
 538 {
 539    return i + (int) ((j - i + 1.0) * RngStream_RandU01 (g));
 540 }
 541
 542 /* Undefine this terms, or supernovae build will fail. */
 543 #undef norm
 544 #undef m1
 545 #undef m2
 546 #undef a12
 547 #undef a13n
 548 #undef a21
 549 #undef a23n
 550 #undef two17
 551 #undef two53
 552 #undef fact