* Bug #8597 fixed - Uncontrolled message of grand/clcg4 should be displayed as warning.
[scilab.git] / scilab / modules / randlib / src / c / clcg4.c
1 /*
2  * 
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
7  * 
8  *  PURPOSE
9  *     clcg4 generator stuff
10  *
11  *  AUTHORS
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).
15  *
16  *  NOTES
17  *     The original code was slightly modified by Bruno Pincon for inclusion
18  *     in Scilab. 
19  *
20  *     list of main modifs :
21  *
22  *       - lot of routine 's names have changed to have some kind of
23  *         uniformity with the others generators 
24  *
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 ...
29  *
30  *       - add a routine advance_state_clcg4 (for compatibility with the
31  *         old package (Scilab used this feature))
32  *
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 &
35  *         Andres) :
36  *
37  *         From the 4 LCG :
38  *
39  *            x_{j,n} = a_j * x_{j,n-1} mod m_j    0 <= j <= 3
40  *
41  *         The output with form (2) (original form in this code) :
42  *         
43  *            z_n = ( sum_j  delta_j * x_{j,n} / m_j ) mod 1
44  *
45  *         have been changed in the form (1) :
46  *
47  *           z_n = ( sum_j  delta_j * x_{j,n} ) mod m_1 (then u_n = z_n / m_1)
48  *         
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.
52  */
53
54
55 /*---------------------------------------------------------------------*/
56 /* clcg4.c   Implementation module                                     */
57 /*---------------------------------------------------------------------*/
58
59 #include "clcg4.h"
60 #include <math.h>             /* for floor */
61 #include "sciprint.h"
62 #include "others_generators.h"
63 #include "localization.h"
64 #include "warningmode.h"
65
66 /***********************************************************************/
67 /* Private part.                                                       */
68 /***********************************************************************/
69
70 #define H   32768               /* = 2^15 : use in MultModM.           */
71
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 };
75
76 static int Ig[4][Maxgen+1], Lg[4][Maxgen+1], Cg[4][Maxgen+1];
77 /* Initial seed, previous seed, and current seed. */
78
79
80 static int  is_init = 0;
81 static int v_default = 31;
82 static int w_default = 41;
83
84
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).                                    */
88 {
89     int R, S0, S1, q, qh, rh, k;
90
91     if (s < 0)  s += M;
92     if (t < 0)  t += M;
93     if (s < H)  { S0 = s;  R = 0; }
94     else
95     {
96         S1 = s/H;  S0 = s - H*S1;
97         qh = M/H;  rh = M - H*qh;
98         if (S1 >= H)
99         {
100             S1 -= H;   k = t/qh;   R = H * (t - k*qh) - k*rh;
101             while (R < 0)  R += M;
102         }
103         else R = 0;
104         if (S1 != 0)
105         {
106             q = M/S1;   k = t/q;   R -= k * (M - S1*q);
107             if (R > 0)  R -= M;
108             R += S1*(t - k*q);
109             while (R < 0)  R += M;
110         }
111         k = R/qh;   R = H * (R - k*qh) - k*rh;
112         while (R < 0) R += M;
113     }
114     if (S0 != 0)
115     {
116         q = M/S0;   k = t/q;   R -= k* (M - S0*q);
117         if (R > 0)  R -= M;
118         R += S0 * (t - k*q);
119         while (R < 0)  R += M;
120     }
121     return R;
122 }
123
124 static void comp_aw_and_avw(int v, int w)
125 {
126     int i, j;
127     for (j = 0; j < 4; j++)
128     {
129         aw [j] = a [j];
130         for (i = 1; i <= w; i++)
131             aw [j]  = MultModM (aw [j], aw [j], m[j]);
132         avw [j] = aw [j];
133         for (i = 1; i <= v; i++)
134             avw [j] = MultModM (avw [j], avw [j], m[j]);
135     }
136 }
137
138 static void init_clcg4(int v, int w)
139 {
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
144     * is_init)
145     */
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]);
149 }
150
151 static int verif_seeds_clcg4(double s0, double s1, double s2, double s3)
152 {
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 )
160         return ( 1 );
161     else
162         return ( 0 );
163 }
164
165 static void display_info_clcg4(void)
166 {
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]"));
169 }
170
171
172 /*---------------------------------------------------------------------*/
173 /* Public part.                                                        */
174 /*---------------------------------------------------------------------*/
175
176
177 int set_seed_clcg4(int g, double s0, double s1, double s2, double s3)
178 {
179     if (! is_init ) {init_clcg4(v_default,w_default); is_init = 1; };
180
181     if ( verif_seeds_clcg4(s0, s1, s2, s3) )
182     {
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())
187         {
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);
189         }
190         return ( 1 );
191     }
192     else
193     {
194         display_info_clcg4();
195         return ( 0 );
196     }
197 }
198
199 void get_state_clcg4(int g, double s[4])
200 {
201     int j;
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];
204 }
205
206 void init_generator_clcg4(int g, SeedType Where)
207 {
208     int j;
209     if (! is_init ) {init_clcg4(v_default,w_default); is_init = 1; };
210     for (j = 0; j < 4; j++)
211     {
212         switch (Where)
213         {
214         case InitialSeed :
215             Lg [j][g] = Ig [j][g];   break;
216         case NewSeed :
217             Lg [j][g] = MultModM (aw [j], Lg [j][g], m [j]);   break;
218         case LastSeed :
219             break;
220         }
221         Cg [j][g] = Lg [j][g];
222     }
223 }
224
225 void advance_state_clcg4(int g, int k)
226 {
227     int b[4];
228     int i, j;
229
230     if (! is_init ) {init_clcg4(v_default,w_default); is_init = 1; };
231
232     for ( j = 0 ; j < 4 ; j++ )
233     {
234         b[j] = a[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] );
238     }
239     init_generator_clcg4(g, InitialSeed);
240 }
241
242 int set_initial_seed_clcg4(double s0, double s1, double s2, double s3)
243 {
244     int g, j;
245
246     if (! is_init )  comp_aw_and_avw(v_default,w_default);
247
248     if ( ! verif_seeds_clcg4(s0, s1, s2, s3) )
249     {
250         display_info_clcg4();
251         return ( 0 );
252     };
253
254     is_init = 1;
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++)
261     {
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);
265     }
266     return ( 1 );
267 }
268
269 unsigned int clcg4(int g)
270 {
271     /* Modif Bruno : the generator have now the form (1) in place of (2) */
272
273     int k,s;
274     double u;
275
276     if (! is_init ) {init_clcg4(v_default,w_default); is_init = 1; };
277
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;
282
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;
286
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;
290
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;
294
295     /*  final step */
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;
302
303     return ((unsigned int) u );
304
305 }