eorng.h

这是linux下的进化计算的源代码。 === === === === === === === === === === === ===== check latest news at http:

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

    /** @brief Print RNG

    @param _os Stream to print RNG on
    */
    void printOn(std::ostream& _os) const
        {
            for (int i = 0; i < N; ++i)
            {
                _os << state[i] << ' ';
            }
            _os << int(next - state) << ' ';
            _os << left << ' ' << cached << ' ' << cacheValue;
        }

    /** @brief Read RNG

    @param _is Stream to read RNG from
    */
    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);

    /** @brief Array for the state */
    uint32_t *state;

    /** Pointer to next available random number */
    uint32_t *next;

    /** Number of random numbers currently left */
    int left;

    /** @brief Is there a valid cached value for the normal distribution? */
    bool cached;

    /** @brief Cached value for normal distribution? */
    double cacheValue;

    /** @brief Size of the state-vector */
    static const int N;

    /** Internal constant */
    static const int M;

    /** @brief Magic constant */
    static const uint32_t K;


    /** @brief Copy constructor

    Private copy ctor and assignment operator to make sure that nobody
    accidentally copies the random number generator. If you want similar RNG's,
    make two RNG's and initialize them with the same seed.

    As it cannot be called, we do not provide an implementation.
    */
    eoRng(const eoRng&);

    /** @brief Assignment operator

    @see Copy constructor eoRng(const eoRng&).
    */
    eoRng& operator=(const eoRng&);
};



namespace eo
{
    /** The one and only global eoRng object */
    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
{
    /** @brief Random function

    This is a convenience function for generating random numbers using the
    global eo::rng object.

    Templatized random function, returns a random double in the range [min, max).
    It works with most basic types such as:
    - char
    - int (short, long, signed and unsigned)
    - float, double

    @param min Minimum for distribution
    @param max Maximum for distribution

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

    /** @brief Random function

    @overload

    This is a convenience function for generating random numbers using the
    global eo::rng object.

    Templatized random function, returns a random double in the range [0, max).
    It works with most basic types such as:
    - char
    - int (short, long, signed and unsigned)
    - float, double

    @param max Maximum for distribution

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

    /** Normal distribution

    This is a convenience function for generating random numbers using the
    global eo::rng object.

    @return ormally distributed random number
    */
    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: