From a26be10ff1d1a4f64fb23ead017e4eff5b8f4e83 Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Wed, 6 Oct 2021 09:45:21 +0100 Subject: [PATCH] BUG: Correct incorrect advance in PCG with emulated int128 Correct incorrect implemetation of carry in PCG64 and PCG64DXSM when advancing more than 2**64 steps closes #20048 --- doc/release/upcoming_changes/20049.change.rst | 5 +++++ numpy/random/src/pcg64/pcg64.c | 3 +-- numpy/random/tests/test_direct.py | 22 +++++++++++++++++++ 3 files changed, 28 insertions(+), 2 deletions(-) create mode 100644 doc/release/upcoming_changes/20049.change.rst diff --git a/doc/release/upcoming_changes/20049.change.rst b/doc/release/upcoming_changes/20049.change.rst new file mode 100644 index 000000000000..e1f08b3437c8 --- /dev/null +++ b/doc/release/upcoming_changes/20049.change.rst @@ -0,0 +1,5 @@ +Corrected ``advance`` in ``PCG64DSXM`` and ``PCG64`` +---------------------------------------------------- +Fixed a bug in the ``advance`` method of ``PCG64DSXM`` and ``PCG64``. The bug only +affects results when the step was larger than :math:`2^{64}` on platforms +that do not support 128-bit integers(e.g., Windows and 32-bit Linux). diff --git a/numpy/random/src/pcg64/pcg64.c b/numpy/random/src/pcg64/pcg64.c index c623c809b02e..b9be1e39d350 100644 --- a/numpy/random/src/pcg64/pcg64.c +++ b/numpy/random/src/pcg64/pcg64.c @@ -109,8 +109,7 @@ pcg128_t pcg_advance_lcg_128(pcg128_t state, pcg128_t delta, pcg128_t cur_mult, cur_plus = pcg128_mult(pcg128_add(cur_mult, PCG_128BIT_CONSTANT(0u, 1u)), cur_plus); cur_mult = pcg128_mult(cur_mult, cur_mult); - delta.low >>= 1; - delta.low += delta.high & 1; + delta.low = (delta.low >> 1) | (delta.high << 63); delta.high >>= 1; } return pcg128_add(pcg128_mult(acc_mult, state), acc_plus); diff --git a/numpy/random/tests/test_direct.py b/numpy/random/tests/test_direct.py index 29054b70b95a..ea1ebacb63c8 100644 --- a/numpy/random/tests/test_direct.py +++ b/numpy/random/tests/test_direct.py @@ -358,6 +358,17 @@ def test_advance_symmetry(self): assert val_neg == val_pos assert val_big == val_pos + def test_advange_large(self): + rs = Generator(self.bit_generator(38219308213743)) + pcg = rs.bit_generator + state = pcg.state["state"] + initial_state = 287608843259529770491897792873167516365 + assert state["state"] == initial_state + pcg.advance(sum(2**i for i in (96, 64, 32, 16, 8, 4, 2, 1))) + state = pcg.state["state"] + advanced_state = 135275564607035429730177404003164635391 + assert state["state"] == advanced_state + class TestPCG64DXSM(Base): @classmethod @@ -386,6 +397,17 @@ def test_advance_symmetry(self): assert val_neg == val_pos assert val_big == val_pos + def test_advange_large(self): + rs = Generator(self.bit_generator(38219308213743)) + pcg = rs.bit_generator + state = pcg.state + initial_state = 287608843259529770491897792873167516365 + assert state["state"]["state"] == initial_state + pcg.advance(sum(2**i for i in (96, 64, 32, 16, 8, 4, 2, 1))) + state = pcg.state["state"] + advanced_state = 277778083536782149546677086420637664879 + assert state["state"] == advanced_state + class TestMT19937(Base): @classmethod