PLplot 5.9.6
|
00001 /* 00002 * A C-program for MT19937, with initialization improved 2002/1/26. 00003 * Coded by Takuji Nishimura and Makoto Matsumoto. 00004 * 00005 * Before using, initialize the state by using init_genrand(seed) 00006 * or init_by_array(init_key, key_length). 00007 * 00008 * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, 00009 * All rights reserved. 00010 * Copyright (C) 2005, Mutsuo Saito, 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.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html 00043 * email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) 00044 */ 00045 00046 #include <stdio.h> 00047 #include "mt19937ar.h" 00048 00049 /* Period parameters */ 00050 #define N 624 00051 #define M 397 00052 #define MATRIX_A 0x9908b0dfUL /* constant vector a */ 00053 #define UPPER_MASK 0x80000000UL /* most significant w-r bits */ 00054 #define LOWER_MASK 0x7fffffffUL /* least significant r bits */ 00055 00056 static unsigned long mt[N]; /* the array for the state vector */ 00057 static int mti = N + 1; /* mti==N+1 means mt[N] is not initialized */ 00058 00059 /* initializes mt[N] with a seed */ 00060 void init_genrand( unsigned long s ) 00061 { 00062 mt[0] = s & 0xffffffffUL; 00063 for ( mti = 1; mti < N; mti++ ) 00064 { 00065 mt[mti] = 00066 ( 1812433253UL * ( mt[mti - 1] ^ ( mt[mti - 1] >> 30 ) ) + mti ); 00067 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ 00068 /* In the previous versions, MSBs of the seed affect */ 00069 /* only MSBs of the array mt[]. */ 00070 /* 2002/01/09 modified by Makoto Matsumoto */ 00071 mt[mti] &= 0xffffffffUL; 00072 /* for >32 bit machines */ 00073 } 00074 } 00075 00076 /* initialize by an array with array-length */ 00077 /* init_key is the array for initializing keys */ 00078 /* key_length is its length */ 00079 /* slight change for C++, 2004/2/26 */ 00080 void init_by_array( unsigned long init_key[], int key_length ) 00081 { 00082 int i, j, k; 00083 init_genrand( 19650218UL ); 00084 i = 1; j = 0; 00085 k = ( N > key_length ? N : key_length ); 00086 for (; k; k-- ) 00087 { 00088 mt[i] = ( mt[i] ^ ( ( mt[i - 1] ^ ( mt[i - 1] >> 30 ) ) * 1664525UL ) ) 00089 + init_key[j] + j; /* non linear */ 00090 mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ 00091 i++; j++; 00092 if ( i >= N ) 00093 { 00094 mt[0] = mt[N - 1]; i = 1; 00095 } 00096 if ( j >= key_length ) 00097 j = 0; 00098 } 00099 for ( k = N - 1; k; k-- ) 00100 { 00101 mt[i] = ( mt[i] ^ ( ( mt[i - 1] ^ ( mt[i - 1] >> 30 ) ) * 1566083941UL ) ) 00102 - i; /* non linear */ 00103 mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ 00104 i++; 00105 if ( i >= N ) 00106 { 00107 mt[0] = mt[N - 1]; i = 1; 00108 } 00109 } 00110 00111 mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ 00112 } 00113 00114 /* generates a random number on [0,0xffffffff]-interval */ 00115 unsigned long genrand_int32( void ) 00116 { 00117 unsigned long y; 00118 static unsigned long mag01[2] = { 0x0UL, MATRIX_A }; 00119 /* mag01[x] = x * MATRIX_A for x=0,1 */ 00120 00121 if ( mti >= N ) /* generate N words at one time */ 00122 { 00123 int kk; 00124 00125 if ( mti == N + 1 ) /* if init_genrand() has not been called, */ 00126 init_genrand( 5489UL ); /* a default initial seed is used */ 00127 00128 for ( kk = 0; kk < N - M; kk++ ) 00129 { 00130 y = ( mt[kk] & UPPER_MASK ) | ( mt[kk + 1] & LOWER_MASK ); 00131 mt[kk] = mt[kk + M] ^ ( y >> 1 ) ^ mag01[y & 0x1UL]; 00132 } 00133 for (; kk < N - 1; kk++ ) 00134 { 00135 y = ( mt[kk] & UPPER_MASK ) | ( mt[kk + 1] & LOWER_MASK ); 00136 mt[kk] = mt[kk + ( M - N )] ^ ( y >> 1 ) ^ mag01[y & 0x1UL]; 00137 } 00138 y = ( mt[N - 1] & UPPER_MASK ) | ( mt[0] & LOWER_MASK ); 00139 mt[N - 1] = mt[M - 1] ^ ( y >> 1 ) ^ mag01[y & 0x1UL]; 00140 00141 mti = 0; 00142 } 00143 00144 y = mt[mti++]; 00145 00146 /* Tempering */ 00147 y ^= ( y >> 11 ); 00148 y ^= ( y << 7 ) & 0x9d2c5680UL; 00149 y ^= ( y << 15 ) & 0xefc60000UL; 00150 y ^= ( y >> 18 ); 00151 00152 return y; 00153 } 00154 00155 /* generates a random number on [0,0x7fffffff]-interval */ 00156 long genrand_int31( void ) 00157 { 00158 return (long) ( genrand_int32() >> 1 ); 00159 } 00160 00161 /* generates a random number on [0,1]-real-interval */ 00162 double genrand_real1( void ) 00163 { 00164 return genrand_int32() * ( 1.0 / 4294967295.0 ); 00165 /* divided by 2^32-1 */ 00166 } 00167 00168 /* generates a random number on [0,1)-real-interval */ 00169 double genrand_real2( void ) 00170 { 00171 return genrand_int32() * ( 1.0 / 4294967296.0 ); 00172 /* divided by 2^32 */ 00173 } 00174 00175 /* generates a random number on (0,1)-real-interval */ 00176 double genrand_real3( void ) 00177 { 00178 return ( ( (double) genrand_int32() ) + 0.5 ) * ( 1.0 / 4294967296.0 ); 00179 /* divided by 2^32 */ 00180 } 00181 00182 /* generates a random number on [0,1) with 53-bit resolution*/ 00183 double genrand_res53( void ) 00184 { 00185 unsigned long a = genrand_int32() >> 5, b = genrand_int32() >> 6; 00186 return ( a * 67108864.0 + b ) * ( 1.0 / 9007199254740992.0 ); 00187 } 00188 /* These real versions are due to Isaku Wada, 2002/01/09 added */