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