Fun with LCG's

lcg.cpp

#include <cstdint>

bool PolyMul(uint32_t from, uint32_t to, uint32_t& out_mul)
{
    for (; ~from & 1; from >>= 1, to >>= 1) {
        if (to & 1)
            return false;
    }

    uint32_t result = 0;

#if 1
    while (uint32_t v = (from * result) ^ to) {
        result |= v & -int32_t(v);
    }
#else
    uint32_t value = 0;

    for (uint32_t i = 0; i < 32; ++i) {
        uint32_t m = UINT32_C(1) << i;

        if ((value ^ to) & m) {
            result |= m;
            value += from << i;
        }
    }
#endif

    out_mul = result;

    return true;
}

void SkippingLCG(uint32_t n, uint32_t& mul, uint32_t& inc)
{
    uint32_t out_mul = 1;
    uint32_t out_inc = 0;

    for (uint32_t i = n, j = mul, k = inc; i; i >>= 1, k += k * j, j *= j) {
        if (i & 1) {
            out_mul *= j;
            out_inc *= j;
            out_inc += k;
        }
    }

    mul = out_mul;
    inc = out_inc;
}

template <uint32_t Mul = 214013, uint32_t Inc = 2531011>
struct LCG {
    uint32_t state;

    LCG(uint32_t seed)
        : state(seed)
    {
    }

    void skip(int32_t n = 1)
    {
        if (n == 0)
            return;

        uint32_t mul = Mul;
        uint32_t inc = Inc;

        if (n < 0) {
            PolyMul(mul, 1, mul);
            inc *= uint32_t(-int32_t(mul));
            n = -n;
        }

        if (n > 1) {
            SkippingLCG(n, mul, inc);
        }

        state = (state * mul) + inc;
    }

    uint32_t operator()(int32_t n = 1)
    {
        uint32_t result = state;
        skip(n);
        return result;
    }
};