Return to mt19937.c CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / engine

Annotation of OpenXM_contrib2/asir2000/engine/mt19937.c, Revision 1.2

1.2     ! noro        1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/mt19937.c,v 1.1.1.1 1999/12/03 07:39:08 noro Exp $
        !             2: */
1.1       noro        3: #include <stdio.h>
                      4:
                      5: /* A C-program for MT19937: Real number version                */
                      6: /*   genrand() generates one pseudorandom real number (double) */
                      7: /* which is uniformly distributed on [0,1]-interval, for each  */
                      8: /* call. sgenrand(seed) set initial values to the working area */
                      9: /* of 624 words. Before genrand(), sgenrand(seed) must be      */
                     10: /* called once. (seed is any 32-bit integer except for 0).     */
                     11: /* Integer generator is obtained by modifying two lines.       */
                     12: /*   Coded by Takuji Nishimura, considering the suggestions by */
                     13: /* Topher Cooper and Marc Rieffel in July-Aug. 1997.           */
                     14:
                     15: /* This library is free software; you can redistribute it and/or   */
                     16: /* modify it under the terms of the GNU Library General Public     */
                     17: /* License as published by the Free Software Foundation; either    */
                     18: /* version 2 of the License, or (at your option) any later         */
                     19: /* version.                                                        */
                     20: /* This library is distributed in the hope that it will be useful, */
                     21: /* but WITHOUT ANY WARRANTY; without even the implied warranty of  */
                     22: /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.            */
                     23: /* See the GNU Library General Public License for more details.    */
                     24: /* You should have received a copy of the GNU Library General      */
                     25: /* Public License along with this library; if not, write to the    */
                     26: /* Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA   */
                     27: /* 02111-1307  USA                                                 */
                     28:
                     29: /* Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.       */
                     30: /* When you use this, send an email to: matumoto@math.keio.ac.jp   */
                     31: /* with an appropriate reference to your work.                     */
                     32:
                     33: #define genrand mt_genrand
                     34: #define sgenrand mt_sgenrand
                     35:
                     36: #include<stdio.h>
                     37:
                     38: /* Period parameters */
                     39: #define N 624
                     40: #define M 397
                     41: #define MATRIX_A 0x9908b0df   /* constant vector a */
                     42: #define UPPER_MASK 0x80000000 /* most significant w-r bits */
                     43: #define LOWER_MASK 0x7fffffff /* least significant r bits */
                     44:
                     45: /* Tempering parameters */
                     46: #define TEMPERING_MASK_B 0x9d2c5680
                     47: #define TEMPERING_MASK_C 0xefc60000
                     48: #define TEMPERING_SHIFT_U(y)  (y >> 11)
                     49: #define TEMPERING_SHIFT_S(y)  (y << 7)
                     50: #define TEMPERING_SHIFT_T(y)  (y << 15)
                     51: #define TEMPERING_SHIFT_L(y)  (y >> 18)
                     52:
                     53: static unsigned long mt[N]; /* the array for the state vector  */
                     54: static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */
                     55:
                     56: int mt_save(name)
                     57: char *name;
                     58: {
                     59:        FILE *fp;
                     60:
                     61:        if ( !(fp = fopen(name,"wb")) )
                     62:                return 0;
                     63:        fwrite(mt,N,sizeof(unsigned long),fp);
                     64:        fwrite(&mti,1,sizeof(int),fp);
                     65:        fclose(fp);
                     66:        return 1;
                     67: }
                     68:
                     69: int mt_load(name)
                     70: char *name;
                     71: {
                     72:        FILE *fp;
                     73:
                     74:        if ( !(fp = fopen(name,"rb")) )
                     75:                return 0;
                     76:        fread(mt,N,sizeof(unsigned long),fp);
                     77:        fread(&mti,1,sizeof(int),fp);
                     78:        fclose(fp);
                     79:        return 1;
                     80: }
                     81:
                     82: /* initializing the array with a NONZERO seed */
                     83: void
                     84: sgenrand(seed)
                     85:     unsigned long seed;
                     86: {
                     87:     /* setting initial seeds to mt[N] using         */
                     88:     /* the generator Line 25 of Table 1 in          */
                     89:     /* [KNUTH 1981, The Art of Computer Programming */
                     90:     /*    Vol. 2 (2nd Ed.), pp102]                  */
                     91:     mt[0]= seed & 0xffffffff;
                     92:     for (mti=1; mti<N; mti++)
                     93:         mt[mti] = (69069 * mt[mti-1]) & 0xffffffff;
                     94: }
                     95:
                     96: #if 0
                     97: double /* generating reals */
                     98: #endif
                     99: unsigned long /* for integer generation */
                    100: genrand()
                    101: {
                    102:     unsigned long y;
                    103:     static unsigned long mag01[2]={0x0, MATRIX_A};
                    104:     /* mag01[x] = x * MATRIX_A  for x=0,1 */
                    105:
                    106:     if (mti >= N) { /* generate N words at one time */
                    107:         int kk;
                    108:
                    109:         if (mti == N+1)   /* if sgenrand() has not been called, */
                    110:             sgenrand(4357); /* a default initial seed is used   */
                    111:
                    112:         for (kk=0;kk<N-M;kk++) {
                    113:             y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
                    114:             mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
                    115:         }
                    116:         for (;kk<N-1;kk++) {
                    117:             y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
                    118:             mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
                    119:         }
                    120:         y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
                    121:         mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
                    122:
                    123:         mti = 0;
                    124:     }
                    125:
                    126:     y = mt[mti++];
                    127:     y ^= TEMPERING_SHIFT_U(y);
                    128:     y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
                    129:     y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
                    130:     y ^= TEMPERING_SHIFT_L(y);
                    131:
                    132: #if 0
                    133:     return ( (double)y / (unsigned long)0xffffffff ); /* reals */
                    134: #endif
                    135:      return y; /* for integer generation */
                    136: }
                    137:
                    138: #if 0
                    139: /* this main() outputs first 1000 generated numbers  */
                    140: main()
                    141: {
                    142:     int j;
                    143:
                    144:     sgenrand(4357); /* any nonzero integer can be used as a seed */
                    145:     for (j=0; j<1000; j++) {
                    146:         printf("%5f ", genrand());
                    147:         if (j%8==7) printf("\n");
                    148:     }
                    149:     printf("\n");
                    150: }
                    151: #endif