Annotation of OpenXM_contrib2/asir2000/engine/mt19937.c, Revision 1.3
1.3 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/mt19937.c,v 1.2 2000/08/21 08:31:28 noro Exp $
1.2 noro 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: {
1.3 ! noro 59: FILE *fp;
1.1 noro 60:
1.3 ! noro 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;
1.1 noro 67: }
68:
69: int mt_load(name)
70: char *name;
71: {
1.3 ! noro 72: FILE *fp;
1.1 noro 73:
1.3 ! noro 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;
1.1 noro 80: }
81:
82: /* initializing the array with a NONZERO seed */
83: void
84: sgenrand(seed)
1.3 ! noro 85: unsigned long seed;
1.1 noro 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
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