src/xbt/RngStream.c

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