Import python venv for stability

2026-02-15 21:24:16 -08:00
parent 1343e93a59
commit 7d784705c9
4997 changed files with 1628270 additions and 0 deletions
@@ -0,0 +1,595 @@
+import os
+import sys
+from os.path import join
+
+import pytest
+
+import numpy as np
+from numpy.random import (
+    MT19937,
+    PCG64,
+    PCG64DXSM,
+    SFC64,
+    Generator,
+    Philox,
+    RandomState,
+    SeedSequence,
+    default_rng,
+)
+from numpy.random._common import interface
+from numpy.testing import (
+    assert_allclose,
+    assert_array_equal,
+    assert_equal,
+    assert_raises,
+)
+
+try:
+    import cffi  # noqa: F401
+
+    MISSING_CFFI = False
+except ImportError:
+    MISSING_CFFI = True
+
+try:
+    import ctypes  # noqa: F401
+
+    MISSING_CTYPES = False
+except ImportError:
+    MISSING_CTYPES = False
+
+if sys.flags.optimize > 1:
+    # no docstrings present to inspect when PYTHONOPTIMIZE/Py_OptimizeFlag > 1
+    # cffi cannot succeed
+    MISSING_CFFI = True
+
+
+pwd = os.path.dirname(os.path.abspath(__file__))
+
+
+def assert_state_equal(actual, target):
+    for key in actual:
+        if isinstance(actual[key], dict):
+            assert_state_equal(actual[key], target[key])
+        elif isinstance(actual[key], np.ndarray):
+            assert_array_equal(actual[key], target[key])
+        else:
+            assert actual[key] == target[key]
+
+
+def uint32_to_float32(u):
+    return ((u >> np.uint32(8)) * (1.0 / 2**24)).astype(np.float32)
+
+
+def uniform32_from_uint64(x):
+    x = np.uint64(x)
+    upper = np.array(x >> np.uint64(32), dtype=np.uint32)
+    lower = np.uint64(0xffffffff)
+    lower = np.array(x & lower, dtype=np.uint32)
+    joined = np.column_stack([lower, upper]).ravel()
+    return uint32_to_float32(joined)
+
+
+def uniform32_from_uint53(x):
+    x = np.uint64(x) >> np.uint64(16)
+    x = np.uint32(x & np.uint64(0xffffffff))
+    return uint32_to_float32(x)
+
+
+def uniform32_from_uint32(x):
+    return uint32_to_float32(x)
+
+
+def uniform32_from_uint(x, bits):
+    if bits == 64:
+        return uniform32_from_uint64(x)
+    elif bits == 53:
+        return uniform32_from_uint53(x)
+    elif bits == 32:
+        return uniform32_from_uint32(x)
+    else:
+        raise NotImplementedError
+
+
+def uniform_from_uint(x, bits):
+    if bits in (64, 63, 53):
+        return uniform_from_uint64(x)
+    elif bits == 32:
+        return uniform_from_uint32(x)
+
+
+def uniform_from_uint64(x):
+    return (x >> np.uint64(11)) * (1.0 / 9007199254740992.0)
+
+
+def uniform_from_uint32(x):
+    out = np.empty(len(x) // 2)
+    for i in range(0, len(x), 2):
+        a = x[i] >> 5
+        b = x[i + 1] >> 6
+        out[i // 2] = (a * 67108864.0 + b) / 9007199254740992.0
+    return out
+
+
+def uniform_from_dsfmt(x):
+    return x.view(np.double) - 1.0
+
+
+def gauss_from_uint(x, n, bits):
+    if bits in (64, 63):
+        doubles = uniform_from_uint64(x)
+    elif bits == 32:
+        doubles = uniform_from_uint32(x)
+    else:  # bits == 'dsfmt'
+        doubles = uniform_from_dsfmt(x)
+    gauss = []
+    loc = 0
+    x1 = x2 = 0.0
+    while len(gauss) < n:
+        r2 = 2
+        while r2 >= 1.0 or r2 == 0.0:
+            x1 = 2.0 * doubles[loc] - 1.0
+            x2 = 2.0 * doubles[loc + 1] - 1.0
+            r2 = x1 * x1 + x2 * x2
+            loc += 2
+
+        f = np.sqrt(-2.0 * np.log(r2) / r2)
+        gauss.append(f * x2)
+        gauss.append(f * x1)
+
+    return gauss[:n]
+
+
+def test_seedsequence():
+    from numpy.random.bit_generator import (
+        ISeedSequence,
+        ISpawnableSeedSequence,
+        SeedlessSeedSequence,
+    )
+
+    s1 = SeedSequence(range(10), spawn_key=(1, 2), pool_size=6)
+    s1.spawn(10)
+    s2 = SeedSequence(**s1.state)
+    assert_equal(s1.state, s2.state)
+    assert_equal(s1.n_children_spawned, s2.n_children_spawned)
+
+    # The interfaces cannot be instantiated themselves.
+    assert_raises(TypeError, ISeedSequence)
+    assert_raises(TypeError, ISpawnableSeedSequence)
+    dummy = SeedlessSeedSequence()
+    assert_raises(NotImplementedError, dummy.generate_state, 10)
+    assert len(dummy.spawn(10)) == 10
+
+
+def test_generator_spawning():
+    """ Test spawning new generators and bit_generators directly.
+    """
+    rng = np.random.default_rng()
+    seq = rng.bit_generator.seed_seq
+    new_ss = seq.spawn(5)
+    expected_keys = [seq.spawn_key + (i,) for i in range(5)]
+    assert [c.spawn_key for c in new_ss] == expected_keys
+
+    new_bgs = rng.bit_generator.spawn(5)
+    expected_keys = [seq.spawn_key + (i,) for i in range(5, 10)]
+    assert [bg.seed_seq.spawn_key for bg in new_bgs] == expected_keys
+
+    new_rngs = rng.spawn(5)
+    expected_keys = [seq.spawn_key + (i,) for i in range(10, 15)]
+    found_keys = [rng.bit_generator.seed_seq.spawn_key for rng in new_rngs]
+    assert found_keys == expected_keys
+
+    # Sanity check that streams are actually different:
+    assert new_rngs[0].uniform() != new_rngs[1].uniform()
+
+
+def test_non_spawnable():
+    from numpy.random.bit_generator import ISeedSequence
+
+    class FakeSeedSequence:
+        def generate_state(self, n_words, dtype=np.uint32):
+            return np.zeros(n_words, dtype=dtype)
+
+    ISeedSequence.register(FakeSeedSequence)
+
+    rng = np.random.default_rng(FakeSeedSequence())
+
+    with pytest.raises(TypeError, match="The underlying SeedSequence"):
+        rng.spawn(5)
+
+    with pytest.raises(TypeError, match="The underlying SeedSequence"):
+        rng.bit_generator.spawn(5)
+
+
+class Base:
+    dtype = np.uint64
+    data2 = data1 = {}
+
+    @classmethod
+    def setup_class(cls):
+        cls.bit_generator = PCG64
+        cls.bits = 64
+        cls.dtype = np.uint64
+        cls.seed_error_type = TypeError
+        cls.invalid_init_types = []
+        cls.invalid_init_values = []
+
+    @classmethod
+    def _read_csv(cls, filename):
+        with open(filename) as csv:
+            seed = csv.readline()
+            seed = seed.split(',')
+            seed = [int(s.strip(), 0) for s in seed[1:]]
+            data = []
+            for line in csv:
+                data.append(int(line.split(',')[-1].strip(), 0))
+            return {'seed': seed, 'data': np.array(data, dtype=cls.dtype)}
+
+    def test_raw(self):
+        bit_generator = self.bit_generator(*self.data1['seed'])
+        uints = bit_generator.random_raw(1000)
+        assert_equal(uints, self.data1['data'])
+
+        bit_generator = self.bit_generator(*self.data1['seed'])
+        uints = bit_generator.random_raw()
+        assert_equal(uints, self.data1['data'][0])
+
+        bit_generator = self.bit_generator(*self.data2['seed'])
+        uints = bit_generator.random_raw(1000)
+        assert_equal(uints, self.data2['data'])
+
+    def test_random_raw(self):
+        bit_generator = self.bit_generator(*self.data1['seed'])
+        uints = bit_generator.random_raw(output=False)
+        assert uints is None
+        uints = bit_generator.random_raw(1000, output=False)
+        assert uints is None
+
+    def test_gauss_inv(self):
+        n = 25
+        rs = RandomState(self.bit_generator(*self.data1['seed']))
+        gauss = rs.standard_normal(n)
+        assert_allclose(gauss,
+                        gauss_from_uint(self.data1['data'], n, self.bits))
+
+        rs = RandomState(self.bit_generator(*self.data2['seed']))
+        gauss = rs.standard_normal(25)
+        assert_allclose(gauss,
+                        gauss_from_uint(self.data2['data'], n, self.bits))
+
+    def test_uniform_double(self):
+        rs = Generator(self.bit_generator(*self.data1['seed']))
+        vals = uniform_from_uint(self.data1['data'], self.bits)
+        uniforms = rs.random(len(vals))
+        assert_allclose(uniforms, vals)
+        assert_equal(uniforms.dtype, np.float64)
+
+        rs = Generator(self.bit_generator(*self.data2['seed']))
+        vals = uniform_from_uint(self.data2['data'], self.bits)
+        uniforms = rs.random(len(vals))
+        assert_allclose(uniforms, vals)
+        assert_equal(uniforms.dtype, np.float64)
+
+    def test_uniform_float(self):
+        rs = Generator(self.bit_generator(*self.data1['seed']))
+        vals = uniform32_from_uint(self.data1['data'], self.bits)
+        uniforms = rs.random(len(vals), dtype=np.float32)
+        assert_allclose(uniforms, vals)
+        assert_equal(uniforms.dtype, np.float32)
+
+        rs = Generator(self.bit_generator(*self.data2['seed']))
+        vals = uniform32_from_uint(self.data2['data'], self.bits)
+        uniforms = rs.random(len(vals), dtype=np.float32)
+        assert_allclose(uniforms, vals)
+        assert_equal(uniforms.dtype, np.float32)
+
+    def test_repr(self):
+        rs = Generator(self.bit_generator(*self.data1['seed']))
+        assert 'Generator' in repr(rs)
+        assert f'{id(rs):#x}'.upper().replace('X', 'x') in repr(rs)
+
+    def test_str(self):
+        rs = Generator(self.bit_generator(*self.data1['seed']))
+        assert 'Generator' in str(rs)
+        assert str(self.bit_generator.__name__) in str(rs)
+        assert f'{id(rs):#x}'.upper().replace('X', 'x') not in str(rs)
+
+    def test_pickle(self):
+        import pickle
+
+        bit_generator = self.bit_generator(*self.data1['seed'])
+        state = bit_generator.state
+        bitgen_pkl = pickle.dumps(bit_generator)
+        reloaded = pickle.loads(bitgen_pkl)
+        reloaded_state = reloaded.state
+        assert_array_equal(Generator(bit_generator).standard_normal(1000),
+                           Generator(reloaded).standard_normal(1000))
+        assert bit_generator is not reloaded
+        assert_state_equal(reloaded_state, state)
+
+        ss = SeedSequence(100)
+        aa = pickle.loads(pickle.dumps(ss))
+        assert_equal(ss.state, aa.state)
+
+    def test_pickle_preserves_seed_sequence(self):
+        # GH 26234
+        # Add explicit test that bit generators preserve seed sequences
+        import pickle
+
+        bit_generator = self.bit_generator(*self.data1['seed'])
+        ss = bit_generator.seed_seq
+        bg_plk = pickle.loads(pickle.dumps(bit_generator))
+        ss_plk = bg_plk.seed_seq
+        assert_equal(ss.state, ss_plk.state)
+        assert_equal(ss.pool, ss_plk.pool)
+
+        bit_generator.seed_seq.spawn(10)
+        bg_plk = pickle.loads(pickle.dumps(bit_generator))
+        ss_plk = bg_plk.seed_seq
+        assert_equal(ss.state, ss_plk.state)
+        assert_equal(ss.n_children_spawned, ss_plk.n_children_spawned)
+
+    def test_invalid_state_type(self):
+        bit_generator = self.bit_generator(*self.data1['seed'])
+        with pytest.raises(TypeError):
+            bit_generator.state = {'1'}
+
+    def test_invalid_state_value(self):
+        bit_generator = self.bit_generator(*self.data1['seed'])
+        state = bit_generator.state
+        state['bit_generator'] = 'otherBitGenerator'
+        with pytest.raises(ValueError):
+            bit_generator.state = state
+
+    def test_invalid_init_type(self):
+        bit_generator = self.bit_generator
+        for st in self.invalid_init_types:
+            with pytest.raises(TypeError):
+                bit_generator(*st)
+
+    def test_invalid_init_values(self):
+        bit_generator = self.bit_generator
+        for st in self.invalid_init_values:
+            with pytest.raises((ValueError, OverflowError)):
+                bit_generator(*st)
+
+    def test_benchmark(self):
+        bit_generator = self.bit_generator(*self.data1['seed'])
+        bit_generator._benchmark(1)
+        bit_generator._benchmark(1, 'double')
+        with pytest.raises(ValueError):
+            bit_generator._benchmark(1, 'int32')
+
+    @pytest.mark.skipif(MISSING_CFFI, reason='cffi not available')
+    def test_cffi(self):
+        bit_generator = self.bit_generator(*self.data1['seed'])
+        cffi_interface = bit_generator.cffi
+        assert isinstance(cffi_interface, interface)
+        other_cffi_interface = bit_generator.cffi
+        assert other_cffi_interface is cffi_interface
+
+    @pytest.mark.skipif(MISSING_CTYPES, reason='ctypes not available')
+    def test_ctypes(self):
+        bit_generator = self.bit_generator(*self.data1['seed'])
+        ctypes_interface = bit_generator.ctypes
+        assert isinstance(ctypes_interface, interface)
+        other_ctypes_interface = bit_generator.ctypes
+        assert other_ctypes_interface is ctypes_interface
+
+    def test_getstate(self):
+        bit_generator = self.bit_generator(*self.data1['seed'])
+        state = bit_generator.state
+        alt_state = bit_generator.__getstate__()
+        assert isinstance(alt_state, tuple)
+        assert_state_equal(state, alt_state[0])
+        assert isinstance(alt_state[1], SeedSequence)
+
+class TestPhilox(Base):
+    @classmethod
+    def setup_class(cls):
+        cls.bit_generator = Philox
+        cls.bits = 64
+        cls.dtype = np.uint64
+        cls.data1 = cls._read_csv(
+            join(pwd, './data/philox-testset-1.csv'))
+        cls.data2 = cls._read_csv(
+            join(pwd, './data/philox-testset-2.csv'))
+        cls.seed_error_type = TypeError
+        cls.invalid_init_types = []
+        cls.invalid_init_values = [(1, None, 1), (-1,), (None, None, 2 ** 257 + 1)]
+
+    def test_set_key(self):
+        bit_generator = self.bit_generator(*self.data1['seed'])
+        state = bit_generator.state
+        keyed = self.bit_generator(counter=state['state']['counter'],
+                                   key=state['state']['key'])
+        assert_state_equal(bit_generator.state, keyed.state)
+
+
+class TestPCG64(Base):
+    @classmethod
+    def setup_class(cls):
+        cls.bit_generator = PCG64
+        cls.bits = 64
+        cls.dtype = np.uint64
+        cls.data1 = cls._read_csv(join(pwd, './data/pcg64-testset-1.csv'))
+        cls.data2 = cls._read_csv(join(pwd, './data/pcg64-testset-2.csv'))
+        cls.seed_error_type = (ValueError, TypeError)
+        cls.invalid_init_types = [(3.2,), ([None],), (1, None)]
+        cls.invalid_init_values = [(-1,)]
+
+    def test_advance_symmetry(self):
+        rs = Generator(self.bit_generator(*self.data1['seed']))
+        state = rs.bit_generator.state
+        step = -0x9e3779b97f4a7c150000000000000000
+        rs.bit_generator.advance(step)
+        val_neg = rs.integers(10)
+        rs.bit_generator.state = state
+        rs.bit_generator.advance(2**128 + step)
+        val_pos = rs.integers(10)
+        rs.bit_generator.state = state
+        rs.bit_generator.advance(10 * 2**128 + step)
+        val_big = rs.integers(10)
+        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
+    def setup_class(cls):
+        cls.bit_generator = PCG64DXSM
+        cls.bits = 64
+        cls.dtype = np.uint64
+        cls.data1 = cls._read_csv(join(pwd, './data/pcg64dxsm-testset-1.csv'))
+        cls.data2 = cls._read_csv(join(pwd, './data/pcg64dxsm-testset-2.csv'))
+        cls.seed_error_type = (ValueError, TypeError)
+        cls.invalid_init_types = [(3.2,), ([None],), (1, None)]
+        cls.invalid_init_values = [(-1,)]
+
+    def test_advance_symmetry(self):
+        rs = Generator(self.bit_generator(*self.data1['seed']))
+        state = rs.bit_generator.state
+        step = -0x9e3779b97f4a7c150000000000000000
+        rs.bit_generator.advance(step)
+        val_neg = rs.integers(10)
+        rs.bit_generator.state = state
+        rs.bit_generator.advance(2**128 + step)
+        val_pos = rs.integers(10)
+        rs.bit_generator.state = state
+        rs.bit_generator.advance(10 * 2**128 + step)
+        val_big = rs.integers(10)
+        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
+    def setup_class(cls):
+        cls.bit_generator = MT19937
+        cls.bits = 32
+        cls.dtype = np.uint32
+        cls.data1 = cls._read_csv(join(pwd, './data/mt19937-testset-1.csv'))
+        cls.data2 = cls._read_csv(join(pwd, './data/mt19937-testset-2.csv'))
+        cls.seed_error_type = ValueError
+        cls.invalid_init_types = []
+        cls.invalid_init_values = [(-1,)]
+
+    def test_seed_float_array(self):
+        assert_raises(TypeError, self.bit_generator, np.array([np.pi]))
+        assert_raises(TypeError, self.bit_generator, np.array([-np.pi]))
+        assert_raises(TypeError, self.bit_generator, np.array([np.pi, -np.pi]))
+        assert_raises(TypeError, self.bit_generator, np.array([0, np.pi]))
+        assert_raises(TypeError, self.bit_generator, [np.pi])
+        assert_raises(TypeError, self.bit_generator, [0, np.pi])
+
+    def test_state_tuple(self):
+        rs = Generator(self.bit_generator(*self.data1['seed']))
+        bit_generator = rs.bit_generator
+        state = bit_generator.state
+        desired = rs.integers(2 ** 16)
+        tup = (state['bit_generator'], state['state']['key'],
+               state['state']['pos'])
+        bit_generator.state = tup
+        actual = rs.integers(2 ** 16)
+        assert_equal(actual, desired)
+        tup = tup + (0, 0.0)
+        bit_generator.state = tup
+        actual = rs.integers(2 ** 16)
+        assert_equal(actual, desired)
+
+
+class TestSFC64(Base):
+    @classmethod
+    def setup_class(cls):
+        cls.bit_generator = SFC64
+        cls.bits = 64
+        cls.dtype = np.uint64
+        cls.data1 = cls._read_csv(
+            join(pwd, './data/sfc64-testset-1.csv'))
+        cls.data2 = cls._read_csv(
+            join(pwd, './data/sfc64-testset-2.csv'))
+        cls.seed_error_type = (ValueError, TypeError)
+        cls.invalid_init_types = [(3.2,), ([None],), (1, None)]
+        cls.invalid_init_values = [(-1,)]
+
+    def test_legacy_pickle(self):
+        # Pickling format was changed in 2.0.x
+        import gzip
+        import pickle
+
+        expected_state = np.array(
+            [
+                9957867060933711493,
+                532597980065565856,
+                14769588338631205282,
+                13
+            ],
+            dtype=np.uint64
+        )
+
+        base_path = os.path.split(os.path.abspath(__file__))[0]
+        pkl_file = os.path.join(base_path, "data", "sfc64_np126.pkl.gz")
+        with gzip.open(pkl_file) as gz:
+            sfc = pickle.load(gz)
+
+        assert isinstance(sfc, SFC64)
+        assert_equal(sfc.state["state"]["state"], expected_state)
+
+
+class TestDefaultRNG:
+    def test_seed(self):
+        for args in [(), (None,), (1234,), ([1234, 5678],)]:
+            rg = default_rng(*args)
+            assert isinstance(rg.bit_generator, PCG64)
+
+    def test_passthrough(self):
+        bg = Philox()
+        rg = default_rng(bg)
+        assert rg.bit_generator is bg
+        rg2 = default_rng(rg)
+        assert rg2 is rg
+        assert rg2.bit_generator is bg
+
+    @pytest.mark.thread_unsafe(
+        reason="np.random.set_bit_generator affects global state"
+    )
+    def test_coercion_RandomState_Generator(self):
+        # use default_rng to coerce RandomState to Generator
+        rs = RandomState(1234)
+        rg = default_rng(rs)
+        assert isinstance(rg.bit_generator, MT19937)
+        assert rg.bit_generator is rs._bit_generator
+
+        # RandomState with a non MT19937 bit generator
+        _original = np.random.get_bit_generator()
+        bg = PCG64(12342298)
+        np.random.set_bit_generator(bg)
+        rs = np.random.mtrand._rand
+        rg = default_rng(rs)
+        assert rg.bit_generator is bg
+
+        # vital to get global state back to original, otherwise
+        # other tests start to fail.
+        np.random.set_bit_generator(_original)
@@ -0,0 +1,131 @@
+import os
+import shutil
+import subprocess
+import sys
+import sysconfig
+import warnings
+from importlib.util import module_from_spec, spec_from_file_location
+
+import pytest
+
+import numpy as np
+from numpy.testing import IS_EDITABLE, IS_WASM
+
+try:
+    import cffi
+except ImportError:
+    cffi = None
+
+if sys.flags.optimize > 1:
+    # no docstrings present to inspect when PYTHONOPTIMIZE/Py_OptimizeFlag > 1
+    # cffi cannot succeed
+    cffi = None
+
+try:
+    with warnings.catch_warnings(record=True) as w:
+        # numba issue gh-4733
+        warnings.filterwarnings('always', '', DeprecationWarning)
+        import numba
+except (ImportError, SystemError):
+    # Certain numpy/numba versions trigger a SystemError due to a numba bug
+    numba = None
+
+try:
+    import cython
+    from Cython.Compiler.Version import version as cython_version
+except ImportError:
+    cython = None
+else:
+    from numpy._utils import _pep440
+    # Note: keep in sync with the one in pyproject.toml
+    required_version = '3.0.6'
+    if _pep440.parse(cython_version) < _pep440.Version(required_version):
+        # too old or wrong cython, skip the test
+        cython = None
+
+
+@pytest.mark.skipif(
+    IS_EDITABLE,
+    reason='Editable install cannot find .pxd headers'
+)
+@pytest.mark.skipif(
+        sys.platform == "win32" and sys.maxsize < 2**32,
+        reason="Failing in 32-bit Windows wheel build job, skip for now"
+)
+@pytest.mark.skipif(IS_WASM, reason="Can't start subprocess")
+@pytest.mark.skipif(cython is None, reason="requires cython")
+@pytest.mark.skipif(sysconfig.get_platform() == 'win-arm64',
+                    reason='Meson unable to find MSVC linker on win-arm64')
+@pytest.mark.slow
+@pytest.mark.thread_unsafe(
+    reason="building cython code in a subprocess doesn't make sense to do in many "
+           "threads and sometimes crashes"
+)
+def test_cython(tmp_path):
+    import glob
+    # build the examples in a temporary directory
+    srcdir = os.path.join(os.path.dirname(__file__), '..')
+    shutil.copytree(srcdir, tmp_path / 'random')
+    build_dir = tmp_path / 'random' / '_examples' / 'cython'
+    target_dir = build_dir / "build"
+    os.makedirs(target_dir, exist_ok=True)
+    # Ensure we use the correct Python interpreter even when `meson` is
+    # installed in a different Python environment (see gh-24956)
+    native_file = str(build_dir / 'interpreter-native-file.ini')
+    with open(native_file, 'w') as f:
+        f.write("[binaries]\n")
+        f.write(f"python = '{sys.executable}'\n")
+        f.write(f"python3 = '{sys.executable}'")
+    if sys.platform == "win32":
+        subprocess.check_call(["meson", "setup",
+                               "--buildtype=release",
+                               "--vsenv", "--native-file", native_file,
+                               str(build_dir)],
+                              cwd=target_dir,
+                              )
+    else:
+        subprocess.check_call(["meson", "setup",
+                               "--native-file", native_file, str(build_dir)],
+                              cwd=target_dir
+                              )
+    subprocess.check_call(["meson", "compile", "-vv"], cwd=target_dir)
+
+    # gh-16162: make sure numpy's __init__.pxd was used for cython
+    # not really part of this test, but it is a convenient place to check
+
+    g = glob.glob(str(target_dir / "*" / "extending.pyx.c"))
+    with open(g[0]) as fid:
+        txt_to_find = 'NumPy API declarations from "numpy/__init__'
+        for line in fid:
+            if txt_to_find in line:
+                break
+        else:
+            assert False, f"Could not find '{txt_to_find}' in C file, wrong pxd used"
+    # import without adding the directory to sys.path
+    suffix = sysconfig.get_config_var('EXT_SUFFIX')
+
+    def load(modname):
+        so = (target_dir / modname).with_suffix(suffix)
+        spec = spec_from_file_location(modname, so)
+        mod = module_from_spec(spec)
+        spec.loader.exec_module(mod)
+        return mod
+
+    # test that the module can be imported
+    load("extending")
+    load("extending_cpp")
+    # actually test the cython c-extension
+    extending_distributions = load("extending_distributions")
+    from numpy.random import PCG64
+    values = extending_distributions.uniforms_ex(PCG64(0), 10, 'd')
+    assert values.shape == (10,)
+    assert values.dtype == np.float64
+
+@pytest.mark.skipif(numba is None or cffi is None,
+                    reason="requires numba and cffi")
+def test_numba():
+    from numpy.random._examples.numba import extending  # noqa: F401
+
+@pytest.mark.skipif(cffi is None, reason="requires cffi")
+def test_cffi():
+    from numpy.random._examples.cffi import extending  # noqa: F401
@@ -0,0 +1,221 @@
+import pytest
+
+import numpy as np
+from numpy.random import MT19937, Generator
+from numpy.testing import assert_, assert_array_equal
+
+
+class TestRegression:
+    def _create_generator(self):
+        return Generator(MT19937(121263137472525314065))
+
+    def test_vonmises_range(self):
+        # Make sure generated random variables are in [-pi, pi].
+        # Regression test for ticket #986.
+        mt19937 = self._create_generator()
+        for mu in np.linspace(-7., 7., 5):
+            r = mt19937.vonmises(mu, 1, 50)
+            assert_(np.all(r > -np.pi) and np.all(r <= np.pi))
+
+    def test_hypergeometric_range(self):
+        # Test for ticket #921
+        mt19937 = self._create_generator()
+        assert_(np.all(mt19937.hypergeometric(3, 18, 11, size=10) < 4))
+        assert_(np.all(mt19937.hypergeometric(18, 3, 11, size=10) > 0))
+
+        # Test for ticket #5623
+        args = (2**20 - 2, 2**20 - 2, 2**20 - 2)  # Check for 32-bit systems
+        assert_(mt19937.hypergeometric(*args) > 0)
+
+    def test_logseries_convergence(self):
+        # Test for ticket #923
+        mt19937 = self._create_generator()
+        N = 1000
+        rvsn = mt19937.logseries(0.8, size=N)
+        # these two frequency counts should be close to theoretical
+        # numbers with this large sample
+        # theoretical large N result is 0.49706795
+        freq = np.sum(rvsn == 1) / N
+        msg = f'Frequency was {freq:f}, should be > 0.45'
+        assert_(freq > 0.45, msg)
+        # theoretical large N result is 0.19882718
+        freq = np.sum(rvsn == 2) / N
+        msg = f'Frequency was {freq:f}, should be < 0.23'
+        assert_(freq < 0.23, msg)
+
+    def test_shuffle_mixed_dimension(self):
+        # Test for trac ticket #2074
+        for t in [[1, 2, 3, None],
+                  [(1, 1), (2, 2), (3, 3), None],
+                  [1, (2, 2), (3, 3), None],
+                  [(1, 1), 2, 3, None]]:
+            mt19937 = Generator(MT19937(12345))
+            shuffled = np.array(t, dtype=object)
+            mt19937.shuffle(shuffled)
+            expected = np.array([t[2], t[0], t[3], t[1]], dtype=object)
+            assert_array_equal(np.array(shuffled, dtype=object), expected)
+
+    def test_call_within_randomstate(self):
+        # Check that custom BitGenerator does not call into global state
+        res = np.array([1, 8, 0, 1, 5, 3, 3, 8, 1, 4])
+        for i in range(3):
+            mt19937 = Generator(MT19937(i))
+            m = Generator(MT19937(4321))
+            # If m.state is not honored, the result will change
+            assert_array_equal(m.choice(10, size=10, p=np.ones(10) / 10.), res)
+
+    def test_multivariate_normal_size_types(self):
+        # Test for multivariate_normal issue with 'size' argument.
+        # Check that the multivariate_normal size argument can be a
+        # numpy integer.
+        mt19937 = self._create_generator()
+        mt19937.multivariate_normal([0], [[0]], size=1)
+        mt19937.multivariate_normal([0], [[0]], size=np.int_(1))
+        mt19937.multivariate_normal([0], [[0]], size=np.int64(1))
+
+    def test_beta_small_parameters(self):
+        # Test that beta with small a and b parameters does not produce
+        # NaNs due to roundoff errors causing 0 / 0, gh-5851
+        mt19937 = self._create_generator()
+        x = mt19937.beta(0.0001, 0.0001, size=100)
+        assert_(not np.any(np.isnan(x)), 'Nans in mt19937.beta')
+
+    def test_beta_very_small_parameters(self):
+        # gh-24203: beta would hang with very small parameters.
+        mt19937 = self._create_generator()
+        mt19937.beta(1e-49, 1e-40)
+
+    def test_beta_ridiculously_small_parameters(self):
+        # gh-24266: beta would generate nan when the parameters
+        # were subnormal or a small multiple of the smallest normal.
+        mt19937 = self._create_generator()
+        tiny = np.finfo(1.0).tiny
+        x = mt19937.beta(tiny / 32, tiny / 40, size=50)
+        assert not np.any(np.isnan(x))
+
+    def test_beta_expected_zero_frequency(self):
+        # gh-24475: For small a and b (e.g. a=0.0025, b=0.0025), beta
+        # would generate too many zeros.
+        mt19937 = self._create_generator()
+        a = 0.0025
+        b = 0.0025
+        n = 1000000
+        x = mt19937.beta(a, b, size=n)
+        nzeros = np.count_nonzero(x == 0)
+        # beta CDF at x = np.finfo(np.double).smallest_subnormal/2
+        # is p = 0.0776169083131899, e.g,
+        #
+        #    import numpy as np
+        #    from mpmath import mp
+        #    mp.dps = 160
+        #    x = mp.mpf(np.finfo(np.float64).smallest_subnormal)/2
+        #    # CDF of the beta distribution at x:
+        #    p = mp.betainc(a, b, x1=0, x2=x, regularized=True)
+        #    n = 1000000
+        #    exprected_freq = float(n*p)
+        #
+        expected_freq = 77616.90831318991
+        assert 0.95 * expected_freq < nzeros < 1.05 * expected_freq
+
+    def test_choice_sum_of_probs_tolerance(self):
+        # The sum of probs should be 1.0 with some tolerance.
+        # For low precision dtypes the tolerance was too tight.
+        # See numpy github issue 6123.
+        mt19937 = self._create_generator()
+        a = [1, 2, 3]
+        counts = [4, 4, 2]
+        for dt in np.float16, np.float32, np.float64:
+            probs = np.array(counts, dtype=dt) / sum(counts)
+            c = mt19937.choice(a, p=probs)
+            assert_(c in a)
+            with pytest.raises(ValueError):
+                mt19937.choice(a, p=probs * 0.9)
+
+    def test_shuffle_of_array_of_different_length_strings(self):
+        # Test that permuting an array of different length strings
+        # will not cause a segfault on garbage collection
+        # Tests gh-7710
+        mt19937 = self._create_generator()
+
+        a = np.array(['a', 'a' * 1000])
+
+        for _ in range(100):
+            mt19937.shuffle(a)
+
+        # Force Garbage Collection - should not segfault.
+        import gc
+        gc.collect()
+
+    def test_shuffle_of_array_of_objects(self):
+        # Test that permuting an array of objects will not cause
+        # a segfault on garbage collection.
+        # See gh-7719
+        mt19937 = self._create_generator()
+        a = np.array([np.arange(1), np.arange(4)], dtype=object)
+
+        for _ in range(1000):
+            mt19937.shuffle(a)
+
+        # Force Garbage Collection - should not segfault.
+        import gc
+        gc.collect()
+
+    def test_permutation_subclass(self):
+
+        class N(np.ndarray):
+            pass
+
+        mt19937 = Generator(MT19937(1))
+        orig = np.arange(3).view(N)
+        perm = mt19937.permutation(orig)
+        assert_array_equal(perm, np.array([2, 0, 1]))
+        assert_array_equal(orig, np.arange(3).view(N))
+
+        class M:
+            a = np.arange(5)
+
+            def __array__(self, dtype=None, copy=None):
+                return self.a
+
+        mt19937 = Generator(MT19937(1))
+        m = M()
+        perm = mt19937.permutation(m)
+        assert_array_equal(perm, np.array([4, 1, 3, 0, 2]))
+        assert_array_equal(m.__array__(), np.arange(5))
+
+    def test_gamma_0(self):
+        mt19937 = self._create_generator()
+        assert mt19937.standard_gamma(0.0) == 0.0
+        assert_array_equal(mt19937.standard_gamma([0.0]), 0.0)
+
+        actual = mt19937.standard_gamma([0.0], dtype='float')
+        expected = np.array([0.], dtype=np.float32)
+        assert_array_equal(actual, expected)
+
+    def test_geometric_tiny_prob(self):
+        # Regression test for gh-17007.
+        # When p = 1e-30, the probability that a sample will exceed 2**63-1
+        # is 0.9999999999907766, so we expect the result to be all 2**63-1.
+        mt19937 = self._create_generator()
+        assert_array_equal(mt19937.geometric(p=1e-30, size=3),
+                           np.iinfo(np.int64).max)
+
+    def test_zipf_large_parameter(self):
+        # Regression test for part of gh-9829: a call such as rng.zipf(10000)
+        # would hang.
+        mt19937 = self._create_generator()
+        n = 8
+        sample = mt19937.zipf(10000, size=n)
+        assert_array_equal(sample, np.ones(n, dtype=np.int64))
+
+    def test_zipf_a_near_1(self):
+        # Regression test for gh-9829: a call such as rng.zipf(1.0000000000001)
+        # would hang.
+        mt19937 = self._create_generator()
+        n = 100000
+        sample = mt19937.zipf(1.0000000000001, size=n)
+        # Not much of a test, but let's do something more than verify that
+        # it doesn't hang.  Certainly for a monotonically decreasing
+        # discrete distribution truncated to signed 64 bit integers, more
+        # than half should be less than 2**62.
+        assert np.count_nonzero(sample < 2**62) > n / 2
@@ -0,0 +1,213 @@
+import sys
+
+import pytest
+
+import numpy as np
+from numpy import random
+from numpy.testing import assert_, assert_array_equal, assert_raises
+
+
+class TestRegression:
+
+    def test_VonMises_range(self):
+        # Make sure generated random variables are in [-pi, pi].
+        # Regression test for ticket #986.
+        for mu in np.linspace(-7., 7., 5):
+            r = random.vonmises(mu, 1, 50)
+            assert_(np.all(r > -np.pi) and np.all(r <= np.pi))
+
+    def test_hypergeometric_range(self):
+        # Test for ticket #921
+        assert_(np.all(random.hypergeometric(3, 18, 11, size=10) < 4))
+        assert_(np.all(random.hypergeometric(18, 3, 11, size=10) > 0))
+
+        # Test for ticket #5623
+        args = [
+            (2**20 - 2, 2**20 - 2, 2**20 - 2),  # Check for 32-bit systems
+        ]
+        is_64bits = sys.maxsize > 2**32
+        if is_64bits and sys.platform != 'win32':
+            # Check for 64-bit systems
+            args.append((2**40 - 2, 2**40 - 2, 2**40 - 2))
+        for arg in args:
+            assert_(random.hypergeometric(*arg) > 0)
+
+    def test_logseries_convergence(self):
+        # Test for ticket #923
+        N = 1000
+        random.seed(0)
+        rvsn = random.logseries(0.8, size=N)
+        # these two frequency counts should be close to theoretical
+        # numbers with this large sample
+        # theoretical large N result is 0.49706795
+        freq = np.sum(rvsn == 1) / N
+        msg = f'Frequency was {freq:f}, should be > 0.45'
+        assert_(freq > 0.45, msg)
+        # theoretical large N result is 0.19882718
+        freq = np.sum(rvsn == 2) / N
+        msg = f'Frequency was {freq:f}, should be < 0.23'
+        assert_(freq < 0.23, msg)
+
+    def test_shuffle_mixed_dimension(self):
+        # Test for trac ticket #2074
+        for t in [[1, 2, 3, None],
+                  [(1, 1), (2, 2), (3, 3), None],
+                  [1, (2, 2), (3, 3), None],
+                  [(1, 1), 2, 3, None]]:
+            rng = random.RandomState(12345)
+            shuffled = list(t)
+            rng.shuffle(shuffled)
+            expected = np.array([t[0], t[3], t[1], t[2]], dtype=object)
+            assert_array_equal(np.array(shuffled, dtype=object), expected)
+
+    def test_call_within_randomstate(self):
+        # Check that custom RandomState does not call into global state
+        m = random.RandomState()
+        res = np.array([0, 8, 7, 2, 1, 9, 4, 7, 0, 3])
+        for i in range(3):
+            random.seed(i)
+            m.seed(4321)
+            # If m.state is not honored, the result will change
+            assert_array_equal(m.choice(10, size=10, p=np.ones(10) / 10.), res)
+
+    def test_multivariate_normal_size_types(self):
+        # Test for multivariate_normal issue with 'size' argument.
+        # Check that the multivariate_normal size argument can be a
+        # numpy integer.
+        random.multivariate_normal([0], [[0]], size=1)
+        random.multivariate_normal([0], [[0]], size=np.int_(1))
+        random.multivariate_normal([0], [[0]], size=np.int64(1))
+
+    def test_beta_small_parameters(self):
+        # Test that beta with small a and b parameters does not produce
+        # NaNs due to roundoff errors causing 0 / 0, gh-5851
+        random.seed(1234567890)
+        x = random.beta(0.0001, 0.0001, size=100)
+        assert_(not np.any(np.isnan(x)), 'Nans in random.beta')
+
+    def test_choice_sum_of_probs_tolerance(self):
+        # The sum of probs should be 1.0 with some tolerance.
+        # For low precision dtypes the tolerance was too tight.
+        # See numpy github issue 6123.
+        random.seed(1234)
+        a = [1, 2, 3]
+        counts = [4, 4, 2]
+        for dt in np.float16, np.float32, np.float64:
+            probs = np.array(counts, dtype=dt) / sum(counts)
+            c = random.choice(a, p=probs)
+            assert_(c in a)
+            assert_raises(ValueError, random.choice, a, p=probs * 0.9)
+
+    def test_shuffle_of_array_of_different_length_strings(self):
+        # Test that permuting an array of different length strings
+        # will not cause a segfault on garbage collection
+        # Tests gh-7710
+        random.seed(1234)
+
+        a = np.array(['a', 'a' * 1000])
+
+        for _ in range(100):
+            random.shuffle(a)
+
+        # Force Garbage Collection - should not segfault.
+        import gc
+        gc.collect()
+
+    def test_shuffle_of_array_of_objects(self):
+        # Test that permuting an array of objects will not cause
+        # a segfault on garbage collection.
+        # See gh-7719
+        random.seed(1234)
+        a = np.array([np.arange(1), np.arange(4)], dtype=object)
+
+        for _ in range(1000):
+            random.shuffle(a)
+
+        # Force Garbage Collection - should not segfault.
+        import gc
+        gc.collect()
+
+    def test_permutation_subclass(self):
+        class N(np.ndarray):
+            pass
+
+        rng = random.RandomState(1)
+        orig = np.arange(3).view(N)
+        perm = rng.permutation(orig)
+        assert_array_equal(perm, np.array([0, 2, 1]))
+        assert_array_equal(orig, np.arange(3).view(N))
+
+        class M:
+            a = np.arange(5)
+
+            def __array__(self, dtype=None, copy=None):
+                return self.a
+
+        rng = random.RandomState(1)
+        m = M()
+        perm = rng.permutation(m)
+        assert_array_equal(perm, np.array([2, 1, 4, 0, 3]))
+        assert_array_equal(m.__array__(), np.arange(5))
+
+    def test_warns_byteorder(self):
+        # GH 13159
+        other_byteord_dt = '<i4' if sys.byteorder == 'big' else '>i4'
+        with pytest.deprecated_call(match='non-native byteorder is not'):
+            random.randint(0, 200, size=10, dtype=other_byteord_dt)
+
+    def test_named_argument_initialization(self):
+        # GH 13669
+        rs1 = np.random.RandomState(123456789)
+        rs2 = np.random.RandomState(seed=123456789)
+        assert rs1.randint(0, 100) == rs2.randint(0, 100)
+
+    def test_choice_retun_dtype(self):
+        # GH 9867, now long since the NumPy default changed.
+        c = np.random.choice(10, p=[.1] * 10, size=2)
+        assert c.dtype == np.dtype(np.long)
+        c = np.random.choice(10, p=[.1] * 10, replace=False, size=2)
+        assert c.dtype == np.dtype(np.long)
+        c = np.random.choice(10, size=2)
+        assert c.dtype == np.dtype(np.long)
+        c = np.random.choice(10, replace=False, size=2)
+        assert c.dtype == np.dtype(np.long)
+
+    @pytest.mark.skipif(np.iinfo('l').max < 2**32,
+                        reason='Cannot test with 32-bit C long')
+    def test_randint_117(self):
+        # GH 14189
+        rng = random.RandomState(0)
+        expected = np.array([2357136044, 2546248239, 3071714933, 3626093760,
+                             2588848963, 3684848379, 2340255427, 3638918503,
+                             1819583497, 2678185683], dtype='int64')
+        actual = rng.randint(2**32, size=10)
+        assert_array_equal(actual, expected)
+
+    def test_p_zero_stream(self):
+        # Regression test for gh-14522.  Ensure that future versions
+        # generate the same variates as version 1.16.
+        rng = random.RandomState(12345)
+        assert_array_equal(rng.binomial(1, [0, 0.25, 0.5, 0.75, 1]),
+                           [0, 0, 0, 1, 1])
+
+    def test_n_zero_stream(self):
+        # Regression test for gh-14522.  Ensure that future versions
+        # generate the same variates as version 1.16.
+        rng = random.RandomState(8675309)
+        expected = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+                             [3, 4, 2, 3, 3, 1, 5, 3, 1, 3]])
+        assert_array_equal(rng.binomial([[0], [10]], 0.25, size=(2, 10)),
+                           expected)
+
+
+def test_multinomial_empty():
+    # gh-20483
+    # Ensure that empty p-vals are correctly handled
+    assert random.multinomial(10, []).shape == (0,)
+    assert random.multinomial(3, [], size=(7, 5, 3)).shape == (7, 5, 3, 0)
+
+
+def test_multinomial_1d_pval():
+    # gh-20483
+    with pytest.raises(TypeError, match="pvals must be a 1-d"):
+        random.multinomial(10, 0.3)
@@ -0,0 +1,175 @@
+import inspect
+import sys
+
+import pytest
+
+import numpy as np
+from numpy import random
+from numpy.testing import IS_PYPY, assert_, assert_array_equal, assert_raises
+
+
+class TestRegression:
+
+    def test_VonMises_range(self):
+        # Make sure generated random variables are in [-pi, pi].
+        # Regression test for ticket #986.
+        for mu in np.linspace(-7., 7., 5):
+            r = random.mtrand.vonmises(mu, 1, 50)
+            assert_(np.all(r > -np.pi) and np.all(r <= np.pi))
+
+    def test_hypergeometric_range(self):
+        # Test for ticket #921
+        assert_(np.all(np.random.hypergeometric(3, 18, 11, size=10) < 4))
+        assert_(np.all(np.random.hypergeometric(18, 3, 11, size=10) > 0))
+
+        # Test for ticket #5623
+        args = [
+            (2**20 - 2, 2**20 - 2, 2**20 - 2),  # Check for 32-bit systems
+        ]
+        is_64bits = sys.maxsize > 2**32
+        if is_64bits and sys.platform != 'win32':
+            # Check for 64-bit systems
+            args.append((2**40 - 2, 2**40 - 2, 2**40 - 2))
+        for arg in args:
+            assert_(np.random.hypergeometric(*arg) > 0)
+
+    def test_logseries_convergence(self):
+        # Test for ticket #923
+        N = 1000
+        np.random.seed(0)
+        rvsn = np.random.logseries(0.8, size=N)
+        # these two frequency counts should be close to theoretical
+        # numbers with this large sample
+        # theoretical large N result is 0.49706795
+        freq = np.sum(rvsn == 1) / N
+        msg = f'Frequency was {freq:f}, should be > 0.45'
+        assert_(freq > 0.45, msg)
+        # theoretical large N result is 0.19882718
+        freq = np.sum(rvsn == 2) / N
+        msg = f'Frequency was {freq:f}, should be < 0.23'
+        assert_(freq < 0.23, msg)
+
+    def test_shuffle_mixed_dimension(self):
+        # Test for trac ticket #2074
+        for t in [[1, 2, 3, None],
+                  [(1, 1), (2, 2), (3, 3), None],
+                  [1, (2, 2), (3, 3), None],
+                  [(1, 1), 2, 3, None]]:
+            rng = np.random.RandomState(12345)
+            shuffled = list(t)
+            rng.shuffle(shuffled)
+            expected = np.array([t[0], t[3], t[1], t[2]], dtype=object)
+            assert_array_equal(np.array(shuffled, dtype=object), expected)
+
+    def test_call_within_randomstate(self):
+        # Check that custom RandomState does not call into global state
+        m = np.random.RandomState()
+        res = np.array([0, 8, 7, 2, 1, 9, 4, 7, 0, 3])
+        for i in range(3):
+            np.random.seed(i)
+            m.seed(4321)
+            # If m.state is not honored, the result will change
+            assert_array_equal(m.choice(10, size=10, p=np.ones(10) / 10.), res)
+
+    def test_multivariate_normal_size_types(self):
+        # Test for multivariate_normal issue with 'size' argument.
+        # Check that the multivariate_normal size argument can be a
+        # numpy integer.
+        np.random.multivariate_normal([0], [[0]], size=1)
+        np.random.multivariate_normal([0], [[0]], size=np.int_(1))
+        np.random.multivariate_normal([0], [[0]], size=np.int64(1))
+
+    def test_beta_small_parameters(self):
+        # Test that beta with small a and b parameters does not produce
+        # NaNs due to roundoff errors causing 0 / 0, gh-5851
+        np.random.seed(1234567890)
+        x = np.random.beta(0.0001, 0.0001, size=100)
+        assert_(not np.any(np.isnan(x)), 'Nans in np.random.beta')
+
+    def test_choice_sum_of_probs_tolerance(self):
+        # The sum of probs should be 1.0 with some tolerance.
+        # For low precision dtypes the tolerance was too tight.
+        # See numpy github issue 6123.
+        np.random.seed(1234)
+        a = [1, 2, 3]
+        counts = [4, 4, 2]
+        for dt in np.float16, np.float32, np.float64:
+            probs = np.array(counts, dtype=dt) / sum(counts)
+            c = np.random.choice(a, p=probs)
+            assert_(c in a)
+            assert_raises(ValueError, np.random.choice, a, p=probs * 0.9)
+
+    def test_shuffle_of_array_of_different_length_strings(self):
+        # Test that permuting an array of different length strings
+        # will not cause a segfault on garbage collection
+        # Tests gh-7710
+        np.random.seed(1234)
+
+        a = np.array(['a', 'a' * 1000])
+
+        for _ in range(100):
+            np.random.shuffle(a)
+
+        # Force Garbage Collection - should not segfault.
+        import gc
+        gc.collect()
+
+    def test_shuffle_of_array_of_objects(self):
+        # Test that permuting an array of objects will not cause
+        # a segfault on garbage collection.
+        # See gh-7719
+        np.random.seed(1234)
+        a = np.array([np.arange(1), np.arange(4)], dtype=object)
+
+        for _ in range(1000):
+            np.random.shuffle(a)
+
+        # Force Garbage Collection - should not segfault.
+        import gc
+        gc.collect()
+
+    def test_permutation_subclass(self):
+        class N(np.ndarray):
+            pass
+
+        rng = np.random.RandomState(1)
+        orig = np.arange(3).view(N)
+        perm = rng.permutation(orig)
+        assert_array_equal(perm, np.array([0, 2, 1]))
+        assert_array_equal(orig, np.arange(3).view(N))
+
+        class M:
+            a = np.arange(5)
+
+            def __array__(self, dtype=None, copy=None):
+                return self.a
+
+        rng = np.random.RandomState(1)
+        m = M()
+        perm = rng.permutation(m)
+        assert_array_equal(perm, np.array([2, 1, 4, 0, 3]))
+        assert_array_equal(m.__array__(), np.arange(5))
+
+    @pytest.mark.skipif(sys.flags.optimize == 2, reason="Python running -OO")
+    @pytest.mark.skipif(IS_PYPY, reason="PyPy does not modify tp_doc")
+    @pytest.mark.parametrize(
+        "cls",
+        [
+            random.Generator,
+            random.MT19937,
+            random.PCG64,
+            random.PCG64DXSM,
+            random.Philox,
+            random.RandomState,
+            random.SFC64,
+            random.BitGenerator,
+            random.SeedSequence,
+            random.bit_generator.SeedlessSeedSequence,
+        ],
+    )
+    def test_inspect_signature(self, cls: type) -> None:
+        assert hasattr(cls, "__text_signature__")
+        try:
+            inspect.signature(cls)
+        except ValueError:
+            pytest.fail(f"invalid signature: {cls.__module__}.{cls.__qualname__}")
@@ -0,0 +1,79 @@
+import numpy as np
+from numpy.random import SeedSequence
+from numpy.testing import assert_array_compare, assert_array_equal
+
+
+def test_reference_data():
+    """ Check that SeedSequence generates data the same as the C++ reference.
+
+    https://gist.github.com/imneme/540829265469e673d045
+    """
+    inputs = [
+        [3735928559, 195939070, 229505742, 305419896],
+        [3668361503, 4165561550, 1661411377, 3634257570],
+        [164546577, 4166754639, 1765190214, 1303880213],
+        [446610472, 3941463886, 522937693, 1882353782],
+        [1864922766, 1719732118, 3882010307, 1776744564],
+        [4141682960, 3310988675, 553637289, 902896340],
+        [1134851934, 2352871630, 3699409824, 2648159817],
+        [1240956131, 3107113773, 1283198141, 1924506131],
+        [2669565031, 579818610, 3042504477, 2774880435],
+        [2766103236, 2883057919, 4029656435, 862374500],
+    ]
+    outputs = [
+        [3914649087, 576849849, 3593928901, 2229911004],
+        [2240804226, 3691353228, 1365957195, 2654016646],
+        [3562296087, 3191708229, 1147942216, 3726991905],
+        [1403443605, 3591372999, 1291086759, 441919183],
+        [1086200464, 2191331643, 560336446, 3658716651],
+        [3249937430, 2346751812, 847844327, 2996632307],
+        [2584285912, 4034195531, 3523502488, 169742686],
+        [959045797, 3875435559, 1886309314, 359682705],
+        [3978441347, 432478529, 3223635119, 138903045],
+        [296367413, 4262059219, 13109864, 3283683422],
+    ]
+    outputs64 = [
+        [2477551240072187391, 9577394838764454085],
+        [15854241394484835714, 11398914698975566411],
+        [13708282465491374871, 16007308345579681096],
+        [15424829579845884309, 1898028439751125927],
+        [9411697742461147792, 15714068361935982142],
+        [10079222287618677782, 12870437757549876199],
+        [17326737873898640088, 729039288628699544],
+        [16644868984619524261, 1544825456798124994],
+        [1857481142255628931, 596584038813451439],
+        [18305404959516669237, 14103312907920476776],
+    ]
+    for seed, expected, expected64 in zip(inputs, outputs, outputs64):
+        expected = np.array(expected, dtype=np.uint32)
+        ss = SeedSequence(seed)
+        state = ss.generate_state(len(expected))
+        assert_array_equal(state, expected)
+        state64 = ss.generate_state(len(expected64), dtype=np.uint64)
+        assert_array_equal(state64, expected64)
+
+
+def test_zero_padding():
+    """ Ensure that the implicit zero-padding does not cause problems.
+    """
+    # Ensure that large integers are inserted in little-endian fashion to avoid
+    # trailing 0s.
+    ss0 = SeedSequence(42)
+    ss1 = SeedSequence(42 << 32)
+    assert_array_compare(
+        np.not_equal,
+        ss0.generate_state(4),
+        ss1.generate_state(4))
+
+    # Ensure backwards compatibility with the original 0.17 release for small
+    # integers and no spawn key.
+    expected42 = np.array([3444837047, 2669555309, 2046530742, 3581440988],
+                          dtype=np.uint32)
+    assert_array_equal(SeedSequence(42).generate_state(4), expected42)
+
+    # Regression test for gh-16539 to ensure that the implicit 0s don't
+    # conflict with spawn keys.
+    assert_array_compare(
+        np.not_equal,
+        SeedSequence(42, spawn_key=(0,)).generate_state(4),
+        expected42)
@@ -0,0 +1,882 @@
+import pickle
+from dataclasses import dataclass
+from functools import partial
+
+import pytest
+
+import numpy as np
+from numpy.random import MT19937, PCG64, PCG64DXSM, SFC64, Generator, Philox
+from numpy.testing import assert_, assert_array_equal, assert_equal
+
+DTYPES_BOOL_INT_UINT = (np.bool, np.int8, np.int16, np.int32, np.int64,
+                        np.uint8, np.uint16, np.uint32, np.uint64)
+
+
+def params_0(f):
+    val = f()
+    assert_(np.isscalar(val))
+    val = f(10)
+    assert_(val.shape == (10,))
+    val = f((10, 10))
+    assert_(val.shape == (10, 10))
+    val = f((10, 10, 10))
+    assert_(val.shape == (10, 10, 10))
+    val = f(size=(5, 5))
+    assert_(val.shape == (5, 5))
+
+
+def params_1(f, bounded=False):
+    a = 5.0
+    b = np.arange(2.0, 12.0)
+    c = np.arange(2.0, 102.0).reshape((10, 10))
+    d = np.arange(2.0, 1002.0).reshape((10, 10, 10))
+    e = np.array([2.0, 3.0])
+    g = np.arange(2.0, 12.0).reshape((1, 10, 1))
+    if bounded:
+        a = 0.5
+        b = b / (1.5 * b.max())
+        c = c / (1.5 * c.max())
+        d = d / (1.5 * d.max())
+        e = e / (1.5 * e.max())
+        g = g / (1.5 * g.max())
+
+    # Scalar
+    f(a)
+    # Scalar - size
+    f(a, size=(10, 10))
+    # 1d
+    f(b)
+    # 2d
+    f(c)
+    # 3d
+    f(d)
+    # 1d size
+    f(b, size=10)
+    # 2d - size - broadcast
+    f(e, size=(10, 2))
+    # 3d - size
+    f(g, size=(10, 10, 10))
+
+
+def comp_state(state1, state2):
+    identical = True
+    if isinstance(state1, dict):
+        for key in state1:
+            identical &= comp_state(state1[key], state2[key])
+    elif type(state1) != type(state2):
+        identical &= type(state1) == type(state2)
+    elif (isinstance(state1, (list, tuple, np.ndarray)) and isinstance(
+            state2, (list, tuple, np.ndarray))):
+        for s1, s2 in zip(state1, state2):
+            identical &= comp_state(s1, s2)
+    else:
+        identical &= state1 == state2
+    return identical
+
+
+def warmup(rg, n=None):
+    if n is None:
+        n = 11 + np.random.randint(0, 20)
+    rg.standard_normal(n)
+    rg.standard_normal(n)
+    rg.standard_normal(n, dtype=np.float32)
+    rg.standard_normal(n, dtype=np.float32)
+    rg.integers(0, 2 ** 24, n, dtype=np.uint64)
+    rg.integers(0, 2 ** 48, n, dtype=np.uint64)
+    rg.standard_gamma(11.0, n)
+    rg.standard_gamma(11.0, n, dtype=np.float32)
+    rg.random(n, dtype=np.float64)
+    rg.random(n, dtype=np.float32)
+
+
+@dataclass
+class RNGData:
+    bit_generator: type[np.random.BitGenerator]
+    advance: int
+    seed: list[int]
+    rg: Generator
+    seed_vector_bits: int
+
+
+class RNG:
+    @classmethod
+    def _create_rng(cls):
+        # Overridden in test classes. Place holder to silence IDE noise
+        bit_generator = PCG64
+        advance = None
+        seed = [12345]
+        rg = Generator(bit_generator(*seed))
+        seed_vector_bits = 64
+        return RNGData(bit_generator, advance, seed, rg, seed_vector_bits)
+
+    def test_init(self):
+        data = self._create_rng()
+        data.rg = Generator(data.bit_generator())
+        state = data.rg.bit_generator.state
+        data.rg.standard_normal(1)
+        data.rg.standard_normal(1)
+        data.rg.bit_generator.state = state
+        new_state = data.rg.bit_generator.state
+        assert_(comp_state(state, new_state))
+
+    def test_advance(self):
+        data = self._create_rng()
+        state = data.rg.bit_generator.state
+        if hasattr(data.rg.bit_generator, 'advance'):
+            data.rg.bit_generator.advance(data.advance)
+            assert_(not comp_state(state, data.rg.bit_generator.state))
+        else:
+            bitgen_name = data.rg.bit_generator.__class__.__name__
+            pytest.skip(f'Advance is not supported by {bitgen_name}')
+
+    def test_jump(self):
+        rg = self._create_rng().rg
+        state = rg.bit_generator.state
+        if hasattr(rg.bit_generator, 'jumped'):
+            bit_gen2 = rg.bit_generator.jumped()
+            jumped_state = bit_gen2.state
+            assert_(not comp_state(state, jumped_state))
+            rg.random(2 * 3 * 5 * 7 * 11 * 13 * 17)
+            rg.bit_generator.state = state
+            bit_gen3 = rg.bit_generator.jumped()
+            rejumped_state = bit_gen3.state
+            assert_(comp_state(jumped_state, rejumped_state))
+        else:
+            bitgen_name = rg.bit_generator.__class__.__name__
+            if bitgen_name not in ('SFC64',):
+                raise AttributeError(f'no "jumped" in {bitgen_name}')
+            pytest.skip(f'Jump is not supported by {bitgen_name}')
+
+    def test_uniform(self):
+        rg = self._create_rng().rg
+        r = rg.uniform(-1.0, 0.0, size=10)
+        assert_(len(r) == 10)
+        assert_((r > -1).all())
+        assert_((r <= 0).all())
+
+    def test_uniform_array(self):
+        rg = self._create_rng().rg
+        r = rg.uniform(np.array([-1.0] * 10), 0.0, size=10)
+        assert_(len(r) == 10)
+        assert_((r > -1).all())
+        assert_((r <= 0).all())
+        r = rg.uniform(np.array([-1.0] * 10),
+                            np.array([0.0] * 10), size=10)
+        assert_(len(r) == 10)
+        assert_((r > -1).all())
+        assert_((r <= 0).all())
+        r = rg.uniform(-1.0, np.array([0.0] * 10), size=10)
+        assert_(len(r) == 10)
+        assert_((r > -1).all())
+        assert_((r <= 0).all())
+
+    def test_random(self):
+        rg = self._create_rng().rg
+        assert_(len(rg.random(10)) == 10)
+        params_0(rg.random)
+
+    def test_standard_normal_zig(self):
+        rg = self._create_rng().rg
+        assert_(len(rg.standard_normal(10)) == 10)
+
+    def test_standard_normal(self):
+        rg = self._create_rng().rg
+        assert_(len(rg.standard_normal(10)) == 10)
+        params_0(rg.standard_normal)
+
+    def test_standard_gamma(self):
+        rg = self._create_rng().rg
+        assert_(len(rg.standard_gamma(10, 10)) == 10)
+        assert_(len(rg.standard_gamma(np.array([10] * 10), 10)) == 10)
+        params_1(rg.standard_gamma)
+
+    def test_standard_exponential(self):
+        rg = self._create_rng().rg
+        assert_(len(rg.standard_exponential(10)) == 10)
+        params_0(rg.standard_exponential)
+
+    def test_standard_exponential_float(self):
+        rg = self._create_rng().rg
+        randoms = rg.standard_exponential(10, dtype='float32')
+        assert_(len(randoms) == 10)
+        assert randoms.dtype == np.float32
+        params_0(partial(rg.standard_exponential, dtype='float32'))
+
+    def test_standard_exponential_float_log(self):
+        rg = self._create_rng().rg
+        randoms = rg.standard_exponential(10, dtype='float32',
+                                               method='inv')
+        assert_(len(randoms) == 10)
+        assert randoms.dtype == np.float32
+        params_0(partial(rg.standard_exponential, dtype='float32',
+                         method='inv'))
+
+    def test_standard_cauchy(self):
+        rg = self._create_rng().rg
+        assert_(len(rg.standard_cauchy(10)) == 10)
+        params_0(rg.standard_cauchy)
+
+    def test_standard_t(self):
+        rg = self._create_rng().rg
+        assert_(len(rg.standard_t(10, 10)) == 10)
+        params_1(rg.standard_t)
+
+    def test_binomial(self):
+        rg = self._create_rng().rg
+        assert_(rg.binomial(10, .5) >= 0)
+        assert_(rg.binomial(1000, .5) >= 0)
+
+    def test_reset_state(self):
+        rg = self._create_rng().rg
+        state = rg.bit_generator.state
+        int_1 = rg.integers(2**31)
+        rg.bit_generator.state = state
+        int_2 = rg.integers(2**31)
+        assert_(int_1 == int_2)
+
+    def test_entropy_init(self):
+        bit_generator = self._create_rng().bit_generator
+        rg = Generator(bit_generator())
+        rg2 = Generator(bit_generator())
+        assert_(not comp_state(rg.bit_generator.state,
+                               rg2.bit_generator.state))
+
+    def test_seed(self):
+        data = self._create_rng()
+        rg = Generator(data.bit_generator(*data.seed))
+        rg2 = Generator(data.bit_generator(*data.seed))
+        rg.random()
+        rg2.random()
+        assert_(comp_state(rg.bit_generator.state, rg2.bit_generator.state))
+
+    def test_reset_state_gauss(self):
+        data = self._create_rng()
+        rg = Generator(data.bit_generator(*data.seed))
+        rg.standard_normal()
+        state = rg.bit_generator.state
+        n1 = rg.standard_normal(size=10)
+        rg2 = Generator(data.bit_generator())
+        rg2.bit_generator.state = state
+        n2 = rg2.standard_normal(size=10)
+        assert_array_equal(n1, n2)
+
+    def test_reset_state_uint32(self):
+        data = self._create_rng()
+        rg = Generator(data.bit_generator(*data.seed))
+        rg.integers(0, 2 ** 24, 120, dtype=np.uint32)
+        state = rg.bit_generator.state
+        n1 = rg.integers(0, 2 ** 24, 10, dtype=np.uint32)
+        rg2 = Generator(data.bit_generator())
+        rg2.bit_generator.state = state
+        n2 = rg2.integers(0, 2 ** 24, 10, dtype=np.uint32)
+        assert_array_equal(n1, n2)
+
+    def test_reset_state_float(self):
+        data = self._create_rng()
+        rg = Generator(data.bit_generator(*data.seed))
+        rg.random(dtype='float32')
+        state = rg.bit_generator.state
+        n1 = rg.random(size=10, dtype='float32')
+        rg2 = Generator(data.bit_generator())
+        rg2.bit_generator.state = state
+        n2 = rg2.random(size=10, dtype='float32')
+        assert_((n1 == n2).all())
+
+    def test_shuffle(self):
+        rg = self._create_rng().rg
+        original = np.arange(200, 0, -1)
+        permuted = rg.permutation(original)
+        assert_((original != permuted).any())
+
+    def test_permutation(self):
+        rg = self._create_rng().rg
+        original = np.arange(200, 0, -1)
+        permuted = rg.permutation(original)
+        assert_((original != permuted).any())
+
+    def test_beta(self):
+        rg = self._create_rng().rg
+        vals = rg.beta(2.0, 2.0, 10)
+        assert_(len(vals) == 10)
+        vals = rg.beta(np.array([2.0] * 10), 2.0)
+        assert_(len(vals) == 10)
+        vals = rg.beta(2.0, np.array([2.0] * 10))
+        assert_(len(vals) == 10)
+        vals = rg.beta(np.array([2.0] * 10), np.array([2.0] * 10))
+        assert_(len(vals) == 10)
+        vals = rg.beta(np.array([2.0] * 10), np.array([[2.0]] * 10))
+        assert_(vals.shape == (10, 10))
+
+    def test_bytes(self):
+        rg = self._create_rng().rg
+        vals = rg.bytes(10)
+        assert_(len(vals) == 10)
+
+    def test_chisquare(self):
+        rg = self._create_rng().rg
+        vals = rg.chisquare(2.0, 10)
+        assert_(len(vals) == 10)
+        params_1(rg.chisquare)
+
+    def test_exponential(self):
+        rg = self._create_rng().rg
+        vals = rg.exponential(2.0, 10)
+        assert_(len(vals) == 10)
+        params_1(rg.exponential)
+
+    def test_f(self):
+        rg = self._create_rng().rg
+        vals = rg.f(3, 1000, 10)
+        assert_(len(vals) == 10)
+
+    def test_gamma(self):
+        rg = self._create_rng().rg
+        vals = rg.gamma(3, 2, 10)
+        assert_(len(vals) == 10)
+
+    def test_geometric(self):
+        rg = self._create_rng().rg
+        vals = rg.geometric(0.5, 10)
+        assert_(len(vals) == 10)
+        params_1(rg.exponential, bounded=True)
+
+    def test_gumbel(self):
+        rg = self._create_rng().rg
+        vals = rg.gumbel(2.0, 2.0, 10)
+        assert_(len(vals) == 10)
+
+    def test_laplace(self):
+        rg = self._create_rng().rg
+        vals = rg.laplace(2.0, 2.0, 10)
+        assert_(len(vals) == 10)
+
+    def test_logitic(self):
+        rg = self._create_rng().rg
+        vals = rg.logistic(2.0, 2.0, 10)
+        assert_(len(vals) == 10)
+
+    def test_logseries(self):
+        rg = self._create_rng().rg
+        vals = rg.logseries(0.5, 10)
+        assert_(len(vals) == 10)
+
+    def test_negative_binomial(self):
+        rg = self._create_rng().rg
+        vals = rg.negative_binomial(10, 0.2, 10)
+        assert_(len(vals) == 10)
+
+    def test_noncentral_chisquare(self):
+        rg = self._create_rng().rg
+        vals = rg.noncentral_chisquare(10, 2, 10)
+        assert_(len(vals) == 10)
+
+    def test_noncentral_f(self):
+        rg = self._create_rng().rg
+        vals = rg.noncentral_f(3, 1000, 2, 10)
+        assert_(len(vals) == 10)
+        vals = rg.noncentral_f(np.array([3] * 10), 1000, 2)
+        assert_(len(vals) == 10)
+        vals = rg.noncentral_f(3, np.array([1000] * 10), 2)
+        assert_(len(vals) == 10)
+        vals = rg.noncentral_f(3, 1000, np.array([2] * 10))
+        assert_(len(vals) == 10)
+
+    def test_normal(self):
+        rg = self._create_rng().rg
+        vals = rg.normal(10, 0.2, 10)
+        assert_(len(vals) == 10)
+
+    def test_pareto(self):
+        rg = self._create_rng().rg
+        vals = rg.pareto(3.0, 10)
+        assert_(len(vals) == 10)
+
+    def test_poisson(self):
+        rg = self._create_rng().rg
+        vals = rg.poisson(10, 10)
+        assert_(len(vals) == 10)
+        vals = rg.poisson(np.array([10] * 10))
+        assert_(len(vals) == 10)
+        params_1(rg.poisson)
+
+    def test_power(self):
+        rg = self._create_rng().rg
+        vals = rg.power(0.2, 10)
+        assert_(len(vals) == 10)
+
+    def test_integers(self):
+        rg = self._create_rng().rg
+        vals = rg.integers(10, 20, 10)
+        assert_(len(vals) == 10)
+
+    def test_rayleigh(self):
+        rg = self._create_rng().rg
+        vals = rg.rayleigh(0.2, 10)
+        assert_(len(vals) == 10)
+        params_1(rg.rayleigh, bounded=True)
+
+    def test_vonmises(self):
+        rg = self._create_rng().rg
+        vals = rg.vonmises(10, 0.2, 10)
+        assert_(len(vals) == 10)
+
+    def test_wald(self):
+        rg = self._create_rng().rg
+        vals = rg.wald(1.0, 1.0, 10)
+        assert_(len(vals) == 10)
+
+    def test_weibull(self):
+        rg = self._create_rng().rg
+        vals = rg.weibull(1.0, 10)
+        assert_(len(vals) == 10)
+
+    def test_zipf(self):
+        rg = self._create_rng().rg
+        vec_1d = np.arange(2.0, 102.0)
+        vec_2d = np.arange(2.0, 102.0)[None, :]
+        mat = np.arange(2.0, 102.0, 0.01).reshape((100, 100))
+        vals = rg.zipf(10, 10)
+        assert_(len(vals) == 10)
+        vals = rg.zipf(vec_1d)
+        assert_(len(vals) == 100)
+        vals = rg.zipf(vec_2d)
+        assert_(vals.shape == (1, 100))
+        vals = rg.zipf(mat)
+        assert_(vals.shape == (100, 100))
+
+    def test_hypergeometric(self):
+        rg = self._create_rng().rg
+        vals = rg.hypergeometric(25, 25, 20)
+        assert_(np.isscalar(vals))
+        vals = rg.hypergeometric(np.array([25] * 10), 25, 20)
+        assert_(vals.shape == (10,))
+
+    def test_triangular(self):
+        rg = self._create_rng().rg
+        vals = rg.triangular(-5, 0, 5)
+        assert_(np.isscalar(vals))
+        vals = rg.triangular(-5, np.array([0] * 10), 5)
+        assert_(vals.shape == (10,))
+
+    def test_multivariate_normal(self):
+        rg = self._create_rng().rg
+        mean = [0, 0]
+        cov = [[1, 0], [0, 100]]  # diagonal covariance
+        x = rg.multivariate_normal(mean, cov, 5000)
+        assert_(x.shape == (5000, 2))
+        x_zig = rg.multivariate_normal(mean, cov, 5000)
+        assert_(x.shape == (5000, 2))
+        x_inv = rg.multivariate_normal(mean, cov, 5000)
+        assert_(x.shape == (5000, 2))
+        assert_((x_zig != x_inv).any())
+
+    def test_multinomial(self):
+        rg = self._create_rng().rg
+        vals = rg.multinomial(100, [1.0 / 3, 2.0 / 3])
+        assert_(vals.shape == (2,))
+        vals = rg.multinomial(100, [1.0 / 3, 2.0 / 3], size=10)
+        assert_(vals.shape == (10, 2))
+
+    def test_dirichlet(self):
+        rg = self._create_rng().rg
+        s = rg.dirichlet((10, 5, 3), 20)
+        assert_(s.shape == (20, 3))
+
+    def test_pickle(self):
+        rg = self._create_rng().rg
+        pick = pickle.dumps(rg)
+        unpick = pickle.loads(pick)
+        assert_(type(rg) == type(unpick))
+        assert_(comp_state(rg.bit_generator.state,
+                           unpick.bit_generator.state))
+
+        pick = pickle.dumps(rg)
+        unpick = pickle.loads(pick)
+        assert_(type(rg) == type(unpick))
+        assert_(comp_state(rg.bit_generator.state,
+                           unpick.bit_generator.state))
+
+    def test_seed_array(self):
+        data = self._create_rng()
+        if data.seed_vector_bits is None:
+            bitgen_name = data.bit_generator.__name__
+            pytest.skip(f'Vector seeding is not supported by {bitgen_name}')
+
+        if data.seed_vector_bits == 32:
+            dtype = np.uint32
+        else:
+            dtype = np.uint64
+        seed = np.array([1], dtype=dtype)
+        bg = data.bit_generator(seed)
+        state1 = bg.state
+        bg = data.bit_generator(1)
+        state2 = bg.state
+        assert_(comp_state(state1, state2))
+
+        seed = np.arange(4, dtype=dtype)
+        bg = data.bit_generator(seed)
+        state1 = bg.state
+        bg = data.bit_generator(seed[0])
+        state2 = bg.state
+        assert_(not comp_state(state1, state2))
+
+        seed = np.arange(1500, dtype=dtype)
+        bg = data.bit_generator(seed)
+        state1 = bg.state
+        bg = data.bit_generator(seed[0])
+        state2 = bg.state
+        assert_(not comp_state(state1, state2))
+
+        seed = 2 ** np.mod(np.arange(1500, dtype=dtype),
+                           data.seed_vector_bits - 1) + 1
+        bg = data.bit_generator(seed)
+        state1 = bg.state
+        bg = data.bit_generator(seed[0])
+        state2 = bg.state
+        assert_(not comp_state(state1, state2))
+
+    def test_uniform_float(self):
+        bit_generator = self._create_rng().bit_generator
+        rg = Generator(bit_generator(12345))
+        warmup(rg)
+        state = rg.bit_generator.state
+        r1 = rg.random(11, dtype=np.float32)
+        rg2 = Generator(bit_generator())
+        warmup(rg2)
+        rg2.bit_generator.state = state
+        r2 = rg2.random(11, dtype=np.float32)
+        assert_array_equal(r1, r2)
+        assert_equal(r1.dtype, np.float32)
+        assert_(comp_state(rg.bit_generator.state, rg2.bit_generator.state))
+
+    def test_gamma_floats(self):
+        bit_generator = self._create_rng().bit_generator
+        rg = Generator(bit_generator())
+        warmup(rg)
+        state = rg.bit_generator.state
+        r1 = rg.standard_gamma(4.0, 11, dtype=np.float32)
+        rg2 = Generator(bit_generator())
+        warmup(rg2)
+        rg2.bit_generator.state = state
+        r2 = rg2.standard_gamma(4.0, 11, dtype=np.float32)
+        assert_array_equal(r1, r2)
+        assert_equal(r1.dtype, np.float32)
+        assert_(comp_state(rg.bit_generator.state, rg2.bit_generator.state))
+
+    def test_normal_floats(self):
+        bit_generator = self._create_rng().bit_generator
+        rg = Generator(bit_generator())
+        warmup(rg)
+        state = rg.bit_generator.state
+        r1 = rg.standard_normal(11, dtype=np.float32)
+        rg2 = Generator(bit_generator())
+        warmup(rg2)
+        rg2.bit_generator.state = state
+        r2 = rg2.standard_normal(11, dtype=np.float32)
+        assert_array_equal(r1, r2)
+        assert_equal(r1.dtype, np.float32)
+        assert_(comp_state(rg.bit_generator.state, rg2.bit_generator.state))
+
+    def test_normal_zig_floats(self):
+        bit_generator = self._create_rng().bit_generator
+        rg = Generator(bit_generator())
+        warmup(rg)
+        state = rg.bit_generator.state
+        r1 = rg.standard_normal(11, dtype=np.float32)
+        rg2 = Generator(bit_generator())
+        warmup(rg2)
+        rg2.bit_generator.state = state
+        r2 = rg2.standard_normal(11, dtype=np.float32)
+        assert_array_equal(r1, r2)
+        assert_equal(r1.dtype, np.float32)
+        assert_(comp_state(rg.bit_generator.state, rg2.bit_generator.state))
+
+    def test_output_fill(self):
+        rg = self._create_rng().rg
+        state = rg.bit_generator.state
+        size = (31, 7, 97)
+        existing = np.empty(size)
+        rg.bit_generator.state = state
+        rg.standard_normal(out=existing)
+        rg.bit_generator.state = state
+        direct = rg.standard_normal(size=size)
+        assert_equal(direct, existing)
+
+        sized = np.empty(size)
+        rg.bit_generator.state = state
+        rg.standard_normal(out=sized, size=sized.shape)
+
+        existing = np.empty(size, dtype=np.float32)
+        rg.bit_generator.state = state
+        rg.standard_normal(out=existing, dtype=np.float32)
+        rg.bit_generator.state = state
+        direct = rg.standard_normal(size=size, dtype=np.float32)
+        assert_equal(direct, existing)
+
+    def test_output_filling_uniform(self):
+        rg = self._create_rng().rg
+        state = rg.bit_generator.state
+        size = (31, 7, 97)
+        existing = np.empty(size)
+        rg.bit_generator.state = state
+        rg.random(out=existing)
+        rg.bit_generator.state = state
+        direct = rg.random(size=size)
+        assert_equal(direct, existing)
+
+        existing = np.empty(size, dtype=np.float32)
+        rg.bit_generator.state = state
+        rg.random(out=existing, dtype=np.float32)
+        rg.bit_generator.state = state
+        direct = rg.random(size=size, dtype=np.float32)
+        assert_equal(direct, existing)
+
+    def test_output_filling_exponential(self):
+        rg = self._create_rng().rg
+        state = rg.bit_generator.state
+        size = (31, 7, 97)
+        existing = np.empty(size)
+        rg.bit_generator.state = state
+        rg.standard_exponential(out=existing)
+        rg.bit_generator.state = state
+        direct = rg.standard_exponential(size=size)
+        assert_equal(direct, existing)
+
+        existing = np.empty(size, dtype=np.float32)
+        rg.bit_generator.state = state
+        rg.standard_exponential(out=existing, dtype=np.float32)
+        rg.bit_generator.state = state
+        direct = rg.standard_exponential(size=size, dtype=np.float32)
+        assert_equal(direct, existing)
+
+    def test_output_filling_gamma(self):
+        rg = self._create_rng().rg
+        state = rg.bit_generator.state
+        size = (31, 7, 97)
+        existing = np.zeros(size)
+        rg.bit_generator.state = state
+        rg.standard_gamma(1.0, out=existing)
+        rg.bit_generator.state = state
+        direct = rg.standard_gamma(1.0, size=size)
+        assert_equal(direct, existing)
+
+        existing = np.zeros(size, dtype=np.float32)
+        rg.bit_generator.state = state
+        rg.standard_gamma(1.0, out=existing, dtype=np.float32)
+        rg.bit_generator.state = state
+        direct = rg.standard_gamma(1.0, size=size, dtype=np.float32)
+        assert_equal(direct, existing)
+
+    def test_output_filling_gamma_broadcast(self):
+        rg = self._create_rng().rg
+        state = rg.bit_generator.state
+        size = (31, 7, 97)
+        mu = np.arange(97.0) + 1.0
+        existing = np.zeros(size)
+        rg.bit_generator.state = state
+        rg.standard_gamma(mu, out=existing)
+        rg.bit_generator.state = state
+        direct = rg.standard_gamma(mu, size=size)
+        assert_equal(direct, existing)
+
+        existing = np.zeros(size, dtype=np.float32)
+        rg.bit_generator.state = state
+        rg.standard_gamma(mu, out=existing, dtype=np.float32)
+        rg.bit_generator.state = state
+        direct = rg.standard_gamma(mu, size=size, dtype=np.float32)
+        assert_equal(direct, existing)
+
+    def test_output_fill_error(self):
+        rg = self._create_rng().rg
+        size = (31, 7, 97)
+        existing = np.empty(size)
+        with pytest.raises(TypeError):
+            rg.standard_normal(out=existing, dtype=np.float32)
+        with pytest.raises(ValueError):
+            rg.standard_normal(out=existing[::3])
+        existing = np.empty(size, dtype=np.float32)
+        with pytest.raises(TypeError):
+            rg.standard_normal(out=existing, dtype=np.float64)
+
+        existing = np.zeros(size, dtype=np.float32)
+        with pytest.raises(TypeError):
+            rg.standard_gamma(1.0, out=existing, dtype=np.float64)
+        with pytest.raises(ValueError):
+            rg.standard_gamma(1.0, out=existing[::3], dtype=np.float32)
+        existing = np.zeros(size, dtype=np.float64)
+        with pytest.raises(TypeError):
+            rg.standard_gamma(1.0, out=existing, dtype=np.float32)
+        with pytest.raises(ValueError):
+            rg.standard_gamma(1.0, out=existing[::3])
+
+    @pytest.mark.parametrize("dtype", DTYPES_BOOL_INT_UINT)
+    def test_integers_broadcast(self, dtype):
+        rg = self._create_rng().rg
+        initial_state = rg.bit_generator.state
+
+        def reset_state(rng):
+            rng.bit_generator.state = initial_state
+
+        if dtype == np.bool:
+            upper = 2
+            lower = 0
+        else:
+            info = np.iinfo(dtype)
+            upper = int(info.max) + 1
+            lower = info.min
+        reset_state(rg)
+        rg.bit_generator.state = initial_state
+        a = rg.integers(lower, [upper] * 10, dtype=dtype)
+        reset_state(rg)
+        b = rg.integers([lower] * 10, upper, dtype=dtype)
+        assert_equal(a, b)
+        reset_state(rg)
+        c = rg.integers(lower, upper, size=10, dtype=dtype)
+        assert_equal(a, c)
+        reset_state(rg)
+        d = rg.integers(np.array(
+            [lower] * 10), np.array([upper], dtype=object), size=10,
+            dtype=dtype)
+        assert_equal(a, d)
+        reset_state(rg)
+        e = rg.integers(
+            np.array([lower] * 10), np.array([upper] * 10), size=10,
+            dtype=dtype)
+        assert_equal(a, e)
+
+        reset_state(rg)
+        a = rg.integers(0, upper, size=10, dtype=dtype)
+        reset_state(rg)
+        b = rg.integers([upper] * 10, dtype=dtype)
+        assert_equal(a, b)
+
+    @pytest.mark.parametrize("dtype", DTYPES_BOOL_INT_UINT)
+    def test_integers_numpy(self, dtype):
+        rg = self._create_rng().rg
+        high = np.array([1])
+        low = np.array([0])
+
+        out = rg.integers(low, high, dtype=dtype)
+        assert out.shape == (1,)
+
+        out = rg.integers(low[0], high, dtype=dtype)
+        assert out.shape == (1,)
+
+        out = rg.integers(low, high[0], dtype=dtype)
+        assert out.shape == (1,)
+
+    @pytest.mark.parametrize("dtype", DTYPES_BOOL_INT_UINT)
+    def test_integers_broadcast_errors(self, dtype):
+        rg = self._create_rng().rg
+        if dtype == np.bool:
+            upper = 2
+            lower = 0
+        else:
+            info = np.iinfo(dtype)
+            upper = int(info.max) + 1
+            lower = info.min
+        with pytest.raises(ValueError):
+            rg.integers(lower, [upper + 1] * 10, dtype=dtype)
+        with pytest.raises(ValueError):
+            rg.integers(lower - 1, [upper] * 10, dtype=dtype)
+        with pytest.raises(ValueError):
+            rg.integers([lower - 1], [upper] * 10, dtype=dtype)
+        with pytest.raises(ValueError):
+            rg.integers([0], [0], dtype=dtype)
+
+
+class TestMT19937(RNG):
+    @classmethod
+    def _create_rng(cls):
+        bit_generator = MT19937
+        advance = None
+        seed = [2 ** 21 + 2 ** 16 + 2 ** 5 + 1]
+        rg = Generator(bit_generator(*seed))
+        seed_vector_bits = 32
+        return RNGData(bit_generator, advance, seed, rg, seed_vector_bits)
+
+    def test_numpy_state(self):
+        rg = self._create_rng().rg
+        nprg = np.random.RandomState()
+        nprg.standard_normal(99)
+        state = nprg.get_state()
+        rg.bit_generator.state = state
+        state2 = rg.bit_generator.state
+        assert_((state[1] == state2['state']['key']).all())
+        assert_(state[2] == state2['state']['pos'])
+
+
+class TestPhilox(RNG):
+    @classmethod
+    def _create_rng(cls):
+        bit_generator = Philox
+        advance = 2**63 + 2**31 + 2**15 + 1
+        seed = [12345]
+        rg = Generator(bit_generator(*seed))
+        seed_vector_bits = 64
+        return RNGData(bit_generator, advance, seed, rg, seed_vector_bits)
+
+
+class TestSFC64(RNG):
+    @classmethod
+    def _create_rng(cls):
+        bit_generator = SFC64
+        advance = None
+        seed = [12345]
+        rg = Generator(bit_generator(*seed))
+        seed_vector_bits = 192
+        return RNGData(bit_generator, advance, seed, rg, seed_vector_bits)
+
+
+class TestPCG64(RNG):
+    @classmethod
+    def _create_rng(cls):
+        bit_generator = PCG64
+        advance = 2**63 + 2**31 + 2**15 + 1
+        seed = [12345]
+        rg = Generator(bit_generator(*seed))
+        seed_vector_bits = 64
+        return RNGData(bit_generator, advance, seed, rg, seed_vector_bits)
+
+
+class TestPCG64DXSM(RNG):
+    @classmethod
+    def _create_rng(cls):
+        bit_generator = PCG64DXSM
+        advance = 2**63 + 2**31 + 2**15 + 1
+        seed = [12345]
+        rg = Generator(bit_generator(*seed))
+        seed_vector_bits = 64
+        return RNGData(bit_generator, advance, seed, rg, seed_vector_bits)
+
+
+class TestDefaultRNG(RNG):
+    @classmethod
+    def _create_rng(cls):
+        # This will duplicate some tests that directly instantiate a fresh
+        # Generator(), but that's okay.
+        bit_generator = PCG64
+        advance = 2**63 + 2**31 + 2**15 + 1
+        seed = [12345]
+        rg = np.random.default_rng(*seed)
+        seed_vector_bits = 64
+        return RNGData(bit_generator, advance, seed, rg, seed_vector_bits)
+
+    def test_default_is_pcg64(self):
+        # In order to change the default BitGenerator, we'll go through
+        # a deprecation cycle to move to a different function.
+        rg = self._create_rng().rg
+        assert_(isinstance(rg.bit_generator, PCG64))
+
+    def test_seed(self):
+        np.random.default_rng()
+        np.random.default_rng(None)
+        np.random.default_rng(12345)
+        np.random.default_rng(0)
+        np.random.default_rng(43660444402423911716352051725018508569)
+        np.random.default_rng([43660444402423911716352051725018508569,
+                               279705150948142787361475340226491943209])
+        with pytest.raises(ValueError):
+            np.random.default_rng(-1)
+        with pytest.raises(ValueError):
+            np.random.default_rng([12345, -1])