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 }