I'd probably add another test covering the issue I originally opened:

    { # You already have this one
        "seed": 0,
        "steps": 10,
        "initial": {"key_md5": "64eaf265d2203179fb5ffb73380cd589", "pos":  9},
        "jumped": {"key_md5": "8cb7b061136efceef5217a9ce2cc9a5a", "pos": 598},
    },
    { # But don't have this one
        "seed": 0,
        "steps": 11,
        "initial": {"key_md5": "64eaf265d2203179fb5ffb73380cd589", "pos": 10},
        "jumped": {"key_md5": "8cb7b061136efceef5217a9ce2cc9a5a", "pos": 599},
    },

The added test ensures that we have the following behaviour:
jump(next(mt19937)) == next(jump(mt19937)) - or something close

The only guarantee that I think we should try and make is that the code is faithful to the original author's code. Whether this is true or not I think is a different issue.

If I understand it correctly, the test you propose will not work (in the current implementation / original authors source code) / tried it out and it didn't work:

jump(next(mt19937)) == next(jump(mt19937)) - or something close

However, if this is the case, then I don't think that this pull request closes #15394, does it?
Or put it the other way around: With the current state of the PR, shouldn't the documentation be adopted too?

Current documentation:

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

On this note:

The only guarantee that I think we should try and make is that the code is faithful to the original author's code. Whether this is true or not I think is a different issue.

To what extend is this the policy? I totally understand that you may not have the time to check that all implementation are "true" - but still I think it's worrisome if the random generator behaves in a way we don't understand. So what should the documentation of jump say? If I exaggerate a little: the original author claims that it is as if 2*(128 * jumps) random numbers have been generated; but when we test it, that it not true`

Another note: Shall we implement the same type of test

jump(next(mt19937)) == next(jump(mt19937)) - or something close

in the tests for PCG64? There it should pass; but nice to have this ensured for the future, too.

Depends on what the point of #15394 is. If it is to check the implementation, then I think this PR closes that. If it is to make jumped work in a way different from the original code, then it does not. I don't think the latter is realistic.

My original point when opening the issue was the following:

As a user of the library I saw that the jump was supposed to offset by a fixed number of iterations. So my expectation and the reason why I opened the issue was jump(next(mt19937)) should always be the same as next(jump(mt19937)) - or at least we might somehow be able to find the same sequence close enough as we should jump very close in both cases. If it is not the case it seems problematic as it looks like the 2**128 skipped iterations are not true so I cannot safely rely on a jump if we don't know where it jumps.

Actually I found this issue while implementing a jump for mersenne twister on a JavaScript library I am maintaining. I already had a working implementation of the mersenne generator but needed a reference one for the jump to compare my results to it. So I compared it to the jump of numpy. In the set of tests, I use on my generators one of them - based on property based testing (see hypothesis for Python) - checks that: jump(next(generator)) == next(jump(generator)). And here is how I found the problem.

🤔

In this case I'd maybe call mt19937_gen before jumping. As in https://github.com/numpy/numpy/blob/28f9cbc5e5510cd780b9d6a18a8e8bc6b45a13d7/numpy/random/src/mt19937/mt19937.h#L27-L43

Reseting pos to 0 when pos is N means that the final generator after jump will be: 0 + 2**128
While if pos was N-2 it would be: N - 2 + 2**128
While if pos was N-1 it would be: N - 1 + 2**128

I think the best we can do here is to have fidelity with the original jump implementation. Going off in other directions seems unnecessarily risky.

I think I agree, except some ease of mind, I honestly do not know why there would be a use-case to generate the same identical stream after different order of operation? The only thing is could be to document it, but since this is the reference implementation...

I totally agree that if the problem comes from the reference implementation then it is another issue.

What I was saying is more that given the documentation of this method, as a user I expect that calling jump will offset the generator by 2**128 iterations which should be the same as calling next 2**128 times in a row. Saying that next(jump) !== jump(next) somehow points that this is not the case 🤔

@bashtage maybe we can just tweak the docs a little to be clear that it only jumps about 2**128 states, 2**128 states from the 0 offset postion (whatever it is)?

People still need a mental model of what's going on. From what you've added, I'd still expect next(jumped()) == jumped(next()) to hold true. If that's not true, then what are the guarantees that we can rely on to reason about our programs?

My guess is that the jump is for the whole state and ignores the offset for the next item. @bashtage Is that the case? If so, it should be possible to work around that.

I'm going to try and contact the author to see if he has any notes to share

@bashtage Did you succeed into getting in touch with the author of the paper? I'd really be interested to know the reason why jump seems to jump at very different locations based on how many times we offset the rng before jumping 🤔

As a reminder my original issue was:

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: 2382575482

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: 2001582899

It seems that the current version of numpy (1.19.4) still has a similiar behaviour. Which seems unexpected given what jump should do 🤔

@bashtage Just bringing up the point of dubzzz again: Did you hear back anything of the original author of the jump code paper?

I recently revisited this issue on pure-rand side with the help of Claude. It is still a work-in-progress but it somehow inspired itself from the numpy implementation but made it (in JavaScript) with the jump>next is equivalent to next>jump. I'll will have to re-read it carefully to get why it suddenly worked but given the tests being executed it seems that it made it through.

PR here.

I'm gonna keep you posted when I understand the diff with its implementation and mine (which was translated from numpy)

Makes sense and is easily readable! Do you think it makes sense to refer to the commit (5055c43) or a maybe rather version number (1.18.1 ?) of randomgen, too?

I suspect that the the title should be in quotes, but I don't know if that makes a difference. I'm not familiar with the citation format we use. @mattip @rkern thoughts?

This may be related to regexp for renaming references in numpydoc... looking into it now

hmm doesn't seem like it, at least not a problem with the section I looked at:

https://github.com/numpy/numpydoc/blob/b4c5fd17e2b85c2416a5e586933eee353b58bf7c/numpydoc/numpydoc.py#L43-L67