pseudorandom.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_PSEUDORANDOM_H
#define THEORETICA_PSEUDORANDOM_H
 
#include <vector>
#include <cstdint>
 
#include "../core/bit_op.h"
#include "../core/error.h"
 
 
namespace theoretica {
 
 
    using pseudorandom_function = uint64_t (*)(uint64_t, std::vector<uint64_t>&);
 
 
    inline uint64_t rand_congruential(
        uint64_t x,
        uint64_t a = 48271,
        uint64_t c = 0,
        uint64_t m = ((uint64_t) 1 << 31) - 1) {
 
        return (a * x + c) % m;
    }
 
 
    inline uint64_t rand_congruential(uint64_t x, std::vector<uint64_t>& state) {
 
        if(state.size() != 3) {
            TH_MATH_ERROR("rand_congruential", state.size(), INVALID_ARGUMENT);
            return 0;
        }
 
        const uint64_t a = state[0];
        const uint64_t c = state[1];
        const uint64_t m = state[2];
 
        if(a > m || c > m) {
            TH_MATH_ERROR("rand_congruential", max(a, c), INVALID_ARGUMENT);
            return 0;
        }
 
        return rand_congruential(x, a, c, m);
    }
 
 
    inline uint64_t rand_xoshiro(uint64_t& a, uint64_t& b, uint64_t& c, uint64_t& d) {
 
        // Add and rotate
        const uint64_t result = bit_rotate(a + d, 23) + a;
        const uint64_t temp = b << 17;
 
        // Shift operations
        c ^= a;
        d ^= b;
        b ^= c;
        a ^= d;
 
        c ^= temp;
        d = bit_rotate(d, 45);
 
        return result;
    }
 
 
    inline uint64_t rand_xoshiro(uint64_t x, std::vector<uint64_t>& state) {
 
        if(state.size() != 4) {
            TH_MATH_ERROR("rand_xoshiro", state.size(), INVALID_ARGUMENT);
            return 0;
        }
 
        return rand_xoshiro(state[0], state[1], state[2], state[3]);
    }
 
 
    inline uint64_t rand_splitmix64(uint64_t x) {
 
        x += 0x9e3779b97f4a7c15;
 
        uint64_t res = x;
        res = (res ^ (res >> 30)) * 0xbf58476d1ce4e5b9;
        res = (res ^ (res >> 27)) * 0x94d049bb133111eb;
 
        return res ^ (res >> 31);
    }
 
 
    inline uint64_t rand_splitmix64(uint64_t x, std::vector<uint64_t>& p) {
        return rand_splitmix64(x);
    }
 
 
    inline uint64_t rand_wyrand(uint64_t& seed, uint64_t p1, uint64_t p2) {
        seed += p1;
        return mix_mum(seed ^ p2, seed);
    }
 
 
    inline uint64_t rand_wyrand(uint64_t x, std::vector<uint64_t>& p) {
 
        if(p.size() != 3) {
            TH_MATH_ERROR("rand_wyrand", p.size(), INVALID_ARGUMENT);
            return 0;
        }
 
        return rand_wyrand(p[0], p[1], p[2]);
    }
 
 
    inline uint64_t rand_middlesquare(uint64_t seed, uint64_t offset = 765872292751861) {
 
        seed += offset;
        uint64_t high, low;
        mul_uint128(seed, seed, low, high);
 
        return (high << 32) | (low >> 32);
    }
 
 
    inline uint64_t rand_middlesquare(uint64_t x, std::vector<uint64_t>& p) {
 
        if(p.size() != 1) {
            TH_MATH_ERROR("rand_middlesquare", p.size(), INVALID_ARGUMENT);
            return 0;
        }
 
        return rand_middlesquare(x, p[0]);
    }
 
}
 
 
#endif