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 }