add square with carry

[lemire]

Stefan Kanthak wrote:

While Lehmer’s multiplicative congruential generator is fast, it does NOT pass PractRand 0.94. In [does it beat the minimal standard](url
http://www.pcg-random.org/posts/does-it-beat-the-minimal-standard.html) and [too big to fail](url
http://www.pcg-random.org/posts/too-big-to-fail.html) Melissa O’Neill tested these generators with PractRand 0.93, where they passed.

Consider an idea from George Marsaglia instead: save and add the “carry” of the multiplication, giving a multiply-with-carry generator:

uint64_t seed = 0, carry = 0x...;
uint64_t square_with_carry() {
__uint128_t mwc = (__uint128_t) seed * 0x... + carry;
seed = mwc;
carry = square >> 64;
return seed;
}

uint64_t seed = 0, carry = 0xb5ad4eceda1ce2a9;
uint64_t square_with_carry() {
__uint128_t square = (__uint128_t) seed * seed + carry;
seed = square;
carry = square >> 64;
return seed;
}

This 64-bit generator passes the PractRand test suite (version 0.94) at least up to 4 TB!
See [source](url
https://godbolt.org/z/qFfRYh) and [source](url
https://godbolt.org/z/4jf1i_) for the code generated by GCC 8.2 and clang 7.0 respectively, plus [VC source](url
https://godbolt.org/z/9RRnBL) for an implementation using Microsofts Visual C compiler (which does not support a 128 bit integer type).
I used the latter for the PractRand test.