eoRNG.h

#ifndef EO_RANDOM_NUMBER_GENERATOR
#define EO_RANDOM_NUMBER_GENERATOR

# if (defined _MSC_VER)

typedef unsigned long uint32_t;
#else
#if (! defined __sun)
// The C99-standard defines uint32_t to be declared in stdint.h, but some
// systems don't have that and implement it in inttypes.h.
#include <stdint.h>
#else
#include <inttypes.h>
#endif
#endif

#include <cmath>
#include <vector>
#include "eoPersistent.h"
#include "eoObject.h"


class eoRng : public eoObject, public eoPersistent
{
public :

    eoRng(uint32_t s)
        : state(0), next(0), left(-1), cached(false)
        {
            state = new uint32_t[N+1];
            initialize(2*s);
        }

    ~eoRng()
        {
            delete [] state;
        }

    void reseed(uint32_t s)
        {
            initialize(2*s);
        }

    /* FIXME remove in next release
    ** Re-initializes the Random Number Generator

    This is the traditional seeding procedure. This version is deprecated and
    only provided for compatibility with old code. In new projects you should
    use reseed.

    @see reseed for details on usage of the seeding value.

    @version old version (deprecated)
    *
    void oldReseed(uint32_t s)
        {
            initialize(s);
        }
    */

    double uniform(double m = 1.0)
        { // random number between [0, m]
            return m * double(rand()) / double(1.0 + rand_max());
        }

    double uniform(double min, double max)
        { // random number between [min, max]
            return min + uniform(max - min);
        }

    uint32_t random(uint32_t m)
        {
            // C++ Standard (4.9 Floatingintegral conversions [conv.fpint])
            // defines floating point to integer conversion as truncation
            // ("rounding towards zero"): "An rvalue of a floating point type
            // can be converted to an rvalue of an integer type. The conversion
            // truncates; that is, the fractional part is discarded"
            return uint32_t(uniform() * double(m));
        }

    bool flip(double bias=0.5)
        {
            return uniform() < bias;
        }

    double normal();

    double normal(double stdev)
        { return stdev * normal(); }

    double normal(double mean, double stdev)
        { return mean + normal(stdev); }

    double negexp(double mean)
        {
            return -mean*log(double(rand()) / rand_max());
        }

    uint32_t rand();

    uint32_t rand_max() const { return uint32_t(0xffffffff); }

    template <typename TYPE>
    int roulette_wheel(const std::vector<TYPE>& vec, TYPE total = 0)
        {
            if (total == 0)
            { // count
                for (unsigned    i = 0; i < vec.size(); ++i)
                    total += vec[i];
            }
            double fortune = uniform() * total;
            int i = 0;
            while (fortune >= 0)
            {
                fortune -= vec[i++];
            }
            return --i;
        }


    template <typename TYPE>
    const TYPE& choice(const std::vector<TYPE>& vec)
        { return vec[random(vec.size())]; }


    template <typename TYPE>
    TYPE& choice(std::vector<TYPE>& vec)
        { return vec[random(vec.size())]; }

    void printOn(std::ostream& _os) const
        {
            for (int i = 0; i < N; ++i)
            {
                _os << state[i] << ' ';
            }
            _os << int(next - state) << ' ';
            _os << left << ' ' << cached << ' ' << cacheValue;
        }

    void readFrom(std::istream& _is)
        {
            for (int i = 0; i < N; ++i)
            {
                _is >> state[i];
            }

            int n;
            _is >> n;
            next = state + n;

            _is >> left;
            _is >> cached;
            _is >> cacheValue;
        }

    std::string className() const { return "Mersenne-Twister"; }

private:

    uint32_t restart();

    /* @brief Initialize state

    We initialize state[0..(N-1)] via the generator

    <tt>x_new = (69069 * x_old) mod 2^32</tt>

    from Line 15 of Table 1, p. 106, Sec. 3.3.4 of Knuth's _The Art of Computer
    Programming_, Volume 2, 3rd ed.

    Notes (SJC): I do not know what the initial state requirements of the
    Mersenne Twister are, but it seems this seeding generator could be better.
    It achieves the maximum period for its modulus (2^30) iff x_initial is odd
    (p. 20-21, Sec. 3.2.1.2, Knuth); if x_initial can be even, you have
    sequences like 0, 0, 0, ...; 2^31, 2^31, 2^31, ...; 2^30, 2^30, 2^30, ...;
    2^29, 2^29 + 2^31, 2^29, 2^29 + 2^31, ..., etc. so I force seed to be odd
    below.

    Even if x_initial is odd, if x_initial is 1 mod 4 then

    the          lowest bit of x is always 1,
    the  next-to-lowest bit of x is always 0,
    the 2nd-from-lowest bit of x alternates      ... 0 1 0 1 0 1 0 1 ... ,
    the 3rd-from-lowest bit of x 4-cycles        ... 0 1 1 0 0 1 1 0 ... ,
    the 4th-from-lowest bit of x has the 8-cycle ... 0 0 0 1 1 1 1 0 ... ,
    ...

    and if x_initial is 3 mod 4 then

    the          lowest bit of x is always 1,
    the  next-to-lowest bit of x is always 1,
    the 2nd-from-lowest bit of x alternates      ... 0 1 0 1 0 1 0 1 ... ,
    the 3rd-from-lowest bit of x 4-cycles        ... 0 0 1 1 0 0 1 1 ... ,
    the 4th-from-lowest bit of x has the 8-cycle ... 0 0 1 1 1 1 0 0 ... ,
    ...

    The generator's potency (min. s>=0 with (69069-1)^s = 0 mod 2^32) is 16,
    which seems to be alright by p. 25, Sec. 3.2.1.3 of Knuth. It also does well
    in the dimension 2..5 spectral tests, but it could be better in dimension 6
    (Line 15, Table 1, p. 106, Sec. 3.3.4, Knuth).

    Note that the random number user does not see the values generated here
    directly since restart() will always munge them first, so maybe none of all
    of this matters. In fact, the seed values made here could even be
    extra-special desirable if the Mersenne Twister theory says so-- that's why
    the only change I made is to restrict to odd seeds.
    */
    void initialize(uint32_t seed);

    uint32_t *state;

    uint32_t *next;

    int left;

    bool cached;

    double cacheValue;

    static const int N;

    static const int M;

    static const uint32_t K;


    eoRng(const eoRng&);

    eoRng& operator=(const eoRng&);
};
namespace eo
{
    extern eoRng rng;
}
using eo::rng;

// Implementation of some eoRng members.... Don't mind the mess, it does work.

#define hiBit(u)       ((u) & 0x80000000U)   // mask all but highest   bit of u
#define loBit(u)       ((u) & 0x00000001U)   // mask all but lowest    bit of u
#define loBits(u)      ((u) & 0x7FFFFFFFU)   // mask     the highest   bit of u
#define mixBits(u, v)  (hiBit(u)|loBits(v))  // move hi bit of u to hi bit of v


inline void eoRng::initialize(uint32_t seed)
{
    left = -1;

    register uint32_t x = (seed | 1U) & 0xFFFFFFFFU, *s = state;
    register int j;

    for(left=0, *s++=x, j=N; --j;
        *s++ = (x*=69069U) & 0xFFFFFFFFU) ;
}



inline uint32_t eoRng::restart()
{
    register uint32_t *p0=state, *p2=state+2, *pM=state+M, s0, s1;
    register int j;

    left=N-1, next=state+1;

    for(s0=state[0], s1=state[1], j=N-M+1; --j; s0=s1, s1=*p2++)
        *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);

    for(pM=state, j=M; --j; s0=s1, s1=*p2++)
        *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);

    s1=state[0], *p0 = *pM ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);
    s1 ^= (s1 >> 11);
    s1 ^= (s1 <<  7) & 0x9D2C5680U;
    s1 ^= (s1 << 15) & 0xEFC60000U;
    return(s1 ^ (s1 >> 18));
}



inline uint32_t eoRng::rand()
{
    if(--left < 0)
        return(restart());
    uint32_t y  = *next++;
    y ^= (y >> 11);
    y ^= (y <<  7) & 0x9D2C5680U;
    y ^= (y << 15) & 0xEFC60000U;
    return(y ^ (y >> 18));
}



inline double eoRng::normal()
{
    if (cached) {
        cached = false;
        return cacheValue;
    }
    double rSquare, var1, var2;
    do {
        var1 = 2.0 * uniform() - 1.0;
        var2 = 2.0 * uniform() - 1.0;
        rSquare = var1 * var1 + var2 * var2;
    } while (rSquare >= 1.0 || rSquare == 0.0);
    double factor = sqrt(-2.0 * log(rSquare) / rSquare);
    cacheValue = var1 * factor;
    cached = true;
    return (var2 * factor);
}



namespace eo
{
    template <typename T>
    inline T random(const T& min, const T& max) {
        return static_cast<T>(rng.uniform() * (max-min)) + min; }

    template <typename T>
    inline T random(const T& max) {
        return static_cast<T>(rng.uniform() * max); }

    inline double normal() { return rng.normal(); }
}


#endif


// Local Variables:
// coding: iso-8859-1
// mode: C++
// c-file-offsets: ((c . 0))
// c-file-style: "Stroustrup"
// fill-column: 80
// End: