diff --git a/demos/Wiener_Filter.ipynb b/demos/Wiener_Filter.ipynb
index b0ef1f8b38a3e372d0c120a45bceaabe287e35c0..ba46b6958af7cb8283eb71f5b1d401f166fc9447 100644
--- a/demos/Wiener_Filter.ipynb
+++ b/demos/Wiener_Filter.ipynb
@@ -140,7 +140,6 @@
    "outputs": [],
    "source": [
     "import numpy as np\n",
-    "np.random.seed(40)\n",
     "import nifty6 as ift\n",
     "import matplotlib.pyplot as plt\n",
     "%matplotlib inline"
diff --git a/demos/bench_gridder.py b/demos/bench_gridder.py
index a2a352bbb56e940694cdc5f0822ae7fcd3ebac59..ed8c5d74990eb077c083477af450a4d5f8e05306 100644
--- a/demos/bench_gridder.py
+++ b/demos/bench_gridder.py
@@ -5,8 +5,6 @@ import numpy as np
 
 import nifty6 as ift
 
-np.random.seed(40)
-
 N0s, a0s, b0s, c0s = [], [], [], []
 
 for ii in range(10, 26):
@@ -15,15 +13,15 @@ for ii in range(10, 26):
     N = int(2**ii)
     print('N = {}'.format(N))
 
-    uv = np.random.rand(N, 2) - 0.5
-    vis = np.random.randn(N) + 1j*np.random.randn(N)
+    rng = ift.random.current_rng()
+    uv = rng.uniform(-.5, .5, (N,2))
+    vis = rng.normal(0., 1., N) + 1j*rng.normal(0., 1., N)
 
     uvspace = ift.RGSpace((nu, nv))
 
     visspace = ift.UnstructuredDomain(N)
 
-    img = np.random.randn(nu*nv)
-    img = img.reshape((nu, nv))
+    img = rng.standard_normal((nu, nv))
     img = ift.makeField(uvspace, img)
 
     t0 = time()
diff --git a/demos/bernoulli_demo.py b/demos/bernoulli_demo.py
index c98dc2a71693c4aaf6e1ba125bf5067e5ec102de..76fec1a116c4f7d804ff108e6935645d6e4e06ef 100644
--- a/demos/bernoulli_demo.py
+++ b/demos/bernoulli_demo.py
@@ -27,8 +27,6 @@ import numpy as np
 import nifty6 as ift
 
 if __name__ == '__main__':
-    np.random.seed(41)
-
     # Set up the position space of the signal
     mode = 2
     if mode == 0:
@@ -62,7 +60,7 @@ if __name__ == '__main__':
     p = R(sky)
     mock_position = ift.from_random('normal', harmonic_space)
     tmp = p(mock_position).val.astype(np.float64)
-    data = np.random.binomial(1, tmp)
+    data = ift.random.current_rng().binomial(1, tmp)
     data = ift.Field.from_raw(R.target, data)
 
     # Compute likelihood and Hamiltonian
diff --git a/demos/getting_started_1.py b/demos/getting_started_1.py
index 886a4ceaa7362f7ae73d90a92221536e45e359d1..3249055044026a8f833f3ab7d2250fc0078cf77c 100644
--- a/demos/getting_started_1.py
+++ b/demos/getting_started_1.py
@@ -46,8 +46,6 @@ def make_random_mask():
 
 
 if __name__ == '__main__':
-    np.random.seed(42)
-
     # Choose space on which the signal field is defined
     if len(sys.argv) == 2:
         mode = int(sys.argv[1])
diff --git a/demos/getting_started_2.py b/demos/getting_started_2.py
index 38b568de571b6a21c74e24f3f486c0816a8c2e67..7856049d9efac5bb6ca51ce3745be05f0b5a446a 100644
--- a/demos/getting_started_2.py
+++ b/demos/getting_started_2.py
@@ -44,8 +44,6 @@ def exposure_2d():
 
 
 if __name__ == '__main__':
-    np.random.seed(42)
-
     # Choose space on which the signal field is defined
     if len(sys.argv) == 2:
         mode = int(sys.argv[1])
@@ -94,7 +92,7 @@ if __name__ == '__main__':
     lamb = R(sky)
     mock_position = ift.from_random('normal', domain)
     data = lamb(mock_position)
-    data = np.random.poisson(data.val.astype(np.float64))
+    data = ift.random.current_rng().poisson(data.val.astype(np.float64))
     data = ift.Field.from_raw(d_space, data)
     likelihood = ift.PoissonianEnergy(data) @ lamb
 
diff --git a/demos/getting_started_3.py b/demos/getting_started_3.py
index c3fdbb83b0206ec923f29366879399bc1910ccbd..21f5eac1c2006d1505475f825b829b76cc897e08 100644
--- a/demos/getting_started_3.py
+++ b/demos/getting_started_3.py
@@ -33,20 +33,18 @@ import nifty6 as ift
 
 
 def random_los(n_los):
-    starts = list(np.random.uniform(0, 1, (n_los, 2)).T)
-    ends = list(np.random.uniform(0, 1, (n_los, 2)).T)
+    starts = list(ift.random.current_rng().random((n_los, 2)).T)
+    ends = list(ift.random.current_rng().random((n_los, 2)).T)
     return starts, ends
 
 
 def radial_los(n_los):
-    starts = list(np.random.uniform(0, 1, (n_los, 2)).T)
-    ends = list(0.5 + 0*np.random.uniform(0, 1, (n_los, 2)).T)
+    starts = list(ift.random.current_rng().random((n_los, 2)).T)
+    ends = list(0.5 + 0*ift.random.current_rng().random((n_los, 2)).T)
     return starts, ends
 
 
 if __name__ == '__main__':
-    np.random.seed(420)
-
     # Choose between random line-of-sight response (mode=0) and radial lines
     # of sight (mode=1)
     if len(sys.argv) == 2:
diff --git a/demos/getting_started_mf.py b/demos/getting_started_mf.py
index d246663394d189c759484825000e78b9d179f38a..555c0623d1cf671902d046c25053af2f8e19df81 100644
--- a/demos/getting_started_mf.py
+++ b/demos/getting_started_mf.py
@@ -44,20 +44,18 @@ class SingleDomain(ift.LinearOperator):
 
 
 def random_los(n_los):
-    starts = list(np.random.uniform(0, 1, (n_los, 2)).T)
-    ends = list(np.random.uniform(0, 1, (n_los, 2)).T)
+    starts = list(ift.random.current_rng().random((n_los, 2)).T)
+    ends = list(ift.random.current_rng().random((n_los, 2)).T)
     return starts, ends
 
 
 def radial_los(n_los):
-    starts = list(np.random.uniform(0, 1, (n_los, 2)).T)
-    ends = list(0.5 + 0*np.random.uniform(0, 1, (n_los, 2)).T)
+    starts = list(ift.random.current_rng().random((n_los, 2)).T)
+    ends = list(0.5 + 0*ift.random.current_rng().random((n_los, 2)).T)
     return starts, ends
 
 
 if __name__ == '__main__':
-    np.random.seed(43)
-
     # Choose between random line-of-sight response (mode=0) and radial lines
     # of sight (mode=1)
     if len(sys.argv) == 2:
diff --git a/demos/polynomial_fit.py b/demos/polynomial_fit.py
index d64b7a1131a088f2dbe26929da5817c702e104e1..2db06700ca9f6067726ed4e8383a95e04c579acf 100644
--- a/demos/polynomial_fit.py
+++ b/demos/polynomial_fit.py
@@ -19,7 +19,6 @@ import matplotlib.pyplot as plt
 import numpy as np
 
 import nifty6 as ift
-np.random.seed(12)
 
 
 def polynomial(coefficients, sampling_points):
@@ -86,7 +85,7 @@ if __name__ == '__main__':
     N_params = 10
     N_samples = 100
     size = (12,)
-    x = np.random.random(size) * 10
+    x = ift.random.current_rng().random(size) * 10
     y = np.sin(x**2) * x**3
     var = np.full_like(y, y.var() / 10)
     var[-2] *= 4
diff --git a/nifty6/__init__.py b/nifty6/__init__.py
index 7ba7583e984d35422c16ca02ea48fdcfbeef6207..e2e8ec76b78d439354f3b6ef83a1b96ec2d21157 100644
--- a/nifty6/__init__.py
+++ b/nifty6/__init__.py
@@ -1,5 +1,7 @@
 from .version import __version__
 
+from . import random
+
 from .domains.domain import Domain
 from .domains.structured_domain import StructuredDomain
 from .domains.unstructured_domain import UnstructuredDomain
diff --git a/nifty6/library/special_distributions.py b/nifty6/library/special_distributions.py
index 9d4f65f86e1fb78fe9f6f0a7dfae456b9edbd67e..5d0d4ad7a64d98695488ac163a6bc428b55df85a 100644
--- a/nifty6/library/special_distributions.py
+++ b/nifty6/library/special_distributions.py
@@ -25,6 +25,7 @@ from ..field import Field
 from ..linearization import Linearization
 from ..operators.operator import Operator
 from ..sugar import makeOp
+from .. import random
 
 
 def _f_on_np(f, arr):
@@ -67,7 +68,8 @@ class _InterpolationOperator(Operator):
         if table_func is not None:
             if inv_table_func is None:
                 raise ValueError
-            a = func(np.random.randn(10))
+# MR FIXME: not sure whether we should have this in production code
+            a = func(random.current_rng().random(10))
             a1 = _f_on_np(lambda x: inv_table_func(table_func(x)), a)
             np.testing.assert_allclose(a, a1)
             self._table = _f_on_np(table_func, self._table)
diff --git a/nifty6/minimization/metric_gaussian_kl_mpi.py b/nifty6/minimization/metric_gaussian_kl_mpi.py
index 6054140f099d7312d94ad718b04129e4386e2ecc..24bd12980654218af6ba1917b475a3062367944c 100644
--- a/nifty6/minimization/metric_gaussian_kl_mpi.py
+++ b/nifty6/minimization/metric_gaussian_kl_mpi.py
@@ -11,7 +11,7 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
-# Copyright(C) 2013-2019 Max-Planck-Society
+# Copyright(C) 2013-2020 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
@@ -113,10 +113,6 @@ class MetricGaussianKL_MPI(Energy):
         preconditioning is done by default.
     _samples : None
         Only a parameter for internal uses. Typically not to be set by users.
-    seed_offset : int
-        A parameter with which one can controll from which seed the samples
-        are drawn. Per default, the seed is different for MPI tasks, but the
-        same every time this class is initialized.
 
     Note
     ----
@@ -132,7 +128,7 @@ class MetricGaussianKL_MPI(Energy):
 
     def __init__(self, mean, hamiltonian, n_samples, constants=[],
                  point_estimates=[], mirror_samples=False,
-                 napprox=0, _samples=None, seed_offset=0,
+                 napprox=0, _samples=None,
                  lh_sampling_dtype=np.float64):
         super(MetricGaussianKL_MPI, self).__init__(mean)
 
@@ -159,10 +155,10 @@ class MetricGaussianKL_MPI(Energy):
             if napprox > 1:
                 met._approximation = makeOp(approximation2endo(met, napprox))
             _samples = []
-            rand_state = np.random.get_state()
+            sseq = random.spawn_sseq(n_samples)
             for i in range(lo, hi):
                 if mirror_samples:
-                    np.random.seed(i//2+seed_offset)
+                    random.push_sseq(sseq[i//2])
                     if (i % 2) and (i-1 >= lo):
                         _samples.append(-_samples[-1])
 
@@ -170,16 +166,16 @@ class MetricGaussianKL_MPI(Energy):
                         _samples.append(((i % 2)*2-1) *
                                         met.draw_sample(from_inverse=True,
                                                         dtype=lh_sampling_dtype))
+                    random.pop_sseq()
                 else:
-                    np.random.seed(i+seed_offset)
+                    random.push_sseq(sseq[i])
                     _samples.append(met.draw_sample(from_inverse=True,
                                                     dtype=lh_sampling_dtype))
-            np.random.set_state(rand_state)
+                    random.pop_sseq()
             _samples = tuple(_samples)
             if mirror_samples:
                 n_samples *= 2
         self._samples = _samples
-        self._seed_offset = seed_offset
         self._n_samples = n_samples
         self._lin = Linearization.make_partial_var(mean, constants)
         v, g = None, None
@@ -205,7 +201,7 @@ class MetricGaussianKL_MPI(Energy):
         return MetricGaussianKL_MPI(
             position, self._hamiltonian, self._n_samples, self._constants,
             self._point_estimates, _samples=self._samples,
-            seed_offset=self._seed_offset, lh_sampling_dtype=self._sampdt)
+            lh_sampling_dtype=self._sampdt)
 
     @property
     def value(self):
@@ -245,11 +241,13 @@ class MetricGaussianKL_MPI(Energy):
             raise NotImplementedError()
         lin = self._lin.with_want_metric()
         samp = full(self._hamiltonian.domain, 0.)
-        rand_state = np.random.get_state()
-        np.random.seed(rank+np.random.randint(99999))
-        for v in self._samples:
+        sseq = random.spawn_sseq(n_samples)
+        for i, v in enumerate(self._samples):
+            # FIXME: this is not yet correct. We need to use the _global_ sample
+            # index, not the local one!
+            random.push_sseq(sseq[i])
             samp = samp + self._hamiltonian(lin+v).metric.draw_sample(from_inverse=False, dtype=dtype)
-        np.random.set_state(rand_state)
+            random.pop_sseq()
         return allreduce_sum_field(samp)
 
     def metric_sample(self, from_inverse=False, dtype=np.float64):
diff --git a/nifty6/random.py b/nifty6/random.py
index 673bf6596444d826f149520a26759f8dd423a603..53dd0faad1e8b9ed699daed867e13ab53e0a3c26 100644
--- a/nifty6/random.py
+++ b/nifty6/random.py
@@ -11,21 +11,50 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
-# Copyright(C) 2013-2019 Max-Planck-Society
+# Copyright(C) 2013-2020 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
 import numpy as np
 
+_sseq = [np.random.SeedSequence(42)]
+_rng = [np.random.default_rng(_sseq[-1])]
+
+
+def spawn_sseq(n, parent=None):
+    if parent is None:
+        global _sseq
+        parent = _sseq[-1]
+    return parent.spawn(n)
+
+
+def current_rng():
+    return _rng[-1]
+
+
+def push_sseq(sseq):
+    _sseq.append(sseq)
+    _rng.append(np.random.default_rng(_sseq[-1]))
+
+
+def push_sseq_from_seed(seed):
+    _sseq.append(np.random.SeedSequence(seed))
+    _rng.append(np.random.default_rng(_sseq[-1]))
+
+
+def pop_sseq():
+    _sseq.pop()
+    _rng.pop()
+
 
 class Random(object):
     @staticmethod
     def pm1(dtype, shape):
         if np.issubdtype(dtype, np.complexfloating):
             x = np.array([1+0j, 0+1j, -1+0j, 0-1j], dtype=dtype)
-            x = x[np.random.randint(4, size=shape)]
+            x = x[_rng[-1].integers(0, 4, size=shape)]
         else:
-            x = 2*np.random.randint(2, size=shape) - 1
+            x = 2*_rng[-1].integers(0, 2, size=shape)-1
         return x.astype(dtype, copy=False)
 
     @staticmethod
@@ -42,10 +71,10 @@ class Random(object):
             raise TypeError("mean must not be complex for a real result field")
         if np.issubdtype(dtype, np.complexfloating):
             x = np.empty(shape, dtype=dtype)
-            x.real = np.random.normal(mean.real, std*np.sqrt(0.5), shape)
-            x.imag = np.random.normal(mean.imag, std*np.sqrt(0.5), shape)
+            x.real = _rng[-1].normal(mean.real, std*np.sqrt(0.5), shape)
+            x.imag = _rng[-1].normal(mean.imag, std*np.sqrt(0.5), shape)
         else:
-            x = np.random.normal(mean, std, shape).astype(dtype, copy=False)
+            x = _rng[-1].normal(mean, std, shape).astype(dtype, copy=False)
         return x
 
     @staticmethod
@@ -57,13 +86,13 @@ class Random(object):
             raise TypeError("low and high must not be complex")
         if np.issubdtype(dtype, np.complexfloating):
             x = np.empty(shape, dtype=dtype)
-            x.real = np.random.uniform(low, high, shape)
-            x.imag = np.random.uniform(low, high, shape)
+            x.real = _rng[-1].uniform(low, high, shape)
+            x.imag = _rng[-1].uniform(low, high, shape)
         elif np.issubdtype(dtype, np.integer):
             if not (np.issubdtype(type(low), np.integer) and
                     np.issubdtype(type(high), np.integer)):
                 raise TypeError("low and high must be integer")
-            x = np.random.randint(low, high+1, shape)
+            x = _rng[-1].integers(low, high+1, shape)
         else:
-            x = np.random.uniform(low, high, shape)
+            x = _rng[-1].uniform(low, high, shape)
         return x.astype(dtype, copy=False)
diff --git a/setup.py b/setup.py
index cd9faeb37b1a8145049de86830d0fbbb9dd8c2f5..159c73e8e14224c64672b649ef7ae5fdde5eed40 100644
--- a/setup.py
+++ b/setup.py
@@ -39,8 +39,8 @@ setup(name="nifty6",
       packages=find_packages(include=["nifty6", "nifty6.*"]),
       zip_safe=True,
       license="GPLv3",
-      setup_requires=['scipy>=1.4.1'],
-      install_requires=['scipy>=1.4.1'],
+      setup_requires=['scipy>=1.4.1', 'numpy>=1.17'],
+      install_requires=['scipy>=1.4.1', 'numpy>=1.17'],
       python_requires='>=3.5',
       classifiers=[
         "Development Status :: 4 - Beta",
diff --git a/test/test_energy_gradients.py b/test/test_energy_gradients.py
index 190162fd41de1f8c626113a85d048cedcd117413..2c8d53d3285c2a811fadaae45bdd5e10f7d57773 100644
--- a/test/test_energy_gradients.py
+++ b/test/test_energy_gradients.py
@@ -37,9 +37,11 @@ pmp = pytest.mark.parametrize
 
 @pytest.fixture(params=PARAMS)
 def field(request):
-    np.random.seed(request.param[0])
+    ift.random.push_sseq_from_seed(request.param[0])
     S = ift.ScalingOperator(request.param[1], 1.)
-    return S.draw_sample()
+    res = S.draw_sample()
+    ift.random.pop_sseq()
+    return res
 
 
 def test_gaussian(field):
@@ -102,7 +104,7 @@ def test_inverse_gamma(field):
         return
     field = field.exp()
     space = field.domain
-    d = np.random.normal(10, size=space.shape)**2
+    d = ift.random.current_rng().normal(10, size=space.shape)**2
     d = ift.Field(space, d)
     energy = ift.InverseGammaLikelihood(d)
     ift.extra.check_jacobian_consistency(energy, field, tol=1e-5)
@@ -113,7 +115,7 @@ def testPoissonian(field):
         return
     field = field.exp()
     space = field.domain
-    d = np.random.poisson(120, size=space.shape)
+    d = ift.random.current_rng().poisson(120, size=space.shape)
     d = ift.Field(space, d)
     energy = ift.PoissonianEnergy(d)
     ift.extra.check_jacobian_consistency(energy, field, tol=1e-7)
@@ -124,7 +126,7 @@ def test_bernoulli(field):
         return
     field = field.sigmoid()
     space = field.domain
-    d = np.random.binomial(1, 0.1, size=space.shape)
+    d = ift.random.current_rng().binomial(1, 0.1, size=space.shape)
     d = ift.Field(space, d)
     energy = ift.BernoulliEnergy(d)
     ift.extra.check_jacobian_consistency(energy, field, tol=1e-5)
diff --git a/test/test_field.py b/test/test_field.py
index 0701d0c6da8bcb2ed3c8dafd2be9df8d167e80c0..03958c54ce3f747b6eeb78d9dfd4524083533c76 100644
--- a/test/test_field.py
+++ b/test/test_field.py
@@ -51,8 +51,6 @@ def _spec2(k):
 ])
 @pmp('space2', [ift.RGSpace((8,), harmonic=True), ift.LMSpace(12)])
 def test_power_synthesize_analyze(space1, space2):
-    np.random.seed(11)
-
     p1 = ift.PowerSpace(space1)
     fp1 = ift.PS_field(p1, _spec1)
     p2 = ift.PowerSpace(space2)
@@ -82,8 +80,6 @@ def test_power_synthesize_analyze(space1, space2):
 ])
 @pmp('space2', [ift.RGSpace((8,), harmonic=True), ift.LMSpace(12)])
 def test_DiagonalOperator_power_analyze2(space1, space2):
-    np.random.seed(11)
-
     fp1 = ift.PS_field(ift.PowerSpace(space1), _spec1)
     fp2 = ift.PS_field(ift.PowerSpace(space2), _spec2)
 
diff --git a/test/test_gaussian_energy.py b/test/test_gaussian_energy.py
index f2a020e83a551b0f3de2616e935aea9734c3bfd6..d7682b07f6627c03182d9241818d4b346e557e44 100644
--- a/test/test_gaussian_energy.py
+++ b/test/test_gaussian_energy.py
@@ -37,7 +37,7 @@ pmp = pytest.mark.parametrize
 @pmp('noise', [1, 1e-2, 1e2])
 @pmp('seed', [4, 78, 23])
 def test_gaussian_energy(space, nonlinearity, noise, seed):
-    np.random.seed(seed)
+    ift.random.push_sseq_from_seed(seed)
     dim = len(space.shape)
     hspace = space.get_default_codomain()
     ht = ift.HarmonicTransformOperator(hspace, target=space)
@@ -45,6 +45,7 @@ def test_gaussian_energy(space, nonlinearity, noise, seed):
     pspace = ift.PowerSpace(hspace, binbounds=binbounds)
     Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
     xi0 = ift.Field.from_random(domain=hspace, random_type='normal')
+    ift.random.pop_sseq()
 
     def pspec(k):
         return 1/(1 + k**2)**dim
diff --git a/test/test_kl.py b/test/test_kl.py
index 31ff187558ca7562aeb119626849f2c310a6b257..78f6a11e156c3f4ffcd6f431927b61a9b9193aeb 100644
--- a/test/test_kl.py
+++ b/test/test_kl.py
@@ -28,7 +28,6 @@ pmp = pytest.mark.parametrize
 @pmp('point_estimates', ([], ['a'], ['b'], ['a', 'b']))
 @pmp('mirror_samples', (True, False))
 def test_kl(constants, point_estimates, mirror_samples):
-    np.random.seed(42)
     dom = ift.RGSpace((12,), (2.12))
     op0 = ift.HarmonicSmoothingOperator(dom, 3)
     op = ift.ducktape(dom, None, 'a')*(op0.ducktape('b'))
diff --git a/test/test_minimizers.py b/test/test_minimizers.py
index e0c6641808cea8db0a95ba936ab8cb3fa1bdcfc0..4aec6fdff7f5dfd14120b742b3473bc7e1521fb9 100644
--- a/test/test_minimizers.py
+++ b/test/test_minimizers.py
@@ -51,7 +51,6 @@ slow_minimizers = ['ift.SteepestDescent(IC)']
      slow_minimizers)
 @pmp('space', spaces)
 def test_quadratic_minimization(minimizer, space):
-    np.random.seed(42)
     starting_point = ift.Field.from_random('normal', domain=space)*10
     covariance_diagonal = ift.Field.from_random('uniform', domain=space) + 0.5
     covariance = ift.DiagonalOperator(covariance_diagonal)
@@ -76,7 +75,6 @@ def test_quadratic_minimization(minimizer, space):
 
 @pmp('space', spaces)
 def test_WF_curvature(space):
-    np.random.seed(42)
     required_result = ift.full(space, 1.)
 
     s = ift.Field.from_random('uniform', domain=space) + 0.5
@@ -121,7 +119,6 @@ def test_rosenbrock(minimizer):
         from scipy.optimize import rosen, rosen_der, rosen_hess_prod
     except ImportError:
         raise SkipTest
-    np.random.seed(42)
     space = ift.DomainTuple.make(ift.UnstructuredDomain((2,)))
     starting_point = ift.Field.from_random('normal', domain=space)*10
 
diff --git a/test/test_operators/test_adjoint.py b/test/test_operators/test_adjoint.py
index c55f6b986c3f76fd250a8134ab15a2e853b95ec7..dcde9aed67c32695011bbc61878fe11f38025f8e 100644
--- a/test/test_operators/test_adjoint.py
+++ b/test/test_operators/test_adjoint.py
@@ -41,15 +41,13 @@ _pow_spaces = [
 pmp = pytest.mark.parametrize
 dtype = list2fixture([np.float64, np.complex128])
 
-np.random.seed(42)
-
 
 @pmp('sp', _p_RG_spaces)
 def testLOSResponse(sp, dtype):
-    starts = np.random.randn(len(sp.shape), 10)
-    ends = np.random.randn(len(sp.shape), 10)
-    sigma_low = 1e-4*np.random.randn(10)
-    sigma_ups = 1e-5*np.random.randn(10)
+    starts = ift.random.current_rng().standard_normal((len(sp.shape), 10))
+    ends = ift.random.current_rng().standard_normal((len(sp.shape), 10))
+    sigma_low = 1e-4*ift.random.current_rng().standard_normal(10)
+    sigma_ups = 1e-5*ift.random.current_rng().standard_normal(10)
     op = ift.LOSResponse(sp, starts, ends, sigma_low, sigma_ups)
     ift.extra.consistency_check(op, dtype, dtype)
 
@@ -70,7 +68,7 @@ def testOperatorCombinations(sp, dtype):
 
 def testLinearInterpolator():
     sp = ift.RGSpace((10, 8), distances=(0.1, 3.5))
-    pos = np.random.rand(2, 23)
+    pos = ift.random.current_rng().random((2, 23))
     pos[0, :] *= 0.9
     pos[1, :] *= 7*3.5
     op = ift.LinearInterpolator(sp, pos)
@@ -250,15 +248,16 @@ def testOuter(fdomain, domain):
 @pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
 @pmp('seed', [12, 3])
 def testValueInserter(sp, seed):
-    np.random.seed(seed)
+    ift.random.push_sseq_from_seed(seed)
     ind = []
     for ss in sp.shape:
         if ss == 1:
             ind.append(0)
         else:
-            ind.append(np.random.randint(0, ss-1))
+            ind.append(int(ift.random.current_rng().integers(0, ss-1)))
     op = ift.ValueInserter(sp, ind)
     ift.extra.consistency_check(op)
+    ift.random.pop_sseq()
 
 
 @pmp('sp', _pow_spaces)
@@ -282,18 +281,19 @@ def testSpecialSum(sp):
 @pmp('sp', [ift.RGSpace(10)])
 @pmp('seed', [12, 3])
 def testMatrixProductOperator(sp, seed):
-    np.random.seed(seed)
-    mat = np.random.randn(*sp.shape, *sp.shape)
+    ift.random.push_sseq_from_seed(seed)
+    mat = ift.random.current_rng().standard_normal((*sp.shape, *sp.shape))
     op = ift.MatrixProductOperator(sp, mat)
     ift.extra.consistency_check(op)
-    mat = mat + 1j*np.random.randn(*sp.shape, *sp.shape)
+    mat = mat + 1j*ift.random.current_rng().standard_normal((*sp.shape, *sp.shape))
     op = ift.MatrixProductOperator(sp, mat)
     ift.extra.consistency_check(op)
+    ift.random.pop_sseq()
 
 
 @pmp('seed', [12, 3])
 def testPartialExtractor(seed):
-    np.random.seed(seed)
+    ift.random.push_sseq_from_seed(seed)
     tgt = {'a': ift.RGSpace(1), 'b': ift.RGSpace(2)}
     dom = tgt.copy()
     dom['c'] = ift.RGSpace(3)
@@ -301,6 +301,7 @@ def testPartialExtractor(seed):
     tgt = ift.MultiDomain.make(tgt)
     op = ift.PartialExtractor(dom, tgt)
     ift.extra.consistency_check(op)
+    ift.random.pop_sseq()
 
 
 @pmp('seed', [12, 3])
diff --git a/test/test_operators/test_fft_operator.py b/test/test_operators/test_fft_operator.py
index 0fc4ccaf699ef1d246f97e6d413ec8e896297734..d8249e6a3ded749d8aec1c4151a072b00bec732e 100644
--- a/test/test_operators/test_fft_operator.py
+++ b/test/test_operators/test_fft_operator.py
@@ -42,7 +42,6 @@ def test_fft1D(d, dtype, op):
     tol = _get_rtol(dtype)
     a = ift.RGSpace(dim1, distances=d)
     b = ift.RGSpace(dim1, distances=1./(dim1*d), harmonic=True)
-    np.random.seed(16)
 
     fft = op(domain=a, target=b)
     inp = ift.Field.from_random(
diff --git a/test/test_operators/test_jacobian.py b/test/test_operators/test_jacobian.py
index 39cca6de94004dc2dc0997ddf825c401091aa03e..bf532d93ba9942be15630889c19eefc57a31b12a 100644
--- a/test/test_operators/test_jacobian.py
+++ b/test/test_operators/test_jacobian.py
@@ -38,18 +38,19 @@ seed = list2fixture([4, 78, 23])
 
 
 def testBasics(space, seed):
-    np.random.seed(seed)
+    ift.random.push_sseq_from_seed(seed)
     S = ift.ScalingOperator(space, 1.)
     s = S.draw_sample()
     var = ift.Linearization.make_var(s)
     model = ift.ScalingOperator(var.target, 6.)
     ift.extra.check_jacobian_consistency(model, var.val)
+    ift.random.pop_sseq()
 
 
 @pmp('type1', ['Variable', 'Constant'])
 @pmp('type2', ['Variable'])
 def testBinary(type1, type2, space, seed):
-    np.random.seed(seed)
+    ift.random.push_sseq_from_seed(seed)
     dom1 = ift.MultiDomain.make({'s1': space})
     dom2 = ift.MultiDomain.make({'s2': space})
     dom = ift.MultiDomain.union((dom1, dom2))
@@ -87,10 +88,11 @@ def testBinary(type1, type2, space, seed):
         model = ift.FFTOperator(space)(select_s1*select_s2)
         pos = ift.from_random("normal", dom)
         ift.extra.check_jacobian_consistency(model, pos, ntries=20)
+    ift.random.pop_sseq()
 
 
 def testSpecialDistributionOps(space, seed):
-    np.random.seed(seed)
+    ift.random.push_sseq_from_seed(seed)
     S = ift.ScalingOperator(space, 1.)
     pos = S.draw_sample()
     alpha = 1.5
@@ -99,11 +101,12 @@ def testSpecialDistributionOps(space, seed):
     ift.extra.check_jacobian_consistency(model, pos, ntries=20)
     model = ift.UniformOperator(space, alpha, q)
     ift.extra.check_jacobian_consistency(model, pos, ntries=20)
+    ift.random.pop_sseq()
 
 
 @pmp('neg', [True, False])
 def testAdder(space, seed, neg):
-    np.random.seed(seed)
+    ift.random.push_sseq_from_seed(seed)
     S = ift.ScalingOperator(space, 1.)
     f = S.draw_sample()
     f1 = S.draw_sample()
@@ -111,6 +114,7 @@ def testAdder(space, seed, neg):
     ift.extra.check_jacobian_consistency(op, f)
     op = ift.Adder(f1.val.ravel()[0], neg=neg, domain=space)
     ift.extra.check_jacobian_consistency(op, f)
+    ift.random.pop_sseq()
 
 
 @pmp('target', [ift.RGSpace(64, distances=.789, harmonic=True),
@@ -119,7 +123,7 @@ def testAdder(space, seed, neg):
 @pmp('causal', [True, False])
 @pmp('minimum_phase', [True, False])
 def testDynamicModel(target, causal, minimum_phase, seed):
-    np.random.seed(seed)
+    ift.random.push_sseq_from_seed(seed)
     dct = {
             'target': target,
             'harmonic_padding': None,
@@ -156,6 +160,7 @@ def testDynamicModel(target, causal, minimum_phase, seed):
         # FIXME I dont know why smaller tol fails for 3D example
         ift.extra.check_jacobian_consistency(
             model, pos, tol=1e-5, ntries=20)
+    ift.random.pop_sseq()
 
 
 @pmp('h_space', _h_spaces)
diff --git a/test/test_operators/test_nft.py b/test/test_operators/test_nft.py
index dcce61bf917795cab321a61f93b67875a9f4d1c8..bd6012dc405b213904519e63adfff395749ce18f 100644
--- a/test/test_operators/test_nft.py
+++ b/test/test_operators/test_nft.py
@@ -21,8 +21,6 @@ from numpy.testing import assert_
 
 import nifty6 as ift
 
-np.random.seed(40)
-
 pmp = pytest.mark.parametrize
 
 
@@ -35,8 +33,9 @@ def _l2error(a, b):
 @pmp('nv', [4, 12, 128])
 @pmp('N', [1, 10, 100])
 def test_gridding(nu, nv, N, eps):
-    uv = np.random.rand(N, 2) - 0.5
-    vis = np.random.randn(N) + 1j*np.random.randn(N)
+    uv = ift.random.current_rng().random((N, 2)) - 0.5
+    vis = (ift.random.current_rng().standard_normal(N)
+           + 1j*ift.random.current_rng().standard_normal(N))
 
     # Nifty
     dom = ift.RGSpace((nu, nv), distances=(0.2, 1.12))
@@ -92,7 +91,7 @@ def test_cartesian():
 @pmp('N', [1, 10, 100])
 def test_build(nu, nv, N, eps):
     dom = ift.RGSpace([nu, nv])
-    uv = np.random.rand(N, 2) - 0.5
+    uv = ift.random.current_rng().random((N, 2)) - 0.5
     GM = ift.GridderMaker(dom, uv=uv, eps=eps)
     R0 = GM.getGridder()
     R1 = GM.getRest()
diff --git a/test/test_operators/test_representation.py b/test/test_operators/test_representation.py
index b817b18be6e5c60f27cc4391ec92cf57ae8cca4d..b6d2f1b3c4b881c1464590f554372f09901ebdfa 100644
--- a/test/test_operators/test_representation.py
+++ b/test/test_operators/test_representation.py
@@ -11,7 +11,7 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
-# Copyright(C) 2013-2019 Max-Planck-Society
+# Copyright(C) 2013-2020 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
@@ -44,10 +44,10 @@ def _check_repr(op):
 
 @pmp('sp', _p_RG_spaces)
 def testLOSResponse(sp, dtype):
-    starts = np.random.randn(len(sp.shape), 10)
-    ends = np.random.randn(len(sp.shape), 10)
-    sigma_low = 1e-4*np.random.randn(10)
-    sigma_ups = 1e-5*np.random.randn(10)
+    starts = ift.random.current_rng().standard_normal((len(sp.shape), 10))
+    ends = ift.random.current_rng().standard_normal((len(sp.shape), 10))
+    sigma_low = 1e-4*ift.random.current_rng().standard_normal(10)
+    sigma_ups = 1e-5*ift.random.current_rng().standard_normal(10)
     _check_repr(ift.LOSResponse(sp, starts, ends, sigma_low, sigma_ups))
 
 
@@ -65,7 +65,7 @@ def testOperatorCombinations(sp, dtype):
 
 def testLinearInterpolator():
     sp = ift.RGSpace((10, 8), distances=(0.1, 3.5))
-    pos = np.random.rand(2, 23)
+    pos = ift.random.current_rng().random((2, 23))
     pos[0, :] *= 0.9
     pos[1, :] *= 7*3.5
     _check_repr(ift.LinearInterpolator(sp, pos))
@@ -206,11 +206,12 @@ def testOuter(fdomain, domain):
 @pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
 @pmp('seed', [12, 3])
 def testValueInserter(sp, seed):
-    np.random.seed(seed)
+    ift.random.push_sseq_from_seed(seed)
     ind = []
     for ss in sp.shape:
         if ss == 1:
             ind.append(0)
         else:
-            ind.append(np.random.randint(0, ss - 1))
+            ind.append(int(ift.random.current_rng().integers(0, ss - 1)))
     _check_repr(ift.ValueInserter(sp, ind))
+    ift.random.pop_sseq()
diff --git a/test/test_operators/test_value_inserter.py b/test/test_operators/test_value_inserter.py
index d420585b357e0bf79ab1440125b24c9ebea77160..1749bf906d71419cdb0f6a81d2bec7ae4c48c561 100644
--- a/test/test_operators/test_value_inserter.py
+++ b/test/test_operators/test_value_inserter.py
@@ -31,10 +31,11 @@ import nifty6 as ift
 ])
 @pytest.mark.parametrize('seed', [13, 2])
 def test_value_inserter(sp, seed):
-    np.random.seed(seed)
-    ind = tuple([np.random.randint(0, ss - 1) for ss in sp.shape])
+    ift.random.push_sseq_from_seed(seed)
+    ind = tuple([int(ift.random.current_rng().integers(0, ss - 1)) for ss in sp.shape])
     op = ift.ValueInserter(sp, ind)
     f = ift.from_random('normal', op.domain)
+    ift.random.pop_sseq()
     inp = f.val
     ret = op(f).val
     assert_(ret[ind] == inp)
diff --git a/test/test_plot.py b/test/test_plot.py
index 3f220c317d5cea55435bef52d3d77d479a29db96..1c8475fa78ac53bfbdff4199bf96ea0be444a7e0 100644
--- a/test/test_plot.py
+++ b/test/test_plot.py
@@ -29,12 +29,12 @@ def test_plots():
 
     fft = ift.FFTOperator(rg_space2)
 
-    field_rg1_1 = ift.Field(rg_space1, np.random.randn(100))
-    field_rg1_2 = ift.Field(rg_space1, np.random.randn(100))
+    field_rg1_1 = ift.Field(rg_space1, ift.random.current_rng().standard_normal(100))
+    field_rg1_2 = ift.Field(rg_space1, ift.random.current_rng().standard_normal(100))
     field_rg2 = ift.Field(
-        rg_space2, np.random.randn(80*60).reshape((80, 60)))
-    field_hp = ift.Field(hp_space, np.random.randn(12*64**2))
-    field_gl = ift.Field(gl_space, np.random.randn(32640))
+        rg_space2, ift.random.current_rng().standard_normal((80,60)))
+    field_hp = ift.Field(hp_space, ift.random.current_rng().standard_normal(12*64**2))
+    field_gl = ift.Field(gl_space, ift.random.current_rng().standard_normal(32640))
     field_ps = ift.power_analyze(fft.times(field_rg2))
 
     plot = ift.Plot()
diff --git a/test/test_spaces/test_gl_space.py b/test/test_spaces/test_gl_space.py
index 76a07cad1eb2c6078142b79d03fff90718e805e4..131bf6ba575574784e6d131ed098c27d6cc33164 100644
--- a/test/test_spaces/test_gl_space.py
+++ b/test/test_spaces/test_gl_space.py
@@ -41,7 +41,6 @@ CONSTRUCTOR_CONFIGS = [[
 
 
 def get_dvol_configs():
-    np.random.seed(42)
     wgt = [2.0943951, 2.0943951]
     # for GLSpace(nlat=2, nlon=3)
     dvol_0 = np.array(
diff --git a/test/test_spaces/test_rg_space.py b/test/test_spaces/test_rg_space.py
index 1cb535c6ba30ef51aab8a3a1722aca0e190ed1ed..b68b0184340243d4868a328f8a954572fbba8abd 100644
--- a/test/test_spaces/test_rg_space.py
+++ b/test/test_spaces/test_rg_space.py
@@ -66,7 +66,6 @@ def get_k_length_array_configs():
 
 
 def get_dvol_configs():
-    np.random.seed(42)
     return [[(11, 11), None, False, 1], [(11, 11), None, False, 1],
             [(11, 11), (1.3, 1.3), False, 1], [(11, 11), (1.3, 1.3), False,
                                                1], [(11, 11), None, True, 1],
diff --git a/test/test_sugar.py b/test/test_sugar.py
index a572c5113784d38885e68f6706a958c5542f4807..660720640ba4b8e0474be8d6761cd27966de7b74 100644
--- a/test/test_sugar.py
+++ b/test/test_sugar.py
@@ -52,7 +52,6 @@ def test_exec_time():
 
 
 def test_calc_pos():
-    np.random.seed(42)
     dom = ift.RGSpace(12, harmonic=True)
     op = ift.HarmonicTransformOperator(dom).exp()
     fld = op(0.1*ift.from_random('normal', op.domain))