@rkern probably knows this best. I guess you are probably right, but doubt there is much you can do about it (except trying to document details). I don't know if e.g. pow has enough round off difference for any problems to be possible (even then I guess they are extremely unlikely) at least on IEEE conform implementations? There are some others such as multivariate normal which also depend on the linear algebra solutions (though I doubt that can cause changes in the stream later on).

The sequences of floats, integers, etc. produced by the base rng should be the same assuming that division conforms to the ieee standard for ieee doubles. Platform variations in library functions such as log may change results for more complicated distributions, depending on the particular function involved, which in turn will affect the following sequence of generated results.

I think you are correct that we cannot make quarantees cross platform, or even between library/compiler versions, when certain of the distributions are involved.

To be a bit more precise, roundoff errors and small function differences will only make small differences in floating results if they are not involved in rejection methods. Differences in functions involved in rejection methods can change the whole following sequence, but that is unlikely unless a very large number of trials are made. We do have tests for repeatability and I'm not aware of any failures at this time.

There's one other source of cross-platform wonkiness that is somewhat easier to realize: internal integer values are C longs. With certain extreme parameter values, 32-bit longs overflow and give different results than on platforms with 64-bit longs. Of course, since they are overflowing, the result is crap anyways, but they can change the stream afterwards for all of the non-crap results that follow.

I don't know if e.g. pow has enough round off difference for any problems to be possible (even then I guess they are extremely unlikely)

Possible, definitely: a single ulp difference in a pow result could be enough to cause a rejection-method inner loop to go for one more iteration, and such single ulp differences are very common: there's wide variation in pow implementations across platforms. In contrast to pow, reproducibility of arithmetic operations on IEEE 754-implementing hardware should be a fairly safe assumption, provided well-known issues associated with x87 usage are avoided, and I'd expect basic libm functions like exp, cos, log, etc. to be mostly correctly rounded most of the time.

I agree that you'd have to be very unlucky to run into any of these cases by accident, so that in practice the guarantee is likely to hold almost always. So perhaps my objection can be dismissed under "practicality beats purity" rules.

Yeah, there are two faces to the compatibility policy: what guarantees users can expect and what goals developers should strive towards. Being compatible "up to numerical rounding errors" across platforms is really an aspirational goal of the latter kind, not a firm guarantee of the former. I think I had the developer-goal in my head when I approved the language. And I was probably thinking of the transformation-type algorithms instead of the rejection-type ones.

The multivariate case mentioned by @seberg is probably the most likely to cause problems as the factorization is only unique up to the sign of the eigenvectors if the singular values are distinct and even less unique if the covariance has degenerate eigenvalues, so differences in implementation of the svd that is used can have a big effect. I believe the Cholesky factorization would be more specific.

Can we close this or do we need to clarify the documentation?

Given that NEP 19 has been accepted, it looks like this issue can be closed.