prng.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_PRNG_H
#define THEORETICA_PRNG_H
 
#include <vector>
#include "./pseudorandom.h"
#include "../core/error.h"
 
 
namespace theoretica {
 
    class PRNG {
        private:
            
            pseudorandom_function f;
 
            uint64_t x;
 
            std::vector<uint64_t> param;
 
        public:
 
            PRNG(pseudorandom_function p,
                uint64_t seed,
                const std::vector<uint64_t>& s) : f(p), x(seed), param(s) {}
 
 
            PRNG(pseudorandom_function p,
                const std::vector<uint64_t>& s) : f(p), x(1), param(s) {}
 
 
            PRNG(pseudorandom_function p, uint64_t seed) : f(p), x(seed) {}
 
 
            PRNG(uint64_t seed) {
                *this = PRNG::xoshiro(seed);
            }
 
 
            inline void seed(uint64_t seed) {
                x = seed;
            }
 
 
            inline uint64_t next() {
                return x = f(x, param);
            }
 
 
            inline uint64_t operator()() {
                return next();
            }
 
 
            inline void discard(uint64_t n) {
                for (uint64_t i = 0; i < n; ++i)
                    next();
            }
 
 
            inline uint64_t last() const {
                return x;
            }
 
 
            inline void set_function(pseudorandom_function p) {
                f = p;
            }
 
            inline pseudorandom_function get_function() const {
                return f;
            }
 
 
            inline void set_param(const std::vector<uint64_t>& v) {
                param = v;
            }
 
            inline void set_param(unsigned int i, uint64_t value) {
                param[i] = value;
            }
 
            inline std::vector<uint64_t> get_param() const {
                return param;
            }
 
 
            inline PRNG& operator>>(uint64_t& n) {
                n = next();
                return *this;
            }
 
 
            inline static PRNG linear_congruential(uint64_t seed = 1) {
 
                if(seed == 0)
                    seed = 1;
 
                return PRNG(
                    randgen_congruential, seed,
                    {48271, 0, ((uint64_t) 1 << 31) - 1}
                );
            }
 
 
            inline static PRNG xoshiro(const std::vector<uint64_t>& p) {
                return PRNG(randgen_xoshiro, 0, p);
            }
 
 
            inline static PRNG xoshiro(uint64_t seed = 1) {
 
                if(seed == 0)
                    seed = 1;
                
                const uint64_t n1 = randgen_splitmix64(seed);
                const uint64_t n2 = randgen_splitmix64(n1);
                const uint64_t n3 = randgen_splitmix64(n2);
                const uint64_t n4 = randgen_splitmix64(n3);
 
                return PRNG(randgen_xoshiro, 0, {n1, n2, n3, n4});
            }
 
 
            inline static PRNG splitmix64(uint64_t seed = 1) {
 
                if(seed == 0)
                    seed = 1;
 
                return PRNG(randgen_splitmix64, seed);
            }
 
 
            inline static PRNG wyrand(uint64_t seed = 1,
                uint64_t p1 = 2549536629329, uint64_t p2 = 136137137) {
 
                if(seed == 0)
                    seed = 1;
 
                if(p2 == 0)
                    p2 = randgen_splitmix64(seed);
 
                return PRNG(randgen_wyrand, 0, {seed, p1, p2});
            }
 
 
            inline static PRNG middlesquare(
                uint64_t seed, uint64_t offset = 765872292751861) {
 
                if(seed == 0)
                    seed = randgen_splitmix64(765872292751861);
 
                if(offset == 0)
                    offset = 765872292751861;
 
                return PRNG(randgen_middlesquare, seed, {offset});
            }
        
    };
 
 
    template<typename Vector>
    inline void shuffle(Vector& v, PRNG& g, unsigned int rounds = 0) {
 
        if(v.size() == 0)
            TH_MATH_ERROR("shuffle", v.size(), INVALID_ARGUMENT);
 
        // The number of rounds defaults to (N - 1)^2
        if(rounds == 0)
            rounds = square(v.size() - 1);
 
        for (unsigned int i = 0; i < rounds; ++i) {
            
            // Generate two random indices
            size_t index1 = g() % v.size();
            size_t index2 = g() % v.size();
 
            // Exchange the two values at the random indices
            const auto x1 = v[index1];
            v[index1] = v[index2];
            v[index2] = x1;
        }
    }
 
}
 
#endif