Which C++ random number engines have a O(1) discard function?

Viewed 834

Since C++11 there are a number of std random number engines. One of the member functions they implement is void discard(int long long z) which skips over z randomly generated numbers. The complexity of this function is given as O(z) on www.cplusplus.com (http://www.cplusplus.com/reference/random/mersenne_twister_engine/discard/)

However, on www.cppreference.com (http://en.cppreference.com/w/cpp/numeric/random/mersenne_twister_engine/discard) there is a note to say that

For some engines, "fast jump" algorithms are known, which advancing the state by many steps (order of millions) without calculating intermediate state transitions.

How do I know for which engines the actual cost of discard is O(1)?

4 Answers

Well, if you use precomputed jump points, O(1) will work for each and every RNG in existence. Please, remember, that there are algorithm which might have better than O(z), but not O(1) - say, O(log2 z).

If we're talking about jump to arbitrary point, things get interesting. For example, for linear congruential generator there is known O(log2 z) jump ahead algorithm, based upon paper by F. Brown, "Random Number Generation with Arbitrary Stride," Trans. Am. Nucl. Soc. (Nov. 1994). Code example is here.

There is LCG RNG in the C++11 standard, not sure how fast jump ahead is done in particular implementation (http://en.cppreference.com/w/cpp/numeric/random/linear_congruential_engine)

PCG family of RNGs share the same property, I believe

Fact is that std::linear_congruential_engine<UIntType,a,c,m>::discard(unsigned long long z) is definitely possible to implement very efficiently. It is mostly equivalent to exponentiation of a to the power of z modulo m (for both zero and non-zero c) - it means most basic software implementation executes in O(log(z % phi(m))) UIntType multiplications (phi(m)=m-1 for prime number and <m in general) and can be implemented to execute much faster with parallel hardware exponentiation algorithm.

^^^NOTE that in a way O(log(z % phi(m))) is O(1) because log2(z % phi(m)) < log2(m) < sizeof(UIntType)*CHAR_BIT - though in practice it is more often like O(log(z)).

Also probably there are efficient algorithms for most other engines' discard functions that would satisfy O(P(size of state)) constraint (P(x) - some low degree polynomial function, most likely 1+epsilon degree or something even smaller like x*log(x), given that log(z) < sizeof(unsigned long long)*CHAR_BIT may be considered as constant).

Somehow for some unknown reason C++ standard (as of ISO/IEC 14882:2017) does not require to implement discard in more efficient way than just z operator()() calls for any PRNG engine including those that definitely allow this.

To me personally it is baffling and JUST MAKES NO SENSE - it brutally violates one of the fundamental C++ language design principles that is to add to the C++ standard only reasonable functionality in terms of performance and practical usefulness.

For example: there is NO SUCH THING AS std::list<T>::operator[](size_type n) even though it is as "easy" as to just call operator++() n times with iterator begin(). And naturally so because O(n) execution time would make this function unreasonable choice in any practical application (a code word for "plain stupid idea"). For this obvious reason a[n] and a.at(n) is not part of mandatory Sequence container requirements (ISO/IEC 14882:2017 26.2.3 Table 87) but instead is a part of Optional sequence container operations (ISO/IEC 14882:2017 26.2.3 Table 88).

SO why in the world then e.discard(z) is part of mandatory Random number engine requirements (ISO/IEC 14882:2017 29.6.1.4 Table 104) with this ridiculous complexity requirement - no worse than the complexity of z consecutive calls e() - instead of some optional operations section entry with adequate complexity requirement like O(size of state) or O(P(size of state))?

Even more baffling was to actually find in my GCC this real-world implementation:

void discard(unsigned long long __z)
{
    for (; __z != 0ULL; --__z)
        (*this)(); //<--  Wait what? Are you kidding me?
}

So once again as it was before we have no other choice than to implement the necessary functionality ourselves...

I don't think such things exists at all. My heuristic conclusion is that O(1)-jump RNG is basically a hash, with all that this implies (e.g. it might not be "good" RNG at all). (See new answer by @AkiSuihkonen and link therein.)

But even if you are asking about O(log z) I don't see that implemented in the STL. In GCC's all the discard functions I was able to grep are all simple loops.

      discard(unsigned long long __z)
      {
        for (; __z != 0ULL; --__z)
          (*this)();
      }

Which is not only sad but also misleading since discard should exists only if there is an efficient way to do it.

The only non trivial one is mersenne (below) but it is still O(z).

    discard(unsigned long long __z)
    {
      while (__z > state_size - _M_p)
    {
      __z -= state_size - _M_p;
      _M_gen_rand();
    }
      _M_p += __z;
    }

Boost's Mersenne, has a skip function but it is only called for skips larger than a (default of) 10000000 (!?). Which already tells me that that the skip is very heavy computationally (even if it is O(log z)). https://www.boost.org/doc/libs/1_72_0/boost/random/mersenne_twister.hpp

Finally, Thrust has an efficient discard for linear congruential apparently, but only in the case c == 0. (Which I am not sure if it makes it less useful as a RNG.) https://thrust.github.io/doc/classthrust_1_1random_1_1linear__congruential__engine_aec05b19d2a85d02f1ff437791ea4dd68.html#aec05b19d2a85d02f1ff437791ea4dd68

All Counter Based Random Number Generators

These work solely by computing a function uint64_t rnd(uint64_t counter) or maybe uint16_t rnd(uint128_t counter). Then the skip function is as easy as

// the method itself is "randomly" generated -- known at least
// by von Neumann to typically generate poor results
struct MyRandomCBRNG {
    uint64_t counter{0};
    void skip(uint64_t a) { counter+=a;}

    uint64_t operator()() {
        uint64_t x = counter++;
        // repeat multiple times if needed
        x = x * 0xdeadbeefcafebabeull ^ (x >> 53) ^ (x << 11) + 0x13459876abcdfdecull;
        return x;
    }
};

One can make even cryptographically strong CBRNGs by hashing the counter concatenated by secret key with something like SHA-512, not to mention Blum-Blum-Shub.

Related