Logo AND Algorithmique Numérique Distribuée

Public GIT Repository
Use xbt_malloc (and friends) instead of direct malloc()
[simgrid.git] / 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