ThreeB 1.1
|
00001 /* 00002 This file is part of B-cubed. 00003 00004 Copyright (C) 2009, 2010, 2011, Edward Rosten and Susan Cox 00005 00006 B-cubed is free software; you can redistribute it and/or 00007 modify it under the terms of the GNU Lesser General Public 00008 License as published by the Free Software Foundation; either 00009 version 3.0 of the License, or (at your option) any later version. 00010 00011 This library is distributed in the hope that it will be useful, 00012 but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00014 Lesser General Public License for more details. 00015 00016 You should have received a copy of the GNU General Public License 00017 along with this program. If not, see <http://www.gnu.org/licenses/> 00018 */ 00019 00020 /************************** mersenne.cpp ********************************** 00021 * Author: Agner Fog 00022 * Date created: 2001 00023 * Last modified: 2008-11-16 00024 * Project: randomc.h 00025 * Platform: Any C++ 00026 * Description: 00027 * Random Number generator of type 'Mersenne Twister' 00028 * 00029 * This random number generator is described in the article by 00030 * M. Matsumoto & T. Nishimura, in: 00031 * ACM Transactions on Modeling and Computer Simulation, 00032 * vol. 8, no. 1, 1998, pp. 3-30. 00033 * Details on the initialization scheme can be found at 00034 * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html 00035 * 00036 * Further documentation: 00037 * The file ran-instructions.pdf contains further documentation and 00038 * instructions. 00039 * 00040 * Copyright 2001-2008 by Agner Fog. 00041 * GNU General Public License http://www.gnu.org/licenses/gpl.html 00042 *******************************************************************************/ 00043 00044 #include "randomc.h" 00045 00046 void CRandomMersenne::Init0(int seed) { 00047 // Seed generator 00048 const uint32_t factor = 1812433253UL; 00049 mt[0]= seed; 00050 for (mti=1; mti < MERS_N; mti++) { 00051 mt[mti] = (factor * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti); 00052 } 00053 } 00054 00055 void CRandomMersenne::RandomInit(int seed) { 00056 // Initialize and seed 00057 Init0(seed); 00058 00059 // Randomize some more 00060 for (int i = 0; i < 37; i++) BRandom(); 00061 } 00062 00063 00064 void CRandomMersenne::RandomInitByArray(int const seeds[], int NumSeeds) { 00065 // Seed by more than 32 bits 00066 int i, j, k; 00067 00068 // Initialize 00069 Init0(19650218); 00070 00071 if (NumSeeds <= 0) return; 00072 00073 // Randomize mt[] using whole seeds[] array 00074 i = 1; j = 0; 00075 k = (MERS_N > NumSeeds ? MERS_N : NumSeeds); 00076 for (; k; k--) { 00077 mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) + (uint32_t)seeds[j] + j; 00078 i++; j++; 00079 if (i >= MERS_N) {mt[0] = mt[MERS_N-1]; i=1;} 00080 if (j >= NumSeeds) j=0;} 00081 for (k = MERS_N-1; k; k--) { 00082 mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) - i; 00083 if (++i >= MERS_N) {mt[0] = mt[MERS_N-1]; i=1;}} 00084 mt[0] = 0x80000000UL; // MSB is 1; assuring non-zero initial array 00085 00086 // Randomize some more 00087 mti = 0; 00088 for (int i = 0; i <= MERS_N; i++) BRandom(); 00089 } 00090 00091 00092 uint32_t CRandomMersenne::BRandom() { 00093 // Generate 32 random bits 00094 uint32_t y; 00095 00096 if (mti >= MERS_N) { 00097 // Generate MERS_N words at one time 00098 const uint32_t LOWER_MASK = (1LU << MERS_R) - 1; // Lower MERS_R bits 00099 const uint32_t UPPER_MASK = 0xFFFFFFFF << MERS_R; // Upper (32 - MERS_R) bits 00100 static const uint32_t mag01[2] = {0, MERS_A}; 00101 00102 int kk; 00103 for (kk=0; kk < MERS_N-MERS_M; kk++) { 00104 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK); 00105 mt[kk] = mt[kk+MERS_M] ^ (y >> 1) ^ mag01[y & 1];} 00106 00107 for (; kk < MERS_N-1; kk++) { 00108 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK); 00109 mt[kk] = mt[kk+(MERS_M-MERS_N)] ^ (y >> 1) ^ mag01[y & 1];} 00110 00111 y = (mt[MERS_N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK); 00112 mt[MERS_N-1] = mt[MERS_M-1] ^ (y >> 1) ^ mag01[y & 1]; 00113 mti = 0; 00114 } 00115 y = mt[mti++]; 00116 00117 // Tempering (May be omitted): 00118 y ^= y >> MERS_U; 00119 y ^= (y << MERS_S) & MERS_B; 00120 y ^= (y << MERS_T) & MERS_C; 00121 y ^= y >> MERS_L; 00122 00123 return y; 00124 } 00125 00126 00127 double CRandomMersenne::Random() { 00128 // Output random float number in the interval 0 <= x < 1 00129 // Multiply by 2^(-32) 00130 return (double)BRandom() * (1./(65536.*65536.)); 00131 } 00132 00133 00134 int CRandomMersenne::IRandom(int min, int max) { 00135 // Output random integer in the interval min <= x <= max 00136 // Relative error on frequencies < 2^-32 00137 if (max <= min) { 00138 if (max == min) return min; else return 0x80000000; 00139 } 00140 // Multiply interval with random and truncate 00141 int r = int((double)(uint32_t)(max - min + 1) * Random() + min); 00142 if (r > max) r = max; 00143 return r; 00144 } 00145 00146 00147 int CRandomMersenne::IRandomX(int min, int max) { 00148 // Output random integer in the interval min <= x <= max 00149 // Each output value has exactly the same probability. 00150 // This is obtained by rejecting certain bit values so that the number 00151 // of possible bit values is divisible by the interval length 00152 if (max <= min) { 00153 if (max == min) return min; else return 0x80000000; 00154 } 00155 #ifdef INT64_SUPPORTED 00156 // 64 bit integers available. Use multiply and shift method 00157 uint32_t interval; // Length of interval 00158 uint64_t longran; // Random bits * interval 00159 uint32_t iran; // Longran / 2^32 00160 uint32_t remainder; // Longran % 2^32 00161 00162 interval = uint32_t(max - min + 1); 00163 if (interval != LastInterval) { 00164 // Interval length has changed. Must calculate rejection limit 00165 // Reject when remainder >= 2^32 / interval * interval 00166 // RLimit will be 0 if interval is a power of 2. No rejection then 00167 RLimit = uint32_t(((uint64_t)1 << 32) / interval) * interval - 1; 00168 LastInterval = interval; 00169 } 00170 do { // Rejection loop 00171 longran = (uint64_t)BRandom() * interval; 00172 iran = (uint32_t)(longran >> 32); 00173 remainder = (uint32_t)longran; 00174 } while (remainder > RLimit); 00175 // Convert back to signed and return result 00176 return (int32_t)iran + min; 00177 00178 #else 00179 // 64 bit integers not available. Use modulo method 00180 uint32_t interval; // Length of interval 00181 uint32_t bran; // Random bits 00182 uint32_t iran; // bran / interval 00183 uint32_t remainder; // bran % interval 00184 00185 interval = uint32_t(max - min + 1); 00186 if (interval != LastInterval) { 00187 // Interval length has changed. Must calculate rejection limit 00188 // Reject when iran = 2^32 / interval 00189 // We can't make 2^32 so we use 2^32-1 and correct afterwards 00190 RLimit = (uint32_t)0xFFFFFFFF / interval; 00191 if ((uint32_t)0xFFFFFFFF % interval == interval - 1) RLimit++; 00192 } 00193 do { // Rejection loop 00194 bran = BRandom(); 00195 iran = bran / interval; 00196 remainder = bran % interval; 00197 } while (iran >= RLimit); 00198 // Convert back to signed and return result 00199 return (int32_t)remainder + min; 00200 00201 #endif 00202 }