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 #ifndef MT19937_H 00021 #define MT19937_H 00022 00023 #include "randomc.h" 00024 #include <iostream> 00025 #include <sstream> 00026 #include <string> 00027 #include <cmath> 00028 #include <iomanip> 00029 00030 ///Useful wrapper for MT19937 random number generator class. 00031 ///@ingroup gUtility 00032 struct MT19937 00033 { 00034 ///Null struct thrown if attempting to load 00035 ///state from stream yields a parse error. 00036 struct ParseError{}; 00037 00038 ///Underlying RNG. 00039 CRandomMersenne rng; 00040 00041 public: 00042 MT19937() 00043 :rng(0) 00044 {} 00045 00046 ///Seed state with a simple RNG 00047 ///@param seed 00048 void simple_seed(int seed) 00049 { 00050 rng.RandomInit(seed); 00051 } 00052 00053 ///Duplicate RNG state 00054 ///@param r RNG to duplicate 00055 void copy_state(const MT19937& r) 00056 { 00057 rng = r.rng; 00058 } 00059 00060 ///Generate a double 00061 double operator()() 00062 { 00063 return rng.Random(); 00064 } 00065 00066 ///Generate an int 00067 uint32_t rand_int() 00068 { 00069 return rng.BRandom(); 00070 } 00071 00072 ///Generate a Gaussian variate 00073 double gaussian() 00074 { 00075 double x1, x2, w; 00076 do { 00077 x1 = 2.0 * (*this)() - 1.0; 00078 x2 = 2.0 * (*this)() - 1.0; 00079 w = x1 * x1 + x2 * x2; 00080 } while ( w >= 1.0 ); 00081 00082 w = std::sqrt( (-2.0 * std::log( w ) ) / w ); 00083 return x1 * w; 00084 //spare so we don't have to save that one extra bit of state y2 = x2 * w; 00085 } 00086 00087 ///Serialise state 00088 ///@param o Stream to serialise to 00089 void write(std::ostream& o) 00090 { 00091 using namespace std; 00092 char f = o.fill(); 00093 ios_base::fmtflags fl = o.flags(); 00094 o << "MT19937 " << hex << setfill('0') << setw(3) << rng.get_index(); 00095 for(int i=0; i < MERS_N; i++) 00096 o << " " << hex << setw(8) << rng.get_state()[i]; 00097 00098 o << setfill(f) << setiosflags(fl); 00099 } 00100 00101 ///De serialise state 00102 ///param is Stream to de-serialise from 00103 void read(std::istream& is) 00104 { 00105 using namespace std; 00106 00107 string ls; 00108 getline(is, ls); 00109 if(ls.size() != 5627) 00110 { 00111 cerr << "MT19937: Expected string of length 5627. Got " << ls.size() << endl; 00112 throw ParseError(); 00113 } 00114 00115 istringstream l(ls); 00116 00117 string s; 00118 uint32_t i; 00119 00120 l >> s; 00121 00122 if(s != "MT19937") 00123 { 00124 cerr << "MT19937: Expected MT19937. Got " << s << endl; 00125 throw ParseError(); 00126 } 00127 00128 for(int n=0; n < MERS_N + 1; n++) 00129 { 00130 l >> hex >> i; 00131 if(l.bad()) 00132 { 00133 cerr << "MT19937: Expected hex number. Got "; 00134 if(l.eof()) 00135 cerr << "EOF" << endl; 00136 else 00137 { 00138 cerr << l.get() << endl; 00139 } 00140 00141 throw ParseError(); 00142 } 00143 00144 if(n==0) 00145 rng.get_index() = i; 00146 else 00147 rng.get_state()[n-1]=i; 00148 00149 } 00150 00151 } 00152 private: 00153 /// Disallow copying, since one almost never wants to do this. 00154 /// Copying has to be explicit via MT19937::copy_state(). 00155 MT19937(const MT19937&); 00156 00157 }; 00158 00159 00160 #endif