jeudi 21 décembre 2017

Concept-Lite for random distribution engine (Cf. UniformRandomBitGenerator for raw bits)

I have some decades-old code for random number generation code that I wrote based on the work of George Marsaglia. I ran some tests and found that the "Mars" code, using his "Monty Python" scheme for standard normal deviates was 3 1/2 times faster than the std::normal_distribution that comes with Visual C++ 2017. I decided to make that code play nicely with std::random. (A generator for random bits is easy. There is an official concept called UniformRandomBitGenerator that fits the bill perfectly. Create custom distribution compatible with C++11 random header Freebie example below.)

My dj::normal_distribution "Monte Python" class is a drop-in replacement for std::normal_distribution, except I did not fancy implementing all the stuff that no-one will ever use.

I have read that just as UniformRandomBitGenerator is a subset of RandomNumberEngine that contains only the useful stuff, there is

[Q] a concept for a distribution generator that does not have all the tedious (and seldom needed) boilerplate of concept RandomNumberDistribution. Is there one?

I have searched in vain. If there is, I would be happy to publish the Monty Python code that satisfies the concept.

Freebie

namespace dj {

// George Marsaglia's long period (2^60) random number generator
// Satisfies UniformRandomBitsGenerator concept
class Mars_URNG {
public:

    Mars_URNG(uint32_t s1 = 0, uint32_t s2 = 0) {
        seed(s1, s2);
    }

    void seed(uint32_t s1 = 0, uint32_t s2 = 0) {
        z = good(s1); w = good(s2);
    }

    // The following four definitions satisfy the concept
    inline uint32_t operator() () {
        return znew() | wnew();
    }
    constexpr uint32_t min() { return 0; }
    constexpr uint32_t max() { return 0xffffffff; }
    using result_type = uint32_t;

private:
    static const uint32_t arby = 0xA5A5A5A5; // Arbitrary default seed
    // Zero is a baaad seed.
    static uint32_t good(uint32_t s) { return s ? s : arby; }
    uint32_t z;
    uint32_t w;
    inline uint32_t znew() { z = 36969 * (z & 0xffff) + (z >> 16); return z << 16; }
    inline uint32_t wnew() { w = 18000 * (w & 0xffff) + (w >> 16); return w & 0xffff; }
};
}




Aucun commentaire:

Enregistrer un commentaire