ranx.hpp

Go to the documentation of this file.

#ifndef RANX_HPP
#define RANX_HPP
#include<vector>
#include <bitset>
/* Thread-safe random number generator */
namespace Random{

    // 3 ints store the state
    struct rands{
        int s[3];
    };

    void init_seeds_ext(rands& sv){
        if ((sv.s[0] & -2) == 0)  sv.s[0] -= 1023;
        if ((sv.s[1] & -8) == 0)  sv.s[1] -= 1023;
        if ((sv.s[2] & -16) == 0) sv.s[2] -= 1023;
    }

    double taus88_ext(rands& sv){
        // Generates a random number between 0 and 1.
        // Reference:
        // L'Ecuyer, P. (1996) `Maximally equidistributed combined Tausworthe
        // generators', Math. of Comput., 65, 203-213.

        // The cycle length is claimed to be about 2^(88) or about 3E+26.
        // Actually - (2^31 - 1).(2^29 - 1).(2^28 - 1).

        uint32_t i1 = sv.s[0];
        uint32_t i2 = sv.s[1];
        uint32_t i3 = sv.s[2];
        uint32_t b  = (((i1 << 13) ^ i1) >> 19);

        i1 = (((i1 &  -2) << 12) ^  b);
        b  = (((i2 <<  2) ^  i2) >> 25);
        i2 = (((i2 &  -8) <<  4) ^  b);
        b  = (((i3 <<  3) ^  i3) >> 11);
        i3 = (((i3 & -16) << 17) ^  b);
        sv.s[0] = i1;
        sv.s[1] = i2;
        sv.s[2] = i3;

        return int( (i1 ^ i2) ^ i3) * 2.3283064365E-10 + 0.5;
    }

    // Global random number pool
    std::vector<rands> sv;

    // Initialize seed
    void init_ranx(int nthreads, int mype){
        // Allocate
        sv = std::vector<rands>(nthreads);

        // Set initial seeds
        for (int i = 0; i<nthreads; i++){
            int gi = mype*nthreads + i;
            sv[i].s[0] = 1234*gi;
            sv[i].s[1] = 2345*gi + 6789;
            sv[i].s[2] = 4321*gi + 10;
            init_seeds_ext( sv[i] );
        }
    }

    // random number generator
    double ranx(int thread_id){
        return taus88_ext( sv[thread_id] );
    }
}

#endif