Logo AND Algorithmique Numérique Distribuée

Public GIT Repository
Remove unused include "simgrid_config.h"
[simgrid.git] / src / rngstreams / 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 "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    if ((*p)->name != NULL)
343       free ((*p)->name);
344    free (*p);
345    *p = NULL;
346 }
347
348 /*-------------------------------------------------------------------------*/
349
350 void RngStream_ResetStartStream (RngStream g)
351 {
352    int i;
353    for (i = 0; i < 6; ++i)
354       g->Cg[i] = g->Bg[i] = g->Ig[i];
355 }
356
357 /*-------------------------------------------------------------------------*/
358
359 void RngStream_ResetNextSubstream (RngStream g)
360 {
361    int i;
362    MatVecModM (A1p76, g->Bg, g->Bg, m1);
363    MatVecModM (A2p76, &g->Bg[3], &g->Bg[3], m2);
364    for (i = 0; i < 6; ++i)
365       g->Cg[i] = g->Bg[i];
366 }
367
368 /*-------------------------------------------------------------------------*/
369
370 void RngStream_ResetStartSubstream (RngStream g)
371 {
372    int i;
373    for (i = 0; i < 6; ++i)
374       g->Cg[i] = g->Bg[i];
375 }
376
377 /*-------------------------------------------------------------------------*/
378
379 int RngStream_SetPackageSeed (unsigned long seed[6])
380 {
381    int i;
382    if (CheckSeed (seed))
383       return -1;                    /* FAILURE */
384    for (i = 0; i < 6; ++i)
385       nextSeed[i] = seed[i];
386    return 0;                       /* SUCCESS */
387 }
388
389 /*-------------------------------------------------------------------------*/
390
391 int RngStream_SetSeed (RngStream g, unsigned long seed[6])
392 {
393    int i;
394    if (CheckSeed (seed))
395       return -1;                    /* FAILURE */
396    for (i = 0; i < 6; ++i)
397       g->Cg[i] = g->Bg[i] = g->Ig[i] = seed[i];
398    return 0;                       /* SUCCESS */ 
399 }
400
401 /*-------------------------------------------------------------------------*/
402
403 void RngStream_AdvanceState (RngStream g, long e, long c)
404 {
405    double B1[3][3], C1[3][3], B2[3][3], C2[3][3];
406
407    if (e > 0) {
408       MatTwoPowModM (A1p0, B1, m1, e);
409       MatTwoPowModM (A2p0, B2, m2, e);
410    } else if (e < 0) {
411       MatTwoPowModM (InvA1, B1, m1, -e);
412       MatTwoPowModM (InvA2, B2, m2, -e);
413    }
414
415    if (c >= 0) {
416       MatPowModM (A1p0, C1, m1, c);
417       MatPowModM (A2p0, C2, m2, c);
418    } else {
419       MatPowModM (InvA1, C1, m1, -c);
420       MatPowModM (InvA2, C2, m2, -c);
421    }
422
423    if (e) {
424       MatMatModM (B1, C1, C1, m1);
425       MatMatModM (B2, C2, C2, m2);
426    }
427
428    MatVecModM (C1, g->Cg, g->Cg, m1);
429    MatVecModM (C2, &g->Cg[3], &g->Cg[3], m2);
430 }
431
432 /*-------------------------------------------------------------------------*/
433
434 void RngStream_GetState (RngStream g, unsigned long seed[6])
435 {
436    int i;
437    for (i = 0; i < 6; ++i)
438       seed[i] = g->Cg[i];
439 }
440
441 /*-------------------------------------------------------------------------*/
442
443 void RngStream_WriteState (RngStream g)
444 {
445    int i;
446    if (g == NULL)
447       return;
448    printf ("The current state of the Rngstream");
449    if (g->name && (strlen (g->name) > 0))
450       printf (" %s", g->name);
451    printf (":\n   Cg = { ");
452
453    for (i = 0; i < 5; i++) {
454       printf ("%lu, ", (unsigned long) g->Cg[i]);
455    }
456    printf ("%lu }\n\n", (unsigned long) g->Cg[5]);
457 }
458
459 /*-------------------------------------------------------------------------*/
460
461 void RngStream_WriteStateFull (RngStream g)
462 {
463    int i;
464    if (g == NULL)
465       return;
466    printf ("The RngStream");
467    if (g->name && (strlen (g->name) > 0))
468       printf (" %s", g->name);
469    printf (":\n   Anti = %s\n", (g->Anti ? "true" : "false"));
470    printf ("   IncPrec = %s\n", (g->IncPrec ? "true" : "false"));
471
472    printf ("   Ig = { ");
473    for (i = 0; i < 5; i++) {
474       printf ("%lu, ", (unsigned long) (g->Ig[i]));
475    }
476    printf ("%lu }\n", (unsigned long) g->Ig[5]);
477
478    printf ("   Bg = { ");
479    for (i = 0; i < 5; i++) {
480       printf ("%lu, ", (unsigned long) (g->Bg[i]));
481    }
482    printf ("%lu }\n", (unsigned long) g->Bg[5]);
483
484    printf ("   Cg = { ");
485    for (i = 0; i < 5; i++) {
486       printf ("%lu, ", (unsigned long) (g->Cg[i]));
487    }
488    printf ("%lu }\n\n", (unsigned long) g->Cg[5]);
489 }
490
491 /*-------------------------------------------------------------------------*/
492
493 void RngStream_IncreasedPrecis (RngStream g, int incp)
494 {
495    g->IncPrec = incp;
496 }
497
498 /*-------------------------------------------------------------------------*/
499
500 void RngStream_SetAntithetic (RngStream g, int a)
501 {
502    g->Anti = a;
503 }
504
505 /*-------------------------------------------------------------------------*/
506
507 double RngStream_RandU01 (RngStream g)
508 {
509    if (g->IncPrec)
510       return U01d (g);
511    else
512       return U01 (g);
513 }
514
515 /*-------------------------------------------------------------------------*/
516
517 int RngStream_RandInt (RngStream g, int i, int j)
518 {
519    return i + (int) ((j - i + 1.0) * RngStream_RandU01 (g));
520 }