random.hpp

/*
 The C++ standard library random stuff can be a bit bonkers(*) at times,
 so we have this custom class which has a decent enough random base generator
 and some simple methods for getting values out as floats etc...

(*) See for example the Microsoft implementation of std::uniform_int_distribution...
Or the Cauchy distribution, which won't allow a scale factor of 0 to be used, while
useful for our audio/music applications as a special case.
*/

struct Xoroshiro128Plus
{
    // have some non-zero init state to avoid the zero init state problem
    // which would cause only zeros to be produced
    uint64_t state[2] = {4294967311, 100007};
    Xoroshiro128Plus()
    {
        // experimentally known that after seeding, useful to advance the state,
        // otoh this might be optimized out by the compiler...?
        operator()();
    }
    Xoroshiro128Plus(uint64_t s1, uint64_t s2) : state{s1, s2} { operator()(); }
    void seed(uint64_t s0, uint64_t s1)
    {
        state[0] = s0;
        state[1] = s1;
        operator()();
    }

    bool isSeeded() { return state[0] || state[1]; }

    static uint64_t rotl(uint64_t x, int k) { return (x << k) | (x >> (64 - k)); }

    uint64_t operator()()
    {
        uint64_t s0 = state[0];
        uint64_t s1 = state[1];
        uint64_t result = s0 + s1;

        s1 ^= s0;
        state[0] = rotl(s0, 55) ^ s1 ^ (s1 << 14);
        state[1] = rotl(s1, 36);

        return result;
    }
    constexpr uint64_t min() const { return 0; }
    constexpr uint64_t max() const { return UINT64_MAX; }
    double nextFloat64() { return (*this)() * 5.421010862427522e-20; }
    uint64_t nextUint64() { return (*this)(); }
    uint32_t nextUint32()
    {
        // Take top 32 bits which has better randomness properties
        return operator()() >> 32;
    }
    float nextFloat() { return nextUint32() * 2.32830629e-10f; }
    float nextFloatInRange(float minvalue, float maxvalue)
    {
        return mapvalue(nextFloat(), 0.0f, 1.0f, minvalue, maxvalue);
    }
    double nextFloat64InRange(double minvalue, double maxvalue)
    {
        return mapvalue(nextFloat64(), 0.0, 1.0, minvalue, maxvalue);
    }
    int nextInt32InRange(int minval, int maxval)
    {
        assert(maxval > minval);
        return minval + (nextUint32() % (maxval - minval));
    }
    double nextCauchy(double location, double scale)
    {
        double z = nextFloat64();
        return location + scale * std::tan(M_PI * (z - 0.5));
    }
    // pretty good substitute for Gauss
    double nextHypCos(double location, double scale)
    {
        // we can't do the final calculation with exactly 0.0 or 1.0, so clamp
        // there might be some other ways to deal with this, but this shall suffice for now
        double z = std::clamp(nextFloat64(), std::numeric_limits<double>::epsilon(),
                              1.0 - std::numeric_limits<double>::epsilon());
        return location + scale * (2.0 / M_PI * std::log(std::tan(M_PI / 2.0 * z)));
    }
};