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