PseudoRandom.h

// PseudoRandom.h
//
// Random number generator.  From Numerical recipes in C

#pragma once
#ifndef Ceda_cxUtils_PseudoRandom_H
#define Ceda_cxUtils_PseudoRandom_H

#include "BasicTypeSizes.h"
#include "cxUtils.h"
#include <mutex>
#include <type_traits>
#include <limits>

/*
todo:  use the following xorshift128+ generator, which uses 128 bits of state and has a maximal 
period of (2^128)-1. It passes BigCrush, but not when reversed:

    #include <stdint.h>

    struct xorshift128p_state 
    {
        uint64_t a, b;
    };

    // The state must be seeded so that it is not all zero
    uint64_t xorshift128p(xorshift128p_state* state)
    {
        uint64_t t = state->a;
        uint64_t const s = state->b;
        state->a = s;
        t ^= t << 23;		// a
        t ^= t >> 17;		// b
        t ^= s ^ (s >> 26);	// c
        state->b = t;
        return t + s;
    }

This generator is one of the fastest generators passing BigCrush.  It is generally superior
to anything in <random>!
*/

namespace ceda
{

class xorshift
{
public:
    xorshift(uint32 seed = 1) : state(seed) {}    

    uint32 get()
    {
	    uint32 x = state;
	    x ^= x << 13;
	    x ^= x >> 17;
	    x ^= x << 5;
	    state = x;
	    return x;
    }

    // get in range [0,n) where n > 0
    int get(int n)
    {
        return get() % (uint32)n;
    }

    // get in range [x1,x2) which must be a nonempty range
    int get(int x1, int x2)
    {
        uint32 n = x2-x1;
        return x1 + (get() % n);
    }

private:
    uint32 state;
};


///////////////////////////////////////////////////////////////////////////////////////////////////
// LinearCongruentialGenerator

class LinearCongruentialGenerator
{
public:
    LinearCongruentialGenerator(uint32 seed = 0) : x(seed) {}    

    uint32 get()
    {
        x = (x * 1664525) + 1013904223;
        return x;
    }

    // get in range [0,n) where n > 0
    int get(int n)
    {
        x = (x * 1664525) + 1013904223;
        return x % (uint32)n;
    }

    // get in range [x1,x2) which must be a nonempty range
    int get(int x1, int x2)
    {
        x = (x * 1664525) + 1013904223;
        uint32 n = x2-x1;
        return x1 + (x % n);
    }

private:
    uint32 x;
};

///////////////////////////////////////////////////////////////////////////////////////////////////
// RandomNumGen

class cxUtils_API RandomNumGen
{
public:
    RandomNumGen(bool seedFromClock = true);
    ~RandomNumGen();

    // Generate a uniformly distributed random number between [0.0, 1.0]
    double GetUniformDistDouble();

    // Generate a uniformly distributed random number in range [x1,x2]
    double GetUniformDistDouble(double x1, double x2);

    // Generate a uniformly distributed random integer in range [x1,x2)
    int32 GetUniformDistInteger(int32 x1, int32 x2);
    int64 GetUniformDistInteger(int64 x1, int64 x2);

    template <typename T, std::enable_if_t<std::is_integral_v<T>, int> = 0>
    T GetUniformDistInteger(T x1, T x2)
    {
        if constexpr (std::numeric_limits<T>::digits <= 32)
        {
            return static_cast<T>(GetUniformDistInteger(static_cast<int32>(x1),static_cast<int32>(x2)));
        }
        else
        {
            return static_cast<T>(GetUniformDistInteger(static_cast<int64>(x1),static_cast<int64>(x2)));
        }
    }

    ssize_t GetUniformDistInteger_ssize_t(ssize_t x1, ssize_t x2)
    {
        return GetUniformDistInteger(x1,x2);
    }

    // Generate a random number which is normally distributed with mean = 0.0 and stddevn = 1.0
    double GetGaussianDistDouble();

private:
    // For uniform number generation
    long ix1, ix2, ix3;
    double* r;

    // For gaussian number generation
    int iset;
    double gset;
};

///////////////////////////////////////////////////////////////////////////////////////////////////
// ThreadSafeRandomNumGen

class ThreadSafeRandomNumGen
{
public:
    ThreadSafeRandomNumGen(bool seedFromClock = true) : rng_(seedFromClock) {}

    double GetUniformDistDouble()
    {
        std::lock_guard<std::mutex> lock(mutex_);
        return rng_.GetUniformDistDouble();
    }

    double GetUniformDistDouble(double x1, double x2)
    {
        std::lock_guard<std::mutex> lock(mutex_);
        return rng_.GetUniformDistDouble(x1,x2);
    }

    int32 GetUniformDistInteger(int32 x1, int32 x2)
    {
        std::lock_guard<std::mutex> lock(mutex_);
        return rng_.GetUniformDistInteger(x1,x2);
    }

    int64 GetUniformDistInteger(int64 x1, int64 x2)
    {
        std::lock_guard<std::mutex> lock(mutex_);
        return rng_.GetUniformDistInteger(x1,x2);
    }

    template <typename T, std::enable_if_t<std::is_integral_v<T>, int> = 0>
    T GetUniformDistInteger(T x1, T x2)
    {
        if constexpr (std::numeric_limits<T>::digits <= 32)
        {
            return static_cast<T>(GetUniformDistInteger(static_cast<int32>(x1),static_cast<int32>(x2)));
        }
        else
        {
            return static_cast<T>(GetUniformDistInteger(static_cast<int64>(x1),static_cast<int64>(x2)));
        }
    }

    ssize_t GetUniformDistInteger_ssize_t(ssize_t x1, ssize_t x2)
    {
        std::lock_guard<std::mutex> lock(mutex_);
        return rng_.GetUniformDistInteger_ssize_t(x1,x2);
    }

    double GetGaussianDistDouble()
    {
        std::lock_guard<std::mutex> lock(mutex_);
        return rng_.GetGaussianDistDouble();
    }

private:
    RandomNumGen rng_;
    mutable std::mutex mutex_;
};

///////////////////////////////////////////////////////////////////////////////////////////////////
// Global, threadsafe functions

// Generate a uniformly distributed random number between [0.0, 1.0]
cxUtils_API double GetUniformDistDouble();

// Generate a uniformly distributed random number in range [x1,x2]
cxUtils_API double GetUniformDistDouble(double x1, double x2);

// Generate a uniformly distributed random integer in range [x1,x2)
cxUtils_API int32 GetUniformDistInteger(int32 x1, int32 x2);
cxUtils_API int64 GetUniformDistInteger(int64 x1, int64 x2);

template <typename T, std::enable_if_t<std::is_integral_v<T>, int> = 0>
T GetUniformDistInteger(T x1, T x2)
{
    if constexpr (std::numeric_limits<T>::digits <= 32)
    {
        return static_cast<T>(GetUniformDistInteger(static_cast<int32>(x1),static_cast<int32>(x2)));
    }
    else
    {
        return static_cast<T>(GetUniformDistInteger(static_cast<int64>(x1),static_cast<int64>(x2)));
    }
}

// Generate a uniformly distributed random integer in range [x1,x2)
cxUtils_API ssize_t GetUniformDistInteger_ssize_t(ssize_t x1, ssize_t x2);

// Generate a random number which is normally distributed with mean = 0.0 and stddevn = 1.0
cxUtils_API double GetGaussianDistDouble();

// todo: make more efficient by using all 64 bits each time
inline void GenRandomBytes( octet_t* dest, ssize_t n )
{
    for( ssize_t i = 0 ; i < n ; ++i, ++dest )
    {
        *dest = (octet_t)GetUniformDistInteger( 0, 255 );
    }    
}

template<typename T>
inline T GetUniformDistVal(T x1, T x2)
{
}

template<>
inline double GetUniformDistVal(double x1, double x2)
{
    return GetUniformDistDouble(x1,x2);
}
template<>
inline int32 GetUniformDistVal(int32 x1, int32 x2)
{
    return GetUniformDistInteger(x1,x2);
}
template<>
inline int64 GetUniformDistVal(int64 x1, int64 x2)
{
    return GetUniformDistInteger(x1,x2);
}

template<typename T>
inline xvector<T> GetUniqueUniformDistVals(T x1, T x2, ssize_t n)
{
    cxAssert( x2 - x1 > n );
    
    xvector<T> ret(n);
    ssize_t i = 0;
    while (i < n)
    {
        T val = GetUniformDistVal<T>(x1, x2);
        
        if (ret.find( val ) == ret.npos)
        {
            ret[i++] = val;
        }
    }
    
    return ret;
}

///////////////////////////////////////////////////////////////////////////////////////////////////
// Random selections with given weightings

cxUtils_API int RandomSelect(RandomNumGen& rng, double w0, double w1);
cxUtils_API int RandomSelect(RandomNumGen& rng, double w0, double w1, double w2);
cxUtils_API int RandomSelect(RandomNumGen& rng, double w0, double w1, double w2,double w3);
cxUtils_API int RandomSelect(RandomNumGen& rng, double w0, double w1, double w2,double w3,double w4);

} // namespace ceda
 
#endif // include guard