cgma
RandomMersenne.cpp
Go to the documentation of this file.
00001 /* 
00002    A C-program for MT19937, with initialization improved 2002/2/10.
00003    Coded by Takuji Nishimura and Makoto Matsumoto.
00004    This is a faster version by taking Shawn Cokus's optimization,
00005    Matthe Bellew's simplification, Isaku Wada's real version.
00006 
00007    Before using, initialize the state by using init_genrand(seed) 
00008    or init_by_array(init_key, key_length).
00009 
00010    Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
00011    All rights reserved.                          
00012 
00013    Redistribution and use in source and binary forms, with or without
00014    modification, are permitted provided that the following conditions
00015    are met:
00016 
00017      1. Redistributions of source code must retain the above copyright
00018         notice, this list of conditions and the following disclaimer.
00019 
00020      2. Redistributions in binary form must reproduce the above copyright
00021         notice, this list of conditions and the following disclaimer in the
00022         documentation and/or other materials provided with the distribution.
00023 
00024      3. The names of its contributors may not be used to endorse or promote 
00025         products derived from this software without specific prior written 
00026         permission.
00027 
00028    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
00029    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
00030    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
00031    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
00032    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00033    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00034    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00035    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00036    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00037    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00038    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00039 
00040 
00041    Any feedback is very welcome.
00042    http://www.math.keio.ac.jp/matumoto/emt.html
00043    email: matumoto@math.keio.ac.jp
00044 */
00045 
00046 #include "RandomMersenne.hpp"
00047 
00048 /* Period parameters */  
00049 #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
00050 #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL))
00051 
00052 /* initializes state[N] with a seed */
00053 void RandomMersenne::init_genrand(unsigned long s)
00054 {
00055   mt[0]= s & 0xffffffffUL;
00056   for (mti=1; mti<N; mti++) {
00057     mt[mti] = 
00058       (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti); 
00059       /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
00060       /* In the previous versions, MSBs of the seed affect   */
00061       /* only MSBs of the array mt[].                        */
00062       /* 2002/01/09 modified by Makoto Matsumoto             */
00063     mt[mti] &= 0xffffffffUL;
00064       /* for >32 bit machines */
00065   }
00066 }
00067 
00068 /* initialize by an array with array-length */
00069 /* init_key is the array for initializing keys */
00070 /* key_length is its length */
00071 void RandomMersenne::init_by_array(unsigned long init_key[], unsigned long key_length)
00072 {
00073   unsigned long i, j, k;
00074   const unsigned long n = (unsigned long)N;
00075   init_genrand(19650218UL);
00076   i=1UL; j=0UL;
00077   k = (n>key_length ? n : key_length);
00078   for (; k; k--) {
00079     mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
00080       + init_key[j] + j; /* non linear */
00081     mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
00082     i++; j++;
00083     if (i>=n) { mt[0] = mt[n-1]; i=1; }
00084     if (j>=key_length) j=0;
00085   }
00086   for (k=n-1; k; k--) {
00087     mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
00088       - i; /* non linear */
00089     mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
00090     i++;
00091     if (i>=n) { mt[0] = mt[n-1]; i=1; }
00092   }
00093 
00094   mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ 
00095 }
00096 
00097 unsigned long RandomMersenne::genrand_int32() 
00098 {
00099   unsigned long y;
00100   static unsigned long mag01[2]={0x0UL, MATRIX_A};
00101     /* mag01[x] = x * MATRIX_A  for x=0,1 */
00102 
00103   if (mti >= N) { /* generate N words at one time */
00104     unsigned kk;
00105 
00106       //if (mti == N+1)   /* if init_genrand() has not been called, */
00107       //init_genrand(5489UL); /* a default initial seed is used */
00108 
00109     for (kk=0;kk<N-M;kk++) {
00110       y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
00111       mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
00112     }
00113     for (;kk<N-1;kk++) {
00114       y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
00115       mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
00116     }
00117     y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
00118     mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
00119 
00120     mti = 0;
00121   }
00122   
00123   y = mt[mti++];
00124 
00125     /* Tempering */
00126   y ^= (y >> 11);
00127   y ^= (y << 7) & 0x9d2c5680UL;
00128   y ^= (y << 15) & 0xefc60000UL;
00129   y ^= (y >> 18);
00130 
00131   return y;
00132 }
00133 
00134 #ifdef TEST_RANDOM
00135 
00136 #include <cstdio>
00137 
00138 int main(void)
00139 {
00140     int i;
00141     unsigned long init[4]={0x123, 0x234, 0x345, 0x456}, length=4;
00142     RandomMersenne rm(init, length);
00143     printf("1000 outputs of genrand_int32()\n");
00144     for (i=0; i<1000; i++) {
00145       printf("%10lu ", rm.genrand_int32());
00146       if (i%5==4) printf("\n");
00147     }
00148     printf("\n1000 outputs of genrand_real2()\n");
00149     for (i=0; i<1000; i++) {
00150       printf("%10.8f ", rm.genrand_real2());
00151       if (i%5==4) printf("\n");
00152     }
00153     return 0;
00154 }
00155 #endif
00156 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines