ThreeB 1.1

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

mt19937.h
Go to the documentation of this file.

    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