Add Lcg128CmDxsm64 generator compatible with NumPy's PCG64DXSM

[adamreichold]

Fixes #1200

[adamreichold]

I see at least two issues here:

• I am not sure if LLVM is clever to enough to make use of the fact that the multiplier is only 64 bits in the step function, so this implementation probably leaves performance on the table. Is there a useful u128xu64 primitive that I should use?
• While I did test compatibility with pcg-cpp's pcg_engines::cm_setseq_dxsm_128_64, I am as of yet unable to prove compatibility with NumPy as I do not have direct access to the state and stream initialization parameters which appear to be "hidden" behind NumPy's SeedSequence mechanism.

[adamreichold]

I am as of yet unable to prove compatibility with NumPy

It does work by explicitly setting the state to what results from our (and pcg-cpp's) initialization, but of course this not so helpful for practical compatibility between programs using NumPy and ones using this crate.

In [1]: import numpy.random as npr

In [2]: rng = npr.PCG64DXSM()

In [3]: rng.state
Out[3]: 
{'bit_generator': 'PCG64DXSM',
 'state': {'state': 265141101324562655654025342627628894790,
  'inc': 109099029642168461922701691248838983415},
 'has_uint32': 0,
 'uinteger': 0}

In [4]: s = rng.state

In [5]: s["state"]["state"] = 2378287639543667446576

In [6]: s["state"]["inc"] = 109

In [7]: rng.state = s

In [8]: rng.random_raw()
Out[8]: 17331114245835578256

In [9]: rng.random_raw()
Out[9]: 10267467544499227306

In [10]: rng.random_raw()
Out[10]: 9726600296081716989

In [11]: rng.random_raw()
Out[11]: 10165951391103677450

In [12]: rng.random_raw()
Out[12]: 12131334649314727261

In [13]: rng.random_raw()
Out[13]: 10134094537930450875

[adamreichold]

I am not sure if LLVM is clever to enough to make use of the fact that the multiplier is only 64 bits in the step function, so this implementation probably leaves performance on the table.

It does look like LLVM exploits the cheap multiplier constant: https://godbolt.org/z/q6xcoaMKr

[vks]

@adamreichold Thanks! Does this PR supersede #1201, or do you think we need the other variants as well?

[vks]

Looks good, I only have minor comments, see above.

It would be nice to have a constructor compatible with NumPy's SeedSequence, but this is very optional and of course not required for this PR. (I don't even know whether anyone needs this compatibility.)

[adamreichold]

Thank you for taking the time to review this!

Does this PR supersede #1201, or do you think we need the other variants as well?

Consider the aim for a minimal implementation for rand_pcg, I think this should supersede #1201: This would make a DXSM-based generator available while making minimal code changes which are also recuperated by achieving NumPy compatibility. I just wanted to make both options available for discussion.

Concerning your inline comments: Since I started by copying the pcg128.rs file, those observations would also apply there. Would prefer if I make the same changes in that file as well or should I limit them to the newly created pcg128cm.rs file?

[adamreichold]

It would be nice to have a constructor compatible with NumPy's SeedSequence, but this is very optional and of course not required for this PR. (I don't even know whether anyone needs this compatibility.)

I will look into this. Due to the value stability guarantees given by NumPy, this should be possible without becoming a maintenance nightmare. But indeed, I would like to make this a separate PR to not block these changes on it (as personally, I would not immediately use that compatibility) and make them as limited as possible.

[vks]

Thanks, looks good to me!

@dhardy I think this is a reasonable addition, and minimal enough given that it is limited to one RNG.

[vks]

@adamreichold

Concerning your inline comments: Since I started by copying the pcg128.rs file, those observations would also apply there. Would prefer if I make the same changes in that file as well or should I limit them to the newly created pcg128cm.rs file?

Sure, feel free to update them as well (in a separate commit if possible).

[adamreichold]

Sure, feel free to update them as well (in a separate commit if possible).

Added a separate commit to do the copy-editing in the existing code.

Also note that according to the original author, this variant is supposed to become the "new" pcg64 eventually which seems to be another good reason to add support for this particular variant IMHO: numpy/numpy#13635 (comment)

[vks]

@dhardy Do you still have comments on this PR?

[dhardy]

No, I don't see much to comment on. Adding benchmarks would be nice but not required. @vks already reviewed so I'll let him merge this.

[adamreichold]

Adding benchmarks would be nice but not required.

I added Pcg64Dxsm variants to the generator benchmarks. The results seem to be on par with Pcg64:

test gen_bytes_pcg32         ... bench:     350,388 ns/iter (+/- 9,087) = 2922 MB/s
test gen_bytes_pcg64         ... bench:     221,418 ns/iter (+/- 4,586) = 4624 MB/s
test gen_bytes_pcg64dxsm     ... bench:     224,971 ns/iter (+/- 3,360) = 4551 MB/s
test gen_bytes_pcg64mcg      ... bench:     167,018 ns/iter (+/- 2,175) = 6131 MB/s
test gen_u32_pcg32           ... bench:       1,200 ns/iter (+/- 13) = 3333 MB/s
test gen_u32_pcg64           ... bench:       1,603 ns/iter (+/- 31) = 2495 MB/s
test gen_u32_pcg64dxsm       ... bench:       1,588 ns/iter (+/- 22) = 2518 MB/s
test gen_u32_pcg64mcg        ... bench:       1,165 ns/iter (+/- 13) = 3433 MB/s
test gen_u64_pcg32           ... bench:       2,316 ns/iter (+/- 40) = 3454 MB/s
test gen_u64_pcg64           ... bench:       1,601 ns/iter (+/- 31) = 4996 MB/s
test gen_u64_pcg64dxsm       ... bench:       1,590 ns/iter (+/- 37) = 5031 MB/s
test gen_u64_pcg64mcg        ... bench:       1,165 ns/iter (+/- 411) = 6866 MB/s
test init_pcg32              ... bench:           7 ns/iter (+/- 0)
test init_pcg64              ... bench:          21 ns/iter (+/- 0)
test init_pcg64dxsm          ... bench:          21 ns/iter (+/- 0)
test init_pcg64mcg           ... bench:           7 ns/iter (+/- 0)

[vks]

Thanks!