09c9f2aa569d1e8895821e5a2e816c2bc50e14f8
[scilab.git] / scilab / modules / randlib / src / c / mt.c
1 /* A C-program for MT19937: Integer version (1999/10/28)          */
2 /*  genrand() generates one pseudorandom unsigned int (32bit) */
3 /* which is uniformly distributed among 0 to 2^32-1  for each     */
4 /* call. sgenrand(seed) sets initial values to the working area   */
5 /* of 624 words. Before genrand(), sgenrand(seed) must be         */
6 /* called once. (seed is any 32-bit integer.)                     */
7 /*   Coded by Takuji Nishimura, considering the suggestions by    */
8 /* Topher Cooper and Marc Rieffel in July-Aug. 1997.              */
9
10 /* This library is free software under the Artistic license:       */
11 /* see the file COPYING distributed together with this code.       */
12 /* For the verification of the code, its output sequence file      */
13 /* mt19937int.out is attached (2001/4/2)                           */
14
15 /* Copyright (C) 1997, 1999 Makoto Matsumoto and Takuji Nishimura. */
16 /* Any feedback is very welcome. For any question, comments,       */
17 /* see http://www.math.keio.ac.jp/matumoto/emt.html or email       */
18 /* matumoto@math.keio.ac.jp                                        */
19
20 /* REFERENCE                                                       */
21 /* M. Matsumoto and T. Nishimura,                                  */
22 /* "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform  */
23 /* Pseudo-Random Number Generator",                                */
24 /* ACM Transactions on Modeling and Computer Simulation,           */
25 /* Vol. 8, No. 1, January 1998, pp 3--30.                          */
26
27
28 /*   
29  *  NOTES
30  *    slightly modified par Bruno Pincon for inclusion in scilab 
31  *      - names have changed (for uniformity with the others genators)
32  *      - add get state routine
33  *      - add a little verif when the state is changed with the simple
34  *        procedure
35  *
36  *     furthers modifications on May 25 2002 :
37  *
38  *     1/ corrections of the followings : 
39  *        
40  *        bug 1 : the complete state was returned at the scilab level
41  *                without the index mti. Now the complete state is a 
42  *                vector of dim 625 with mti as the first component
43  *        bug 2 : the set_state doesn't work if the generator was not
44  *                initialised => add a is_init var and returned
45  *                the state given with the default initialisation.
46  *
47  *     2/ Following the modif in the new version of this generator I have
48  *       changed the simple initialisation (but not put the init via array)
49  *
50  *     Sept 2005 : fix for bug 1568
51  */
52
53 #include <math.h>
54 #include "grand.h"            /* to check prototypes */
55 #include "sciprint.h"
56 #include "others_generators.h"
57 #include "localization.h"
58
59 /* Period parameters */  
60 #define N 624
61 #define M 397
62 #define MATRIX_A 0x9908b0df   /* constant vector a */
63 #define UPPER_MASK 0x80000000 /* most significant w-r bits */
64 #define LOWER_MASK 0x7fffffff /* least significant r bits */
65
66 /* Tempering parameters */   
67 #define TEMPERING_MASK_B 0x9d2c5680
68 #define TEMPERING_MASK_C 0xefc60000
69 #define TEMPERING_SHIFT_U(y)  (y >> 11)
70 #define TEMPERING_SHIFT_S(y)  (y << 7)
71 #define TEMPERING_SHIFT_T(y)  (y << 15)
72 #define TEMPERING_SHIFT_L(y)  (y >> 18)
73
74 static unsigned long mt[N]; /* the array for the state vector  */
75 static int mti=N;
76 static int is_init=0;  
77 static double DEFAULT_SEED=5489.0;
78
79 unsigned long randmt()
80 {
81     unsigned long y;
82     static unsigned long mag01[2]={0x0, MATRIX_A};
83     /* mag01[x] = x * MATRIX_A  for x=0,1 */
84
85     if (mti >= N) { /* generate N words at one time */
86         int kk;
87
88         if ( ! is_init )
89           set_state_mt_simple(DEFAULT_SEED);
90
91         for (kk=0;kk<N-M;kk++) {
92             y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
93             mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
94         }
95         for (;kk<N-1;kk++) {
96             y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
97             mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
98         }
99         y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
100         mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
101
102         mti = 0;
103     }
104   
105     y = mt[mti++];
106     y ^= TEMPERING_SHIFT_U(y);
107     y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
108     y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
109     y ^= TEMPERING_SHIFT_L(y);
110
111     return ( y );
112 }
113
114
115 int set_state_mt_simple(double s)
116 {
117   /*   set the initial state with the simple procedure  */
118   unsigned long seed;
119
120   if ( s == floor(s) && 0.0 <= s && s <= 4294967295.0)
121     {
122       seed = (unsigned long) s;
123         mt[0]= seed & 0xffffffff;
124       for (mti=1; mti<N; mti++)
125         {
126           mt[mti] =  (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti); 
127           /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
128           /* In the previous versions, MSBs of the seed affect   */
129           /* only MSBs of the array mt[].                        */
130           /* 2002/01/09 modified by Makoto Matsumoto             */
131           mt[mti] &= 0xffffffffUL;   /* for >32 bit machines */
132         }
133       is_init = 1;
134       return ( 1 );
135     }
136   else
137     {
138       sciprint(_("Bad seed for mt, must be an int in [0, 2^32-1]\n"));
139       return ( 0 );
140     }
141 }
142
143
144 /* 
145  *  Initialization by "set_state_simple_mt()" is an example. Theoretically,
146  *  there are 2^19937-1 possible states as an intial state.           
147  *   This function allows to choose any of 2^19937-1 ones.            
148  *   Essential bits in "seed_array[]" is following 19937 bits:        
149  *      (seed_array[0]&UPPER_MASK), seed_array[1], ..., seed_array[N-1].
150  *      (seed_array[0]&LOWER_MASK) is discarded.                         
151  *
152  *   Theoretically,                                                   
153  *      (seed_array[0]&UPPER_MASK), seed_array[1], ..., seed_array[N-1] 
154  *   can take any values except all zeros.                            
155  */
156
157 int set_state_mt(double seed_array[])
158
159 {
160     int i, mti_try;
161
162     mti_try = (int) seed_array[0];
163     if (mti_try < 1  ||  mti_try > 624)
164       {
165         sciprint(_("The first component of the mt state mt, must be an int in [1, 624]\n"));
166         return ( 0 );
167       }
168     is_init = 1;
169     mti = mti_try;
170     for (i=0;i<N;i++) 
171       mt[i] = ((unsigned long) seed_array[i+1]) & 0xffffffff;
172     return ( 1 );
173 }
174
175
176 /*  To return the state at the scilab level  */
177 void get_state_mt(double state[])
178 {
179     int i;
180
181     if ( ! is_init )
182       set_state_mt_simple(DEFAULT_SEED);
183     
184     state[0] = (double) mti;
185     for (i=0;i<N;i++) 
186       state[i+1] = (double) mt[i];
187 }
188
189
190
191
192
193
194
195
196
197