Revert "* Bug #8597 fixed - Uncontrolled message of grand/clcg4 should be displayed...
[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
65 /***********************************************************************/
66 /* Private part.                                                       */
67 /***********************************************************************/
68
69 #define H   32768               /* = 2^15 : use in MultModM.           */
70
71 static int aw[4], avw[4],      /*   a[j]^{2^w} et a[j]^{2^{v+w}}.     */
72 a[4] = { 45991, 207707, 138556, 49689 },
73 m[4] = { 2147483647, 2147483543, 2147483423, 2147483323 };
74
75 static int Ig[4][Maxgen+1], Lg[4][Maxgen+1], Cg[4][Maxgen+1];
76 /* Initial seed, previous seed, and current seed. */
77
78
79 static int  is_init = 0;
80 static int v_default = 31;
81 static int w_default = 41;
82
83
84 static int MultModM (int s, int t, int M)
85 /* Returns (s*t) MOD M.  Assumes that -M < s < M and -M < t < M.    */
86 /* See L'Ecuyer and Cote (1991).                                    */
87 {
88     int R, S0, S1, q, qh, rh, k;
89
90     if (s < 0)  s += M;
91     if (t < 0)  t += M;
92     if (s < H)  { S0 = s;  R = 0; }
93     else
94     {
95         S1 = s/H;  S0 = s - H*S1;
96         qh = M/H;  rh = M - H*qh;
97         if (S1 >= H)
98         {
99             S1 -= H;   k = t/qh;   R = H * (t - k*qh) - k*rh;
100             while (R < 0)  R += M;
101         }
102         else R = 0;
103         if (S1 != 0)
104         {
105             q = M/S1;   k = t/q;   R -= k * (M - S1*q);
106             if (R > 0)  R -= M;
107             R += S1*(t - k*q);
108             while (R < 0)  R += M;
109         }
110         k = R/qh;   R = H * (R - k*qh) - k*rh;
111         while (R < 0) R += M;
112     }
113     if (S0 != 0)
114     {
115         q = M/S0;   k = t/q;   R -= k* (M - S0*q);
116         if (R > 0)  R -= M;
117         R += S0 * (t - k*q);
118         while (R < 0)  R += M;
119     }
120     return R;
121 }
122
123 static void comp_aw_and_avw(int v, int w)
124 {
125     int i, j;
126     for (j = 0; j < 4; j++)
127     {
128         aw [j] = a [j];
129         for (i = 1; i <= w; i++)
130             aw [j]  = MultModM (aw [j], aw [j], m[j]);
131         avw [j] = aw [j];
132         for (i = 1; i <= v; i++)
133             avw [j] = MultModM (avw [j], avw [j], m[j]);
134     }
135 }
136
137 static void init_clcg4(int v, int w)
138 {
139     /* currently the scilab interface don't let the user chooses
140     * v and w (always v_default and w_default) so this routine
141     * is in the "private" part (also because initialisation is
142     * always perform inside this module, depending of the var
143     * is_init)
144     */
145     double sd[4] = {11111111., 22222222., 33333333., 44444444.};
146     comp_aw_and_avw(v, w);
147     set_initial_seed_clcg4(sd[0], sd[1], sd[2], sd[3]);
148 }
149
150 static int verif_seeds_clcg4(double s0, double s1, double s2, double s3)
151 {
152     /* verify that the seeds are "integers" and are in the good range */
153     if ( s0 == floor(s0) && s1 == floor(s1) &&
154         s2 == floor(s2) && s3 == floor(s3) &&
155         1 <= s0  &&  s0 <= 2147483646      &&
156         1 <= s1  &&  s1 <= 2147483542      &&
157         1 <= s2  &&  s2 <= 2147483422      &&
158         1 <= s3  &&  s3 <= 2147483322 )
159         return ( 1 );
160     else
161         return ( 0 );
162 }
163
164 static void display_info_clcg4(void)
165 {
166     /* display the seeds range (in case of error) */
167     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]"));
168 }
169
170
171 /*---------------------------------------------------------------------*/
172 /* Public part.                                                        */
173 /*---------------------------------------------------------------------*/
174
175
176 int set_seed_clcg4(int g, double s0, double s1, double s2, double s3)
177 {
178     if (! is_init ) {init_clcg4(v_default,w_default); is_init = 1; };
179
180     if ( verif_seeds_clcg4(s0, s1, s2, s3) )
181     {
182         Ig [0][g] = (int) s0; Ig [1][g] = (int) s1;
183         Ig [2][g] = (int) s2; Ig [3][g] = (int) s3;
184         init_generator_clcg4(g, InitialSeed);
185         sciprint(_("\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);
186         return ( 1 );
187     }
188     else
189     {
190         display_info_clcg4();
191         return ( 0 );
192     }
193 }
194
195 void get_state_clcg4(int g, double s[4])
196 {
197     int j;
198     if (! is_init ) {init_clcg4(v_default,w_default); is_init = 1; };
199     for (j = 0; j < 4; j++)  s [j] = (double) Cg [j][g];
200 }
201
202 void init_generator_clcg4(int g, SeedType Where)
203 {
204     int j;
205     if (! is_init ) {init_clcg4(v_default,w_default); is_init = 1; };
206     for (j = 0; j < 4; j++)
207     {
208         switch (Where)
209         {
210         case InitialSeed :
211             Lg [j][g] = Ig [j][g];   break;
212         case NewSeed :
213             Lg [j][g] = MultModM (aw [j], Lg [j][g], m [j]);   break;
214         case LastSeed :
215             break;
216         }
217         Cg [j][g] = Lg [j][g];
218     }
219 }
220
221 void advance_state_clcg4(int g, int k)
222 {
223     int b[4];
224     int i, j;
225
226     if (! is_init ) {init_clcg4(v_default,w_default); is_init = 1; };
227
228     for ( j = 0 ; j < 4 ; j++ )
229     {
230         b[j] = a[j];
231         for ( i = 1 ; i <= k ; i++ )
232             b[j] = MultModM( b[j], b[j], m[j]);
233         Ig[j][g] = MultModM ( b[j], Cg[j][g], m[j] );
234     }
235     init_generator_clcg4(g, InitialSeed);
236 }
237
238 int set_initial_seed_clcg4(double s0, double s1, double s2, double s3)
239 {
240     int g, j;
241
242     if (! is_init )  comp_aw_and_avw(v_default,w_default);
243
244     if ( ! verif_seeds_clcg4(s0, s1, s2, s3) )
245     {
246         display_info_clcg4();
247         return ( 0 );
248     };
249
250     is_init = 1;
251     Ig [0][0] = (int) s0;
252     Ig [1][0] = (int) s1;
253     Ig [2][0] = (int) s2;
254     Ig [3][0] = (int) s3;
255     init_generator_clcg4(0, InitialSeed);
256     for (g = 1; g <= Maxgen; g++)
257     {
258         for (j = 0; j < 4; j++)
259             Ig [j][g] = MultModM (avw [j], Ig [j][g-1], m [j]);
260         init_generator_clcg4(g, InitialSeed);
261     }
262     return ( 1 );
263 }
264
265 unsigned int clcg4(int g)
266 {
267     /* Modif Bruno : the generator have now the form (1) in place of (2) */
268
269     int k,s;
270     double u;
271
272     if (! is_init ) {init_clcg4(v_default,w_default); is_init = 1; };
273
274     /*  advance the 4 LCG */
275     s = Cg [0][g];  k = s / 46693;
276     s = 45991 * (s - k * 46693) - k * 25884;
277     if (s < 0) s = s + 2147483647;  Cg [0][g] = s;
278
279     s = Cg [1][g];  k = s / 10339;
280     s = 207707 * (s - k * 10339) - k * 870;
281     if (s < 0) s = s + 2147483543;  Cg [1][g] = s;
282
283     s = Cg [2][g];  k = s / 15499;
284     s = 138556 * (s - k * 15499) - k * 3979;
285     if (s < 0) s = s + 2147483423;  Cg [2][g] = s;
286
287     s = Cg [3][g];  k = s / 43218;
288     s = 49689 * (s - k * 43218) - k * 24121;
289     if (s < 0) s = s + 2147483323;  Cg [3][g] = s;
290
291     /*  final step */
292     u = (double)(Cg[0][g] - Cg[1][g]) + (double)(Cg[2][g] - Cg[3][g]);
293     /*  we must do  u mod 2147483647 with u in [- 4294966863 ; 4294967066 ] : */
294     if (u < 0) u += 2147483647;
295     if (u < 0) u += 2147483647;
296     if (u >= 2147483647) u -= 2147483647;
297     if (u >= 2147483647) u -= 2147483647;
298
299     return ((unsigned int) u );
300
301 }