Annotation of OpenXM_contrib2/asir2000/engine/mt19937.c, Revision 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
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to mt19937.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine