3 * Copyright (C) 2010 - DIGITEO - Michael Baudin
4 * Copyright (C) 2004 - Bruno Pincon
5 * Copyright (C) 1998 - Terry H. Andres
6 * Copyright (C) 1998 - Pierre Lecuyer
9 * clcg4 generator stuff
12 * The following code is from L'Ecuyer and Andres "A Randow Number based
13 * on the combinaison of Four LCG" (distributed at the Pierre L'Ecuyer
14 * home page with a corresponding paper).
17 * The original code was slightly modified by Bruno Pincon for inclusion
20 * list of main modifs :
22 * - lot of routine 's names have changed to have some kind of
23 * uniformity with the others generators
25 * - add a var is_init so that initialisation is performed inside
26 * this module (to simplify the interface). And bring modif in
27 * the different routines :
28 * if (!is_init) then proceed to initialisation ...
30 * - add a routine advance_state_clcg4 (for compatibility with the
31 * old package (Scilab used this feature))
33 * - I have change the generator (clcg4 routine) so as it has the
34 * form (1) in place of (2) (see the joined paper of L'Ecuyer &
39 * x_{j,n} = a_j * x_{j,n-1} mod m_j 0 <= j <= 3
41 * The output with form (2) (original form in this code) :
43 * z_n = ( sum_j delta_j * x_{j,n} / m_j ) mod 1
45 * have been changed in the form (1) :
47 * z_n = ( sum_j delta_j * x_{j,n} ) mod m_1 (then u_n = z_n / m_1)
49 * to have some "uniformity" with all the others generators (which
50 * gives integers). Also it is better for the uin(a,b) generation
51 * to start from integers.
55 /*---------------------------------------------------------------------*/
56 /* clcg4.c Implementation module */
57 /*---------------------------------------------------------------------*/
60 #include <math.h> /* for floor */
62 #include "others_generators.h"
63 #include "localization.h"
64 #include "warningmode.h"
66 /***********************************************************************/
68 /***********************************************************************/
70 #define H 32768 /* = 2^15 : use in MultModM. */
72 static int aw[4], avw[4], /* a[j]^{2^w} et a[j]^{2^{v+w}}. */
73 a[4] = { 45991, 207707, 138556, 49689 },
74 m[4] = { 2147483647, 2147483543, 2147483423, 2147483323 };
76 static int Ig[4][Maxgen+1], Lg[4][Maxgen+1], Cg[4][Maxgen+1];
77 /* Initial seed, previous seed, and current seed. */
80 static int is_init = 0;
81 static int v_default = 31;
82 static int w_default = 41;
85 static int MultModM (int s, int t, int M)
86 /* Returns (s*t) MOD M. Assumes that -M < s < M and -M < t < M. */
87 /* See L'Ecuyer and Cote (1991). */
89 int R, S0, S1, q, qh, rh, k;
93 if (s < H) { S0 = s; R = 0; }
96 S1 = s/H; S0 = s - H*S1;
97 qh = M/H; rh = M - H*qh;
100 S1 -= H; k = t/qh; R = H * (t - k*qh) - k*rh;
101 while (R < 0) R += M;
106 q = M/S1; k = t/q; R -= k * (M - S1*q);
109 while (R < 0) R += M;
111 k = R/qh; R = H * (R - k*qh) - k*rh;
112 while (R < 0) R += M;
116 q = M/S0; k = t/q; R -= k* (M - S0*q);
119 while (R < 0) R += M;
124 static void comp_aw_and_avw(int v, int w)
127 for (j = 0; j < 4; j++)
130 for (i = 1; i <= w; i++)
131 aw [j] = MultModM (aw [j], aw [j], m[j]);
133 for (i = 1; i <= v; i++)
134 avw [j] = MultModM (avw [j], avw [j], m[j]);
138 static void init_clcg4(int v, int w)
140 /* currently the scilab interface don't let the user chooses
141 * v and w (always v_default and w_default) so this routine
142 * is in the "private" part (also because initialisation is
143 * always perform inside this module, depending of the var
146 double sd[4] = {11111111., 22222222., 33333333., 44444444.};
147 comp_aw_and_avw(v, w);
148 set_initial_seed_clcg4(sd[0], sd[1], sd[2], sd[3]);
151 static int verif_seeds_clcg4(double s0, double s1, double s2, double s3)
153 /* verify that the seeds are "integers" and are in the good range */
154 if ( s0 == floor(s0) && s1 == floor(s1) &&
155 s2 == floor(s2) && s3 == floor(s3) &&
156 1 <= s0 && s0 <= 2147483646 &&
157 1 <= s1 && s1 <= 2147483542 &&
158 1 <= s2 && s2 <= 2147483422 &&
159 1 <= s3 && s3 <= 2147483322 )
165 static void display_info_clcg4(void)
167 /* display the seeds range (in case of error) */
168 sciprint(_("\n bad seeds for clcg4, must be integers with s1 in [1, 2147483646]\n s2 in [1, 2147483542]\n s3 in [1, 2147483422]\n s4 in [1, 2147483322]"));
172 /*---------------------------------------------------------------------*/
174 /*---------------------------------------------------------------------*/
177 int set_seed_clcg4(int g, double s0, double s1, double s2, double s3)
179 if (! is_init ) {init_clcg4(v_default,w_default); is_init = 1; };
181 if ( verif_seeds_clcg4(s0, s1, s2, s3) )
183 Ig [0][g] = (int) s0; Ig [1][g] = (int) s1;
184 Ig [2][g] = (int) s2; Ig [3][g] = (int) s3;
185 init_generator_clcg4(g, InitialSeed);
186 if (getWarningMode())
188 sciprint(_("WARNING: %s\n"), _("be aware that you have may lost synchronization\n between the virtual gen %d and the others !\n use grand(\"setall\", s1, s2, s3, s4) if you want recover it."), g);
194 display_info_clcg4();
199 void get_state_clcg4(int g, double s[4])
202 if (! is_init ) {init_clcg4(v_default,w_default); is_init = 1; };
203 for (j = 0; j < 4; j++) s [j] = (double) Cg [j][g];
206 void init_generator_clcg4(int g, SeedType Where)
209 if (! is_init ) {init_clcg4(v_default,w_default); is_init = 1; };
210 for (j = 0; j < 4; j++)
215 Lg [j][g] = Ig [j][g]; break;
217 Lg [j][g] = MultModM (aw [j], Lg [j][g], m [j]); break;
221 Cg [j][g] = Lg [j][g];
225 void advance_state_clcg4(int g, int k)
230 if (! is_init ) {init_clcg4(v_default,w_default); is_init = 1; };
232 for ( j = 0 ; j < 4 ; j++ )
235 for ( i = 1 ; i <= k ; i++ )
236 b[j] = MultModM( b[j], b[j], m[j]);
237 Ig[j][g] = MultModM ( b[j], Cg[j][g], m[j] );
239 init_generator_clcg4(g, InitialSeed);
242 int set_initial_seed_clcg4(double s0, double s1, double s2, double s3)
246 if (! is_init ) comp_aw_and_avw(v_default,w_default);
248 if ( ! verif_seeds_clcg4(s0, s1, s2, s3) )
250 display_info_clcg4();
255 Ig [0][0] = (int) s0;
256 Ig [1][0] = (int) s1;
257 Ig [2][0] = (int) s2;
258 Ig [3][0] = (int) s3;
259 init_generator_clcg4(0, InitialSeed);
260 for (g = 1; g <= Maxgen; g++)
262 for (j = 0; j < 4; j++)
263 Ig [j][g] = MultModM (avw [j], Ig [j][g-1], m [j]);
264 init_generator_clcg4(g, InitialSeed);
269 unsigned int clcg4(int g)
271 /* Modif Bruno : the generator have now the form (1) in place of (2) */
276 if (! is_init ) {init_clcg4(v_default,w_default); is_init = 1; };
278 /* advance the 4 LCG */
279 s = Cg [0][g]; k = s / 46693;
280 s = 45991 * (s - k * 46693) - k * 25884;
281 if (s < 0) s = s + 2147483647; Cg [0][g] = s;
283 s = Cg [1][g]; k = s / 10339;
284 s = 207707 * (s - k * 10339) - k * 870;
285 if (s < 0) s = s + 2147483543; Cg [1][g] = s;
287 s = Cg [2][g]; k = s / 15499;
288 s = 138556 * (s - k * 15499) - k * 3979;
289 if (s < 0) s = s + 2147483423; Cg [2][g] = s;
291 s = Cg [3][g]; k = s / 43218;
292 s = 49689 * (s - k * 43218) - k * 24121;
293 if (s < 0) s = s + 2147483323; Cg [3][g] = s;
296 u = (double)(Cg[0][g] - Cg[1][g]) + (double)(Cg[2][g] - Cg[3][g]);
297 /* we must do u mod 2147483647 with u in [- 4294966863 ; 4294967066 ] : */
298 if (u < 0) u += 2147483647;
299 if (u < 0) u += 2147483647;
300 if (u >= 2147483647) u -= 2147483647;
301 if (u >= 2147483647) u -= 2147483647;
303 return ((unsigned int) u );