Logo AND Algorithmique Numérique Distribuée

Public GIT Repository
rename the plugins from the command line, and document it
[simgrid.git] / src / xbt / RngStream.c
1 /* Copyright (c) 2012-2018. 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 }