Hello,

I observed what seems to be an unexpected behaviour in the jumpedmethod of the PRNG MT19937.

According to the official documentation calling MT19937.jumped() is supposed to:

Returns a new bit generator with the state jumped
The state of the returned big generator is jumped as-if 2**(128 * jumps) random numbers have been generated.

According to this spec I would expect jump > random_raw and random_raw > jump to produce the same generators. Indeed jump > random_raw should offest the generator by 2**128 then by 1, and random_raw > jump should offset the generator by 1 then by 2**128 which should be equivalent.

But when running the snippet below it does not seem to be the case.

Reproducing code example:

from numpy.random import MT19937

rng = MT19937(42)
rng.random_raw()   # value given state offset of 0
rng.random_raw()   # value given state offset of 1
rng = rng.jumped()
rng.random_raw()   # value given state offset of 2**128 +2
# Got: 1016764315

rng = MT19937(42)
rng.random_raw()   # value given state offset of 0
rng = rng.jumped()
rng.random_raw()   # value given state offset of 2**128 +1
rng.random_raw()   # value given state offset of 2**128 +2
# Got: 1769577901

NOTE: I called random_raw before calling jumped in order not to fall into the strange case covered by the comment https://github.com/numpy/numpy/blob/master/numpy/random/src/mt19937/mt19937-jump.c#L193

NOTE: I have not checked - for the moment - whether the jump implementation at http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/JUMP/index.html has the issue or not.

Error message:

Numpy/Python version information:

1.18.1 3.8.1 (tags/v3.8.1:1b293b6, Dec 18 2019, 22:39:24) [MSC v.1916 32 bit (Intel)]