mtrand.h

#ifndef UTIL_MTRAND_H
#define UTIL_MTRAND_H

// mtrand.h
// C++ include file for MT19937, with initialization improved 2002/1/26.
// Coded by Takuji Nishimura and Makoto Matsumoto.
// Ported to C++ by Jasper Bedaux 2003/1/1 (see http://www.bedaux.net/mtrand/).
// Modified by David Morse 2011/12/1
// The generators returning floating point numbers are based on
// a version by Isaku Wada, 2002/01/09
//
// Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions
// are met:
//
// 1. Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//
// 3. The names of its contributors may not be used to endorse or promote
//    products derived from this software without specific prior written
//    permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Any feedback is very welcome.
// http://www.math.keio.ac.jp/matumoto/emt.html
// email: matumoto@math.keio.ac.jp
//
// Feedback about the C++ port should be sent to Jasper Bedaux,
// see http://www.bedaux.net/mtrand/ for e-mail address and info.
//
// Modifications by David Morse:
//
// 1) Made state variables non-static members of MTRand_int32,
//    rather than static members.
//
// 2) Added a constructor to construct MTRand from an MTRand_int32,
//    to allow use of original test, in which two generators shared
//    state through shared static variables.
//
// 3) Removed private declarations of copy constructors and assigment,
//    allowing use of the default versions.

namespace Util 
{

   class MTRand_int32
   {
   public:

      MTRand_int32();

      MTRand_int32(unsigned long s);

      MTRand_int32(const unsigned long* array, int size);

      MTRand_int32(const MTRand_int32&);

      void operator=(const MTRand_int32&);

      virtual ~MTRand_int32() {} // destructor

      void seed(unsigned long);

      void seed(const unsigned long*, int size);

      unsigned long operator()();

      template <class Archive>
      void serialize(Archive& ar, const unsigned int version);

   protected:

      unsigned long rand_int32(); // generate 32 bit random integer

   private:

      static const int n = 624, m = 397; // compile time constants

      // The state array and variable p were static members in the C++
      // port by J. Bedaux. Made non-static by David Morse, Dec. 2011

      // state vector array
      unsigned long state[n];

      // position in state array
      int  p;

      // true if init function is called
      bool init;

      // Used by gen_state()
      unsigned long twiddle(unsigned long, unsigned long);

      // Generate new state
      void gen_state();

   };

   // Inline methods.

   inline
   unsigned long MTRand_int32::twiddle(unsigned long u, unsigned long v)
   {
      return (((u & 0x80000000UL) | (v & 0x7FFFFFFFUL)) >> 1)
        ^ ((v & 1UL) ? 0x9908B0DFUL : 0x0UL);
   }

   /*
   * Generate 32 bit random int, private method.
   */
   inline unsigned long MTRand_int32::rand_int32()
   {
      if (p == n) gen_state(); // new state vector needed
      // gen_state() is split off to be non-inline, because it is only called 
      // once in every 624 calls and otherwise irand() would become too big 
      // to get inlined
      unsigned long x = state[p++];
      x ^= (x >> 11);
      x ^= (x << 7) & 0x9D2C5680UL;
      x ^= (x << 15) & 0xEFC60000UL;
      return x ^ (x >> 18);
   }

   /*
   * Serialize to/from an archive.
   */
   template <class Archive>
   void MTRand_int32::serialize(Archive& ar, const unsigned int version)
   {
      for (int i=0; i < n; ++i) {
         ar & state[i];
      }
      ar & p;
      ar & init;
   }

   inline unsigned long MTRand_int32::operator()() { return rand_int32(); }

   class MTRand : public MTRand_int32
   {
   public:

      MTRand() : MTRand_int32() {}

      MTRand(unsigned long seed) : MTRand_int32(seed) {}

      MTRand(const unsigned long* seed, int size) : MTRand_int32(seed, size) {}

      MTRand(const MTRand_int32& other) : MTRand_int32(other) {}

      ~MTRand() {}

      double operator()()
      {  return static_cast<double>(rand_int32()) * (1. / 4294967296.); }
      // divided by 2^32

   //private:

      //MTRand(const MTRand&); // copy constructor not defined
      //void operator=(const MTRand&); // assignment operator not defined

   };

   class MTRand_closed : public MTRand_int32
   {
   public:

      MTRand_closed()
       : MTRand_int32()
      {}

      MTRand_closed(unsigned long seed)
       : MTRand_int32(seed)
      {}

      MTRand_closed(const unsigned long* seed, int size)
       : MTRand_int32(seed, size)
      {}

      ~MTRand_closed() {}

      double operator()()
      {
        return static_cast<double>(rand_int32()) * (1. / 4294967295.); }
        // divided by 2^32 - 1

   //private:

      //MTRand_closed(const MTRand_closed&); // copy constructor not defined

      //void operator=(const MTRand_closed&); // assignment operator not defined
   };

   class MTRand_open : public MTRand_int32
   {

   public:

      MTRand_open()
       : MTRand_int32()
      {}

      MTRand_open(unsigned long seed)
       : MTRand_int32(seed)
      {}

      MTRand_open(const unsigned long* seed, int size)
       : MTRand_int32(seed, size)
      {}

      ~MTRand_open() {}

      double operator()()
      {
        return (static_cast<double>(rand_int32()) + .5)*(1./4294967296.);
        // divided by 2^32
      }

   //private:
      //MTRand_open(const MTRand_open&); // copy constructor not defined
      //void operator=(const MTRand_open&); // assignment operator not defined
   };

   class MTRand53 : public MTRand_int32
   {
   public:

      MTRand53()
       : MTRand_int32()
      {}

      MTRand53(unsigned long seed)
       : MTRand_int32(seed)
      {}

      MTRand53(const unsigned long* seed, int size)
       : MTRand_int32(seed, size)
      {}

      ~MTRand53() {}
      double operator()()
      {
        return (static_cast<double>(rand_int32() >> 5) * 67108864. +
          static_cast<double>(rand_int32() >> 6)) * (1. / 9007199254740992.);
      }

   //private:

      //MTRand53(const MTRand53&); // copy constructor not defined
      //void operator=(const MTRand53&); // assignment operator not defined

   };

}
#endif // MTRAND_H