Return to mt19937.c CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / engine

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

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