Repository: sympy/sympy
Pull request: create core.random
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->

Fixes #12628
Fixes #14589
closes #16967 as alternate

#### Brief description of what is fixed or changed

locate random-related items to the core and provide needed random generators that all run from a single instance of python's `random.Random()`
#### Other comments

#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
NO ENTRY
<!-- END RELEASE NOTES -->


Problem statement:
as_numer_denom non-deterministically influences the state of Python’s random number generator
When I apply `as_numer_denom` to a power expression, this changes the state of Python’s random number generator. This happens in a non-deterministic fashion, i.e., it isn’t just that random numbers are generated for whatever reason, but the state seems to be reset using some external information (like the system time).

Here is a minimal non-working example:

``` Python

import random
from sympy.abc import x

results = set()
for _ in range(10):
	random.seed(42)
	(x**2).as_numer_denom()
	results.add(random.random())
assert len(results)==1
```

interaction of sympy and random
I noticed the following strange behaviour when generating a random function
`import sympy, random

x = sympy.symbols('x')

random.seed(10)

y = 0

for i in range(4): y += sympy.sin(random.uniform(-10,10) * x)

random.seed(10)

z = 0

for i in range(4): z += sympy.sin(random.uniform(-10,10) * x)
`
I get different values for y and z.

Hints from issue discussion:
We should be using a random.Random instance so that it doesn't affect the random seed for other applications. 

As a workaround, I suggest doing the same in your own code, so that nothing else can mess with random like this. 
> As a workaround, I suggest doing the same in your own code, so that nothing else can mess with random like this.

That’s what I eventually did. I stumbled upon this when checking something quickly and dirtily and since I didn’t expect SymPy to use random numbers for this (let alone in a way that is not even deterministic), I was rather baffled and it took me a while to find out what happened.
For what it's worth I just spent a long time discovering the same thing by process of elimination in my own code.
Same here. Spent about two hours until finding the issue.

I've written another short example because I found something that might be useful: the issue happens only at the first call of `my_test`. Subsequent calls do not affect python's random state.

```python
import sympy as sp
import random

def my_test():
    random.seed('HolaAmigo')
    # Should print -1735515812038708070 (may vary on your laptop)
    print(hash(random.getstate()))
    for i in range(100):
        x = random.randrange(2) * sp.Symbol('x')
    # Should print 3791330111857148225 (may vary on your laptop)
    print(hash(random.getstate()))
    return

my_test() # fails
my_test() # ok
my_test() # ok
my_test() # ok
```
You may comment out the sp.Symbol to check what the value would be if no intervention from sympy occurred.
I ran into a similar problem, using `sympy.acos`.
The following program shows different results for variable `after` for consecutive runs.

    import sympy
    import random
    random.seed(28)
    
    before = random.uniform(0,1)
    print('before: ',before)
    m0, m1= sympy.symbols('m_0 m_1', real=True)
    result = sympy.acos(-m0/m1)
    after = random.uniform(0,1)
    print('after: ',after)

What about adding a warning in the documentation about this behavior? I think the section `Gotchas and Pitfalls` would be good for this.
This is a bug. SymPy should be using its own independent random state wherever it is using random numbers. Of course, the simplest workaround is to also do the same in your own code (use `random.Random()` instead of `random`). That's recommended practice regardless, because you never know what other code will be calling the global random number generator.
The values I get do not look different:
```
>>> random.seed(10)
>>> y = 0
>>> for i in range(4): y += sympy.sin(random.uniform(-10,10) * x)
... 
>>> y
-sin(1.42221890649771*x) + sin(1.42805189379827*x) + sin(1.56182602268941*x) - sin(5.87803535720997*x)
>>> random.seed(10)
>>> z = 0
>>> for i in range(4): z += sympy.sin(random.uniform(-10,10) * x)
... 
>>> z
-sin(1.42221890649771*x) + sin(1.42805189379827*x) + sin(1.56182602268941*x) - sin(5.87803535720997*x)
>>> 
```
Here is a screenshot :
 


    Le lundi 2 avril 2018 à 13:21:38 UTC+2, Kalevi Suominen <notifications@github.com> a écrit :  
 
 
The values I get do not look different:
>>> random.seed(10)
>>> y = 0
>>> for i in range(4): y += sympy.sin(random.uniform(-10,10) * x)
... 
>>> y
-sin(1.42221890649771*x) + sin(1.42805189379827*x) + sin(1.56182602268941*x) - sin(5.87803535720997*x)
>>> random.seed(10)
>>> z = 0
>>> for i in range(4): z += sympy.sin(random.uniform(-10,10) * x)
... 
>>> z
-sin(1.42221890649771*x) + sin(1.42805189379827*x) + sin(1.56182602268941*x) - sin(5.87803535720997*x)
>>> 

—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub, or mute the thread.
    
I get

```
>>> y
sin(1.42805189379827*x) + sin(1.55778912915404*x) + sin(5.54243799539346*x) + sin(8.0835978129377*x)
>>> z
-sin(1.42221890649771*x) + sin(1.42805189379827*x) + sin(3.06291637379069*x) - sin(7.68145060132939*x)
```

and I get those same two values each time. 
Actually after a second run, I get the same values as @jksuom. So it's probably related to the cache somehow. 
@esandier you should be able to get consistent results by using `random.Random` instances, which are seeded independently (this is always a good practice anyway). 
Indeed this works. Evaluating the loop within "with evaluate(False)" works as well.
Similar issue: https://github.com/sympy/sympy/issues/12628
Why did you close this? I don't think it was fixed. 
Sorry,

I'm not familiar with how issues are handled, I thought I was supposed to close it after a while

Etienne
 

    Le mardi 1 mai 2018 à 21:33:31 UTC+2, Aaron Meurer <notifications@github.com> a écrit :  
 
 
Why did you close this? I don't think it was fixed.

—
You are receiving this because you modified the open/close state.
Reply to this email directly, view it on GitHub, or mute the thread.
    
@esandier : Bonjour, pourriez vous m'écrire ? Je me suis inscrit en follower sur votre projet pywims. Je souhaiterais l'intégrer dans un de mes sites.
Merci.

Commit message: add seed function to core.random

Patch to review:
diff --git a/sympy/benchmarks/bench_symbench.py b/sympy/benchmarks/bench_symbench.py
index 2ec37f7bb512..8ea700b44b67 100644
--- a/sympy/benchmarks/bench_symbench.py
+++ b/sympy/benchmarks/bench_symbench.py
@@ -1,5 +1,5 @@
 #!/usr/bin/env python
-from random import random
+from sympy.core.random import random
 from sympy.core.numbers import (I, Integer, pi)
 from sympy.core.symbol import Symbol
 from sympy.core.sympify import sympify
diff --git a/sympy/combinatorics/partitions.py b/sympy/combinatorics/partitions.py
index 2742c3cb0683..77266fc3e0e2 100644
--- a/sympy/combinatorics/partitions.py
+++ b/sympy/combinatorics/partitions.py
@@ -598,7 +598,7 @@ def random_integer_partition(n, seed=None):
     >>> random_integer_partition(1)
     [1]
     """
-    from sympy.testing.randtest import _randint
+    from sympy.core.random import _randint
 
     n = as_int(n)
     if n < 1:
diff --git a/sympy/combinatorics/perm_groups.py b/sympy/combinatorics/perm_groups.py
index c865b12eadc6..c0b91c85c9b3 100644
--- a/sympy/combinatorics/perm_groups.py
+++ b/sympy/combinatorics/perm_groups.py
@@ -1,4 +1,4 @@
-from random import randrange, choice
+from sympy.core.random import randrange, choice
 from math import log
 from sympy.ntheory import primefactors
 from sympy.core.symbol import Symbol
@@ -15,7 +15,7 @@
 from sympy.functions.combinatorial.factorials import factorial
 from sympy.ntheory import sieve
 from sympy.utilities.iterables import has_variety, is_sequence, uniq
-from sympy.testing.randtest import _randrange
+from sympy.core.random import _randrange
 from itertools import islice
 from sympy.core.sympify import _sympify
 rmul = Permutation.rmul_with_af
diff --git a/sympy/combinatorics/tests/test_schur_number.py b/sympy/combinatorics/tests/test_schur_number.py
index 97f41471147c..a2b9348122cd 100644
--- a/sympy/combinatorics/tests/test_schur_number.py
+++ b/sympy/combinatorics/tests/test_schur_number.py
@@ -1,6 +1,6 @@
 from sympy.core import S, Rational
 from sympy.combinatorics.schur_number import schur_partition, SchurNumber
-from sympy.testing.randtest import _randint
+from sympy.core.random import _randint
 from sympy.testing.pytest import raises
 from sympy.core.symbol import symbols
 
diff --git a/sympy/combinatorics/tests/test_testutil.py b/sympy/combinatorics/tests/test_testutil.py
index 393746ce0bd8..e13f4d5b9913 100644
--- a/sympy/combinatorics/tests/test_testutil.py
+++ b/sympy/combinatorics/tests/test_testutil.py
@@ -5,7 +5,7 @@
     _verify_normal_closure
 from sympy.combinatorics.permutations import Permutation
 from sympy.combinatorics.perm_groups import PermutationGroup
-from random import shuffle
+from sympy.core.random import shuffle
 
 
 def test_cmp_perm_lists():
diff --git a/sympy/core/assumptions.py b/sympy/core/assumptions.py
index cf1f5e04d261..11ac8fd9e451 100644
--- a/sympy/core/assumptions.py
+++ b/sympy/core/assumptions.py
@@ -214,7 +214,7 @@
 from .core import BasicMeta
 from .sympify import sympify
 
-from random import shuffle
+from sympy.core.random import shuffle
 
 
 _assume_rules = FactRules([
diff --git a/sympy/core/expr.py b/sympy/core/expr.py
index c9c9e8150c6d..2c71539ee6a9 100644
--- a/sympy/core/expr.py
+++ b/sympy/core/expr.py
@@ -491,13 +491,13 @@ def _random(self, n=None, re_min=-1, im_min=-1, re_max=1, im_max=1):
         See Also
         ========
 
-        sympy.testing.randtest.random_complex_number
+        sympy.core.random.random_complex_number
         """
 
         free = self.free_symbols
         prec = 1
         if free:
-            from sympy.testing.randtest import random_complex_number
+            from sympy.core.random import random_complex_number
             a, c, b, d = re_min, re_max, im_min, im_max
             reps = dict(list(zip(free, [random_complex_number(a, b, c, d, rational=True)
                            for zi in free])))
diff --git a/sympy/core/random.py b/sympy/core/random.py
new file mode 100644
index 000000000000..a3477b0d52cb
--- /dev/null
+++ b/sympy/core/random.py
@@ -0,0 +1,209 @@
+"""When you need to use random, import it from here so there is only
+one generator working for SymPy. Imports from here should behave the
+same as if they were being imported from Python's random module. But
+only the routines currently used in SymPy are included here. To use
+others import `rng` and access the method directly. For example, to
+capture the current state of the generator use `rng.getstate()`.
+
+There is intentionally no Random to import from here. If you want
+to control the state of the generator, import `seed` and call it
+with or without an argument to set the state.
+
+EXAMPLES
+========
+
+>>> from sympy.core.random import random, seed
+>>> assert random() < 1
+>>> seed(1); a = random()
+>>> b = random()
+>>> seed(1); c = random()
+>>> assert a == c
+>>> assert a != b  # remote possibility this will fail
+"""
+
+import random as _random
+rng = _random.Random()
+
+seed = lambda s=None: rng.seed(s)
+
+choice = rng.choice
+random = rng.random
+randint = rng.randint
+randrange = rng.randrange
+shuffle = rng.shuffle
+uniform = rng.uniform
+
+def random_complex_number(a=2, b=-1, c=3, d=1, rational=False, tolerance=None):
+    """
+    Return a random complex number.
+
+    To reduce chance of hitting branch cuts or anything, we guarantee
+    b <= Im z <= d, a <= Re z <= c
+
+    When rational is True, a rational approximation to a random number
+    is obtained within specified tolerance, if any.
+    """
+    from sympy.core.numbers import I
+    from sympy.simplify.simplify import nsimplify
+    A, B = uniform(a, c), uniform(b, d)
+    if not rational:
+        return A + I*B
+    return (nsimplify(A, rational=True, tolerance=tolerance) +
+        I*nsimplify(B, rational=True, tolerance=tolerance))
+
+
+def verify_numerically(f, g, z=None, tol=1.0e-6, a=2, b=-1, c=3, d=1):
+    """
+    Test numerically that f and g agree when evaluated in the argument z.
+
+    If z is None, all symbols will be tested. This routine does not test
+    whether there are Floats present with precision higher than 15 digits
+    so if there are, your results may not be what you expect due to round-
+    off errors.
+
+    Examples
+    ========
+
+    >>> from sympy import sin, cos
+    >>> from sympy.abc import x
+    >>> from sympy.core.random import verify_numerically as tn
+    >>> tn(sin(x)**2 + cos(x)**2, 1, x)
+    True
+    """
+    from sympy.core.symbol import Symbol
+    from sympy.core.numbers import comp
+    from sympy.core.containers import Tuple
+    f, g, z = Tuple(f, g, z)
+    z = [z] if isinstance(z, Symbol) else (f.free_symbols | g.free_symbols)
+    reps = list(zip(z, [random_complex_number(a, b, c, d) for _ in z]))
+    z1 = f.subs(reps).n()
+    z2 = g.subs(reps).n()
+    return comp(z1, z2, tol)
+
+
+def test_derivative_numerically(f, z, tol=1.0e-6, a=2, b=-1, c=3, d=1):
+    """
+    Test numerically that the symbolically computed derivative of f
+    with respect to z is correct.
+
+    This routine does not test whether there are Floats present with
+    precision higher than 15 digits so if there are, your results may
+    not be what you expect due to round-off errors.
+
+    Examples
+    ========
+
+    >>> from sympy import sin
+    >>> from sympy.abc import x
+    >>> from sympy.core.random import test_derivative_numerically as td
+    >>> td(sin(x), x)
+    True
+    """
+    from sympy.core.numbers import comp
+    from sympy.core.function import Derivative
+    z0 = random_complex_number(a, b, c, d)
+    f1 = f.diff(z).subs(z, z0)
+    f2 = Derivative(f, z).doit_numerically(z0)
+    return comp(f1.n(), f2.n(), tol)
+
+
+def _randrange(seed=None):
+    """Return a randrange generator. ``seed`` can be
+        o None - return randomly seeded generator
+        o int - return a generator seeded with the int
+        o list - the values to be returned will be taken from the list
+          in the order given; the provided list is not modified.
+
+    Examples
+    ========
+
+    >>> from sympy.core.random import _randrange
+    >>> rr = _randrange()
+    >>> rr(1000) # doctest: +SKIP
+    999
+    >>> rr = _randrange(3)
+    >>> rr(1000) # doctest: +SKIP
+    238
+    >>> rr = _randrange([0, 5, 1, 3, 4])
+    >>> rr(3), rr(3)
+    (0, 1)
+    """
+    if seed is None:
+        return random.randrange
+    elif isinstance(seed, int):
+        random.seed(seed)
+        return random.randrange
+    elif is_sequence(seed):
+        seed = list(seed)  # make a copy
+        seed.reverse()
+
+        def give(a, b=None, seq=seed):
+            if b is None:
+                a, b = 0, a
+            a, b = as_int(a), as_int(b)
+            w = b - a
+            if w < 1:
+                raise ValueError('_randrange got empty range')
+            try:
+                x = seq.pop()
+            except IndexError:
+                raise ValueError('_randrange sequence was too short')
+            if a <= x < b:
+                return x
+            else:
+                return give(a, b, seq)
+        return give
+    else:
+        raise ValueError('_randrange got an unexpected seed')
+
+
+def _randint(seed=None):
+    """Return a randint generator. ``seed`` can be
+        o None - return randomly seeded generator
+        o int - return a generator seeded with the int
+        o list - the values to be returned will be taken from the list
+          in the order given; the provided list is not modified.
+
+    Examples
+    ========
+
+    >>> from sympy.core.random import _randint
+    >>> ri = _randint()
+    >>> ri(1, 1000) # doctest: +SKIP
+    999
+    >>> ri = _randint(3)
+    >>> ri(1, 1000) # doctest: +SKIP
+    238
+    >>> ri = _randint([0, 5, 1, 2, 4])
+    >>> ri(1, 3), ri(1, 3)
+    (1, 2)
+    """
+    if seed is None:
+        return random.randint
+    elif isinstance(seed, int):
+        random.seed(seed)
+        return random.randint
+    elif is_sequence(seed):
+        seed = list(seed)  # make a copy
+        seed.reverse()
+
+        def give(a, b, seq=seed):
+            a, b = as_int(a), as_int(b)
+            w = b - a
+            if w < 0:
+                raise ValueError('_randint got empty range')
+            try:
+                x = seq.pop()
+            except IndexError:
+                raise ValueError('_randint sequence was too short')
+            if a <= x <= b:
+                return x
+            else:
+                return give(a, b, seq)
+        return give
+    else:
+        raise ValueError('_randint got an unexpected seed')
+
+
+from sympy.utilities.iterables import is_sequence
+from sympy.utilities.misc import as_int
diff --git a/sympy/core/sympify.py b/sympy/core/sympify.py
index ae1f95bff9a2..86deba6bf541 100644
--- a/sympy/core/sympify.py
+++ b/sympy/core/sympify.py
@@ -6,7 +6,7 @@
 
 from inspect import getmro
 import string
-from random import choice
+from sympy.core.random import choice
 
 from .parameters import global_parameters
 
diff --git a/sympy/core/tests/test_arit.py b/sympy/core/tests/test_arit.py
index 1cc63610c734..ac95ead473e3 100644
--- a/sympy/core/tests/test_arit.py
+++ b/sympy/core/tests/test_arit.py
@@ -21,7 +21,7 @@
 from sympy.core.expr import unchanged
 from sympy.utilities.iterables import cartes
 from sympy.testing.pytest import XFAIL, raises, warns_deprecated_sympy
-from sympy.testing.randtest import verify_numerically
+from sympy.core.random import verify_numerically
 from sympy.functions.elementary.trigonometric import asin
 
 
diff --git a/sympy/core/tests/test_expand.py b/sympy/core/tests/test_expand.py
index d041b209048d..25029d4d8e72 100644
--- a/sympy/core/tests/test_expand.py
+++ b/sympy/core/tests/test_expand.py
@@ -11,7 +11,7 @@
 from sympy.core.function import expand, expand_multinomial, expand_power_base
 
 from sympy.testing.pytest import raises
-from sympy.testing.randtest import verify_numerically
+from sympy.core.random import verify_numerically
 
 from sympy.abc import x, y, z
 
diff --git a/sympy/core/tests/test_function.py b/sympy/core/tests/test_function.py
index b83806537e5b..ec89f9d2e127 100644
--- a/sympy/core/tests/test_function.py
+++ b/sympy/core/tests/test_function.py
@@ -623,8 +623,7 @@ def ok(a):
 
 
 def test_derivative_numerically():
-    from random import random
-    z0 = random() + I*random()
+    z0 = x._random()
     assert abs(Derivative(sin(x), x).doit_numerically(z0) - cos(z0)) < 1e-15
 
 
diff --git a/sympy/core/tests/test_random.py b/sympy/core/tests/test_random.py
new file mode 100644
index 000000000000..16b768e49953
--- /dev/null
+++ b/sympy/core/tests/test_random.py
@@ -0,0 +1,30 @@
+import random
+from sympy.core.symbol import Symbol, symbols
+from sympy.functions.elementary.trigonometric import sin, acos
+from sympy.abc import x
+
+def test_random():
+    random.seed(42)
+    a = random.random()
+    random.seed(42)
+    Symbol('z').is_finite
+    b = random.random()
+    assert a == b
+
+    got = set()
+    for i in range(2):
+        random.seed(28)
+        m0, m1 = symbols('m_0 m_1', real=True)
+        _ = acos(-m0/m1)
+        got.add(random.uniform(0,1))
+    assert len(got) == 1
+
+    random.seed(10)
+    y = 0
+    for i in range(4):
+        y += sin(random.uniform(-10,10) * x)
+    random.seed(10)
+    z = 0
+    for i in range(4):
+        z += sin(random.uniform(-10,10) * x)
+    assert y == z
diff --git a/sympy/core/tests/test_relational.py b/sympy/core/tests/test_relational.py
index eeb5f9afec4d..a28da387c0f3 100644
--- a/sympy/core/tests/test_relational.py
+++ b/sympy/core/tests/test_relational.py
@@ -379,7 +379,7 @@ def test_new_relational():
     assert (x < 0) != StrictLessThan(x, 1)
 
     # finally, some fuzz testing
-    from random import randint
+    from sympy.core.random import randint
     for i in range(100):
         while 1:
             strtype, length = (chr, 65535) if randint(0, 1) else (chr, 255)
diff --git a/sympy/crypto/crypto.py b/sympy/crypto/crypto.py
index 821398d0aaf9..05d5c1dcefeb 100644
--- a/sympy/crypto/crypto.py
+++ b/sympy/crypto/crypto.py
@@ -27,7 +27,7 @@
 from sympy.polys.polytools import gcd, Poly
 from sympy.utilities.misc import as_int, filldedent, translate
 from sympy.utilities.iterables import uniq, multiset
-from sympy.testing.randtest import _randrange, _randint
+from sympy.core.random import _randrange, _randint
 
 
 class NonInvertibleCipherWarning(RuntimeWarning):
@@ -2571,7 +2571,7 @@ def elgamal_private_key(digit=10, seed=None):
     =====
 
     For testing purposes, the ``seed`` parameter may be set to control
-    the output of this routine. See sympy.testing.randtest._randrange.
+    the output of this routine. See sympy.core.random._randrange.
 
     Examples
     ========
@@ -2656,7 +2656,7 @@ def encipher_elgamal(i, key, seed=None):
     =====
 
     For testing purposes, the ``seed`` parameter may be set to control
-    the output of this routine. See sympy.testing.randtest._randrange.
+    the output of this routine. See sympy.core.random._randrange.
 
     Examples
     ========
@@ -2764,7 +2764,7 @@ def dh_private_key(digit=10, seed=None):
     =====
 
     For testing purposes, the ``seed`` parameter may be set to control
-    the output of this routine. See sympy.testing.randtest._randrange.
+    the output of this routine. See sympy.core.random._randrange.
 
     Examples
     ========
diff --git a/sympy/crypto/tests/test_crypto.py b/sympy/crypto/tests/test_crypto.py
index a38af2297dd2..36067bfd0181 100644
--- a/sympy/crypto/tests/test_crypto.py
+++ b/sympy/crypto/tests/test_crypto.py
@@ -22,7 +22,7 @@
 
 from sympy.testing.pytest import raises, warns
 
-from random import randrange
+from sympy.core.random import randrange
 
 def test_encipher_railfence():
     assert encipher_railfence("hello world",2) == "hlowrdel ol"
diff --git a/sympy/functions/combinatorial/tests/test_comb_numbers.py b/sympy/functions/combinatorial/tests/test_comb_numbers.py
index 607bb0362a57..64a30a48c1dc 100644
--- a/sympy/functions/combinatorial/tests/test_comb_numbers.py
+++ b/sympy/functions/combinatorial/tests/test_comb_numbers.py
@@ -492,7 +492,7 @@ def test_nC_nP_nT():
         nP, nC, nT, stirling, _stirling1, _stirling2, _multiset_histogram, _AOP_product)
 
     from sympy.combinatorics.permutations import Permutation
-    from random import choice
+    from sympy.core.random import choice
 
     c = string.ascii_lowercase
     for i in range(100):
diff --git a/sympy/functions/elementary/tests/test_miscellaneous.py b/sympy/functions/elementary/tests/test_miscellaneous.py
index 15a61a455eba..a09dde334499 100644
--- a/sympy/functions/elementary/tests/test_miscellaneous.py
+++ b/sympy/functions/elementary/tests/test_miscellaneous.py
@@ -396,7 +396,7 @@ def test_issue_11099():
     assert Max(x, y).evalf(subs=fixed_test_data) == \
         Max(x, y).subs(fixed_test_data).evalf()
     # randomly generate some test data
-    from random import randint
+    from sympy.core.random import randint
     for i in range(20):
         random_test_data = {x: randint(-100, 100), y: randint(-100, 100)}
         assert Min(x, y).evalf(subs=random_test_data) == \
diff --git a/sympy/functions/special/tests/test_bessel.py b/sympy/functions/special/tests/test_bessel.py
index 9904823a55cf..a34247545864 100644
--- a/sympy/functions/special/tests/test_bessel.py
+++ b/sympy/functions/special/tests/test_bessel.py
@@ -18,7 +18,7 @@
 from sympy.series.series import series
 from sympy.functions.special.bessel import (airyai, airybi,
                                             airyaiprime, airybiprime, marcumq)
-from sympy.testing.randtest import (random_complex_number as randcplx,
+from sympy.core.random import (random_complex_number as randcplx,
                                       verify_numerically as tn,
                                       test_derivative_numerically as td,
                                       _randint)
@@ -483,7 +483,7 @@ def test_branching():
     assert besseli(n, polar_lift(x)) == besseli(n, x)
 
     def tn(func, s):
-        from random import uniform
+        from sympy.core.random import uniform
         c = uniform(1, 5)
         expr = func(s, c*exp_polar(I*pi)) - func(s, c*exp_polar(-I*pi))
         eps = 1e-15
diff --git a/sympy/functions/special/tests/test_elliptic_integrals.py b/sympy/functions/special/tests/test_elliptic_integrals.py
index 54b783c0d6a7..52f8c09422cb 100644
--- a/sympy/functions/special/tests/test_elliptic_integrals.py
+++ b/sympy/functions/special/tests/test_elliptic_integrals.py
@@ -10,7 +10,7 @@
 from sympy.series.order import O
 from sympy.functions.special.elliptic_integrals import (elliptic_k as K,
     elliptic_f as F, elliptic_e as E, elliptic_pi as P)
-from sympy.testing.randtest import (test_derivative_numerically as td,
+from sympy.core.random import (test_derivative_numerically as td,
                                       random_complex_number as randcplx,
                                       verify_numerically as tn)
 from sympy.abc import z, m, n
diff --git a/sympy/functions/special/tests/test_error_functions.py b/sympy/functions/special/tests/test_error_functions.py
index 71e84db48a15..3d8200148ba0 100644
--- a/sympy/functions/special/tests/test_error_functions.py
+++ b/sympy/functions/special/tests/test_error_functions.py
@@ -355,7 +355,7 @@ def test_erf2inv():
 
 
 def mytn(expr1, expr2, expr3, x, d=0):
-    from sympy.testing.randtest import verify_numerically, random_complex_number
+    from sympy.core.random import verify_numerically, random_complex_number
     subs = {}
     for a in expr1.free_symbols:
         if a != x:
@@ -365,7 +365,7 @@ def mytn(expr1, expr2, expr3, x, d=0):
 
 
 def mytd(expr1, expr2, x):
-    from sympy.testing.randtest import test_derivative_numerically, \
+    from sympy.core.random import test_derivative_numerically, \
         random_complex_number
     subs = {}
     for a in expr1.free_symbols:
@@ -375,7 +375,7 @@ def mytd(expr1, expr2, x):
 
 
 def tn_branch(func, s=None):
-    from random import uniform
+    from sympy.core.random import uniform
 
     def fn(x):
         if s is None:
@@ -508,7 +508,7 @@ def test__eis():
 
 def tn_arg(func):
     def test(arg, e1, e2):
-        from random import uniform
+        from sympy.core.random import uniform
         v = uniform(1, 5)
         v1 = func(arg*x).subs(x, v).n()
         v2 = func(e1*v + e2*1e-15).n()
@@ -776,7 +776,7 @@ def test_fresnel():
         meijerg(((), (1,)), ((Rational(1, 4),),
         (Rational(3, 4), 0)), -pi**2*z**4/16)/(2*(-z)**Rational(1, 4)*(z**2)**Rational(1, 4))
 
-    from sympy.testing.randtest import verify_numerically
+    from sympy.core.random import verify_numerically
 
     verify_numerically(re(fresnels(z)), fresnels(z).as_real_imag()[0], z)
     verify_numerically(im(fresnels(z)), fresnels(z).as_real_imag()[1], z)
diff --git a/sympy/functions/special/tests/test_gamma_functions.py b/sympy/functions/special/tests/test_gamma_functions.py
index 4856db15bbdd..0a200c534be3 100644
--- a/sympy/functions/special/tests/test_gamma_functions.py
+++ b/sympy/functions/special/tests/test_gamma_functions.py
@@ -18,7 +18,7 @@
 from sympy.core.expr import unchanged
 from sympy.core.function import ArgumentIndexError
 from sympy.testing.pytest import raises
-from sympy.testing.randtest import (test_derivative_numerically as td,
+from sympy.core.random import (test_derivative_numerically as td,
                                       random_complex_number as randcplx,
                                       verify_numerically as tn)
 
@@ -109,7 +109,7 @@ def test_gamma_series():
 
 
 def tn_branch(s, func):
-    from random import uniform
+    from sympy.core.random import uniform
     c = uniform(1, 5)
     expr = func(s, c*exp_polar(I*pi)) - func(s, c*exp_polar(-I*pi))
     eps = 1e-15
diff --git a/sympy/functions/special/tests/test_hyper.py b/sympy/functions/special/tests/test_hyper.py
index e8113319f764..716deac580db 100644
--- a/sympy/functions/special/tests/test_hyper.py
+++ b/sympy/functions/special/tests/test_hyper.py
@@ -12,7 +12,7 @@
 from sympy.abc import x, z, k
 from sympy.series.limits import limit
 from sympy.testing.pytest import raises, slow
-from sympy.testing.randtest import (
+from sympy.core.random import (
     random_complex_number as randcplx,
     verify_numerically as tn,
     test_derivative_numerically as td)
diff --git a/sympy/functions/special/tests/test_zeta_functions.py b/sympy/functions/special/tests/test_zeta_functions.py
index e2d98c073432..90dbeb5855e1 100644
--- a/sympy/functions/special/tests/test_zeta_functions.py
+++ b/sympy/functions/special/tests/test_zeta_functions.py
@@ -11,7 +11,7 @@
 from sympy.core.function import ArgumentIndexError
 from sympy.functions.combinatorial.numbers import bernoulli, factorial
 from sympy.testing.pytest import raises
-from sympy.testing.randtest import (test_derivative_numerically as td,
+from sympy.core.random import (test_derivative_numerically as td,
                       random_complex_number as randcplx, verify_numerically)
 
 x = Symbol('x')
diff --git a/sympy/geometry/tests/test_entity.py b/sympy/geometry/tests/test_entity.py
index e279f86efb43..cecfdb785506 100644
--- a/sympy/geometry/tests/test_entity.py
+++ b/sympy/geometry/tests/test_entity.py
@@ -6,8 +6,6 @@
 from sympy.geometry.entity import scale, GeometryEntity
 from sympy.testing.pytest import raises
 
-from random import random
-
 
 def test_entity():
     x = Symbol('x', real=True)
@@ -85,8 +83,10 @@ def test_reflect_entity_overrides():
     assert c.area == -cr.area
 
     pent = RegularPolygon((1, 2), 1, 5)
-    l = Line(pent.vertices[1],
-        slope=Rational(random() - .5, random() - .5))
+    slope = S.ComplexInfinity
+    while slope is S.ComplexInfinity:
+        slope = Rational(*(x._random()/2).as_real_imag())
+    l = Line(pent.vertices[1], slope=slope)
     rpent = pent.reflect(l)
     assert rpent.center == pent.center.reflect(l)
     rvert = [i.reflect(l) for i in pent.vertices]
diff --git a/sympy/geometry/tests/test_polygon.py b/sympy/geometry/tests/test_polygon.py
index ab9216af869a..4f0507ee44d2 100644
--- a/sympy/geometry/tests/test_polygon.py
+++ b/sympy/geometry/tests/test_polygon.py
@@ -9,7 +9,7 @@
                             Polygon, Ray, RegularPolygon, Segment, Triangle,
                             are_similar, convex_hull, intersection, Line, Ray2D)
 from sympy.testing.pytest import raises, slow, warns
-from sympy.testing.randtest import verify_numerically
+from sympy.core.random import verify_numerically
 from sympy.geometry.polygon import rad, deg
 from sympy.integrals.integrals import integrate
 
diff --git a/sympy/geometry/tests/test_util.py b/sympy/geometry/tests/test_util.py
index 45e1a7bcf230..51e119edb452 100644
--- a/sympy/geometry/tests/test_util.py
+++ b/sympy/geometry/tests/test_util.py
@@ -73,7 +73,7 @@ def test_centroid():
 
 
 def test_farthest_points_closest_points():
-    from random import randint
+    from sympy.core.random import randint
     from sympy.utilities.iterables import subsets
 
     for how in (min, max):
diff --git a/sympy/integrals/rubi/utility_function.py b/sympy/integrals/rubi/utility_function.py
index b64bac06f30d..a28d59f850e8 100644
--- a/sympy/integrals/rubi/utility_function.py
+++ b/sympy/integrals/rubi/utility_function.py
@@ -44,7 +44,7 @@
 from sympy.simplify.simplify import fraction, simplify, cancel, powsimp, nsimplify
 from sympy.utilities.decorator import doctest_depends_on
 from sympy.utilities.iterables import flatten
-from random import randint
+from sympy.core.random import randint
 
 
 class rubi_unevaluated_expr(UnevaluatedExpr):
diff --git a/sympy/integrals/tests/test_integrals.py b/sympy/integrals/tests/test_integrals.py
index b38398aa55d0..5d0634245fec 100644
--- a/sympy/integrals/tests/test_integrals.py
+++ b/sympy/integrals/tests/test_integrals.py
@@ -39,7 +39,7 @@
 from sympy.physics import units
 from sympy.testing.pytest import (raises, slow, skip, ON_TRAVIS,
     warns_deprecated_sympy)
-from sympy.testing.randtest import verify_numerically
+from sympy.core.random import verify_numerically
 
 
 x, y, a, t, x_1, x_2, z, s, b = symbols('x y a t x_1 x_2 z s b')
diff --git a/sympy/integrals/tests/test_meijerint.py b/sympy/integrals/tests/test_meijerint.py
index 7e2961ebb606..1f3f136b8320 100644
--- a/sympy/integrals/tests/test_meijerint.py
+++ b/sympy/integrals/tests/test_meijerint.py
@@ -16,7 +16,7 @@
     meijerint_indefinite, _inflate_g, _create_lookup_table,
     meijerint_definite, meijerint_inversion)
 from sympy.testing.pytest import slow
-from sympy.testing.randtest import (verify_numerically,
+from sympy.core.random import (verify_numerically,
         random_complex_number as randcplx)
 from sympy.abc import x, y, a, b, c, d, s, t, z
 
@@ -338,7 +338,7 @@ def test_inversion_exp_real_nonreal_shift():
 
 @slow
 def test_lookup_table():
-    from random import uniform, randrange
+    from sympy.core.random import uniform, randrange
     from sympy.core.add import Add
     from sympy.integrals.meijerint import z as z_dummy
     table = {}
diff --git a/sympy/ntheory/residue_ntheory.py b/sympy/ntheory/residue_ntheory.py
index 101c08a7e23c..4f629b0ebc95 100644
--- a/sympy/ntheory/residue_ntheory.py
+++ b/sympy/ntheory/residue_ntheory.py
@@ -6,7 +6,7 @@
 from .primetest import isprime
 from .factor_ import factorint, trailing, totient, multiplicity
 from sympy.utilities.misc import as_int
-from random import randint, Random
+from sympy.core.random import _randint, randint
 
 from itertools import cycle, product
 
@@ -1196,13 +1196,11 @@ def _discrete_log_pollard_rho(n, a, b, order=None, retries=10, rseed=None):
 
     if order is None:
         order = n_order(b, n)
-    prng = Random()
-    if rseed is not None:
-        prng.seed(rseed)
+    randint = _randint(rseed)
 
     for i in range(retries):
-        aa = prng.randint(1, order - 1)
-        ba = prng.randint(1, order - 1)
+        aa = randint(1, order - 1)
+        ba = randint(1, order - 1)
         xa = pow(b, aa, n) * pow(a, ba, n) % n
 
         c = xa % 3
diff --git a/sympy/ntheory/tests/test_bbp_pi.py b/sympy/ntheory/tests/test_bbp_pi.py
index 806945a6a550..c18188e3c8a3 100644
--- a/sympy/ntheory/tests/test_bbp_pi.py
+++ b/sympy/ntheory/tests/test_bbp_pi.py
@@ -1,4 +1,4 @@
-from random import randint
+from sympy.core.random import randint
 
 from sympy.ntheory.bbp_pi import pi_hex_digits
 from sympy.testing.pytest import raises
diff --git a/sympy/ntheory/tests/test_egyptian_fraction.py b/sympy/ntheory/tests/test_egyptian_fraction.py
index 10e2ef3e214e..a9a9fac578d9 100644
--- a/sympy/ntheory/tests/test_egyptian_fraction.py
+++ b/sympy/ntheory/tests/test_egyptian_fraction.py
@@ -2,7 +2,7 @@
 from sympy.ntheory.egyptian_fraction import egyptian_fraction
 from sympy.core.add import Add
 from sympy.testing.pytest import raises
-from sympy.testing.randtest import random_complex_number
+from sympy.core.random import random_complex_number
 
 
 def test_egyptian_fraction():
diff --git a/sympy/physics/quantum/circuitutils.py b/sympy/physics/quantum/circuitutils.py
index 10a634daf8b6..95fa3ce15e10 100644
--- a/sympy/physics/quantum/circuitutils.py
+++ b/sympy/physics/quantum/circuitutils.py
@@ -406,7 +406,7 @@ def random_reduce(circuit, gate_ids, seed=None):
         given by the list
 
     """
-    from sympy.testing.randtest import _randrange
+    from sympy.core.random import _randrange
 
     if not gate_ids:
         return circuit
@@ -460,7 +460,7 @@ def random_insert(circuit, choices, seed=None):
     Indices for insertion should be [0, n] if n is the length of the
     circuit.
     """
-    from sympy.testing.randtest import _randrange
+    from sympy.core.random import _randrange
 
     if not choices:
         return circuit
diff --git a/sympy/physics/quantum/identitysearch.py b/sympy/physics/quantum/identitysearch.py
index 0dac602c7451..04b00ef208e4 100644
--- a/sympy/physics/quantum/identitysearch.py
+++ b/sympy/physics/quantum/identitysearch.py
@@ -1,5 +1,5 @@
 from collections import deque
-from random import randint
+from sympy.core.random import randint
 
 from sympy.external import import_module
 from sympy.core.basic import Basic
diff --git a/sympy/physics/quantum/tests/test_matrixutils.py b/sympy/physics/quantum/tests/test_matrixutils.py
index 8726abfd30fe..af581f7aa236 100644
--- a/sympy/physics/quantum/tests/test_matrixutils.py
+++ b/sympy/physics/quantum/tests/test_matrixutils.py
@@ -1,4 +1,4 @@
-from random import randint
+from sympy.core.random import randint
 
 from sympy.core.numbers import Integer
 from sympy.matrices.dense import (Matrix, ones, zeros)
diff --git a/sympy/polys/factortools.py b/sympy/polys/factortools.py
index 0b7da8e2d26b..c322bef488a2 100644
--- a/sympy/polys/factortools.py
+++ b/sympy/polys/factortools.py
@@ -977,7 +977,7 @@ def dmp_zz_wang(f, u, K, mod=None, seed=None):
     .. [2] [Geddes92]_
 
     """
-    from sympy.testing.randtest import _randint
+    from sympy.core.random import _randint
 
     randint = _randint(seed)
 
diff --git a/sympy/polys/galoistools.py b/sympy/polys/galoistools.py
index ea52353d0581..9ffe6b5dbb45 100644
--- a/sympy/polys/galoistools.py
+++ b/sympy/polys/galoistools.py
@@ -1,7 +1,7 @@
 """Dense univariate polynomials with coefficients in Galois fields. """
 
 
-from random import uniform
+from sympy.core.random import uniform
 from math import ceil as _ceil, sqrt as _sqrt
 
 from sympy.core.mul import prod
diff --git a/sympy/polys/tests/test_polyroots.py b/sympy/polys/tests/test_polyroots.py
index fa5900298e01..18a389a41a4e 100644
--- a/sympy/polys/tests/test_polyroots.py
+++ b/sympy/polys/tests/test_polyroots.py
@@ -23,7 +23,7 @@
 from sympy.polys.polyutils import _nsort
 
 from sympy.testing.pytest import raises, slow
-from sympy.testing.randtest import verify_numerically
+from sympy.core.random import verify_numerically
 import mpmath
 from itertools import product
 
diff --git a/sympy/simplify/tests/test_cse.py b/sympy/simplify/tests/test_cse.py
index 9a0cfa323090..efb39bc2d3ea 100644
--- a/sympy/simplify/tests/test_cse.py
+++ b/sympy/simplify/tests/test_cse.py
@@ -454,7 +454,7 @@ def test_issue_11230():
     assert not any(i.is_Mul for a in C for i in a.args)
 
     # random tests for the issue
-    from random import choice
+    from sympy.core.random import choice
     from sympy.core.function import expand_mul
     s = symbols('a:m')
     # 35 Mul tests, none of which should ever fail
diff --git a/sympy/simplify/tests/test_fu.py b/sympy/simplify/tests/test_fu.py
index 94cd0c8cc368..3bf4f74c522c 100644
--- a/sympy/simplify/tests/test_fu.py
+++ b/sympy/simplify/tests/test_fu.py
@@ -12,7 +12,7 @@
     TR111, TR2, TR2i, TR3, TR5, TR6, TR7, TR8, TR9, TRmorrie, _TR56 as T,
     TRpower, hyper_as_trig, fu, process_common_addends, trig_split,
     as_f_sign_1)
-from sympy.testing.randtest import verify_numerically
+from sympy.core.random import verify_numerically
 from sympy.abc import a, b, c, x, y, z
 
 
diff --git a/sympy/simplify/tests/test_hyperexpand.py b/sympy/simplify/tests/test_hyperexpand.py
index f770865e84d3..146dadc183eb 100644
--- a/sympy/simplify/tests/test_hyperexpand.py
+++ b/sympy/simplify/tests/test_hyperexpand.py
@@ -1,4 +1,4 @@
-from random import randrange
+from sympy.core.random import randrange
 
 from sympy.simplify.hyperexpand import (ShiftA, ShiftB, UnShiftA, UnShiftB,
                        MeijerShiftA, MeijerShiftB, MeijerShiftC, MeijerShiftD,
@@ -20,7 +20,7 @@
 from sympy.functions.special.hyper import (hyper, meijerg)
 from sympy.abc import z, a, b, c
 from sympy.testing.pytest import XFAIL, raises, slow, ON_TRAVIS, skip
-from sympy.testing.randtest import verify_numerically as tn
+from sympy.core.random import verify_numerically as tn
 
 from sympy.core.numbers import (Rational, pi)
 from sympy.functions.elementary.exponential import (exp, exp_polar, log)
diff --git a/sympy/simplify/tests/test_trigsimp.py b/sympy/simplify/tests/test_trigsimp.py
index a960b3b9bfdd..3fda27955c95 100644
--- a/sympy/simplify/tests/test_trigsimp.py
+++ b/sympy/simplify/tests/test_trigsimp.py
@@ -379,7 +379,7 @@ def test_issue_15129_trigsimp_methods():
 
 def test_exptrigsimp():
     def valid(a, b):
-        from sympy.testing.randtest import verify_numerically as tn
+        from sympy.core.random import verify_numerically as tn
         if not (tn(a, b) and a == b):
             return False
         return True
diff --git a/sympy/solvers/diophantine/tests/test_diophantine.py b/sympy/solvers/diophantine/tests/test_diophantine.py
index d98971a94435..2f1744e98d56 100644
--- a/sympy/solvers/diophantine/tests/test_diophantine.py
+++ b/sympy/solvers/diophantine/tests/test_diophantine.py
@@ -658,7 +658,7 @@ def test_sum_of_three_squares():
 
 
 def test_sum_of_four_squares():
-    from random import randint
+    from sympy.core.random import randint
 
     # this should never fail
     n = randint(1, 100000000000000)
diff --git a/sympy/solvers/ode/tests/test_riccati.py b/sympy/solvers/ode/tests/test_riccati.py
index 24428d0a52d8..1c9487afed93 100644
--- a/sympy/solvers/ode/tests/test_riccati.py
+++ b/sympy/solvers/ode/tests/test_riccati.py
@@ -1,4 +1,4 @@
-from random import randint
+from sympy.core.random import randint
 from sympy.core.function import Function
 from sympy.core.mul import Mul
 from sympy.core.numbers import (I, Rational, oo)
diff --git a/sympy/solvers/tests/test_solvers.py b/sympy/solvers/tests/test_solvers.py
index d32e2eb05a53..942810fb63c9 100644
--- a/sympy/solvers/tests/test_solvers.py
+++ b/sympy/solvers/tests/test_solvers.py
@@ -37,7 +37,7 @@
 from sympy.polys.rootoftools import CRootOf
 
 from sympy.testing.pytest import slow, XFAIL, SKIP, raises
-from sympy.testing.randtest import verify_numerically as tn
+from sympy.core.random import verify_numerically as tn
 
 from sympy.abc import a, b, c, d, e, k, h, p, x, y, z, t, q, m, R
 
diff --git a/sympy/solvers/tests/test_solveset.py b/sympy/solvers/tests/test_solveset.py
index 5b230469083f..236a56dd7c3b 100644
--- a/sympy/solvers/tests/test_solveset.py
+++ b/sympy/solvers/tests/test_solveset.py
@@ -32,7 +32,7 @@
 from sympy.utilities.iterables import numbered_symbols
 
 from sympy.testing.pytest import (XFAIL, raises, skip, slow, SKIP, _both_exp_pow)
-from sympy.testing.randtest import verify_numerically as tn
+from sympy.core.random import verify_numerically as tn
 from sympy.physics.units import cm
 
 from sympy.solvers import solve
diff --git a/sympy/stats/tests/test_continuous_rv.py b/sympy/stats/tests/test_continuous_rv.py
index a9eb50845ec7..d74525c87629 100644
--- a/sympy/stats/tests/test_continuous_rv.py
+++ b/sympy/stats/tests/test_continuous_rv.py
@@ -45,7 +45,7 @@
 from sympy.stats.compound_rv import CompoundPSpace
 from sympy.stats.symbolic_probability import Probability
 from sympy.testing.pytest import raises, XFAIL, slow, ignore_warnings
-from sympy.testing.randtest import verify_numerically as tn
+from sympy.core.random import verify_numerically as tn
 
 oo = S.Infinity
 
diff --git a/sympy/tensor/tests/test_tensor_operators.py b/sympy/tensor/tests/test_tensor_operators.py
index a3aca8dc7898..996d0af06070 100644
--- a/sympy/tensor/tests/test_tensor_operators.py
+++ b/sympy/tensor/tests/test_tensor_operators.py
@@ -9,7 +9,7 @@
 from sympy.matrices.dense import diag
 from sympy.tensor.array import Array
 
-from random import randint
+from sympy.core.random import randint
 
 
 L = TensorIndexType("L")
diff --git a/sympy/testing/randtest.py b/sympy/testing/randtest.py
index 4d5232acd521..588496e9be05 100644
--- a/sympy/testing/randtest.py
+++ b/sympy/testing/randtest.py
@@ -1,174 +1,26 @@
-""" Helpers for randomized testing """
-
-from random import uniform, Random, randrange, randint
-
-from sympy.core.containers import Tuple
-from sympy.core.function import Derivative
-from sympy.core.numbers import comp, I
-from sympy.core.symbol import Symbol
-from sympy.simplify.simplify import nsimplify
-from sympy.utilities.iterables import is_sequence
-from sympy.utilities.misc import as_int
-
-
-def random_complex_number(a=2, b=-1, c=3, d=1, rational=False, tolerance=None):
-    """
-    Return a random complex number.
-
-    To reduce chance of hitting branch cuts or anything, we guarantee
-    b <= Im z <= d, a <= Re z <= c
-
-    When rational is True, a rational approximation to a random number
-    is obtained within specified tolerance, if any.
-    """
-    A, B = uniform(a, c), uniform(b, d)
-    if not rational:
-        return A + I*B
-    return (nsimplify(A, rational=True, tolerance=tolerance) +
-        I*nsimplify(B, rational=True, tolerance=tolerance))
-
-
-def verify_numerically(f, g, z=None, tol=1.0e-6, a=2, b=-1, c=3, d=1):
-    """
-    Test numerically that f and g agree when evaluated in the argument z.
-
-    If z is None, all symbols will be tested. This routine does not test
-    whether there are Floats present with precision higher than 15 digits
-    so if there are, your results may not be what you expect due to round-
-    off errors.
-
-    Examples
-    ========
-
-    >>> from sympy import sin, cos
-    >>> from sympy.abc import x
-    >>> from sympy.testing.randtest import verify_numerically as tn
-    >>> tn(sin(x)**2 + cos(x)**2, 1, x)
-    True
-    """
-    f, g, z = Tuple(f, g, z)
-    z = [z] if isinstance(z, Symbol) else (f.free_symbols | g.free_symbols)
-    reps = list(zip(z, [random_complex_number(a, b, c, d) for _ in z]))
-    z1 = f.subs(reps).n()
-    z2 = g.subs(reps).n()
-    return comp(z1, z2, tol)
-
-
-def test_derivative_numerically(f, z, tol=1.0e-6, a=2, b=-1, c=3, d=1):
-    """
-    Test numerically that the symbolically computed derivative of f
-    with respect to z is correct.
-
-    This routine does not test whether there are Floats present with
-    precision higher than 15 digits so if there are, your results may
-    not be what you expect due to round-off errors.
-
-    Examples
-    ========
-
-    >>> from sympy import sin
-    >>> from sympy.abc import x
-    >>> from sympy.testing.randtest import test_derivative_numerically as td
-    >>> td(sin(x), x)
-    True
-    """
-    z0 = random_complex_number(a, b, c, d)
-    f1 = f.diff(z).subs(z, z0)
-    f2 = Derivative(f, z).doit_numerically(z0)
-    return comp(f1.n(), f2.n(), tol)
-
-
-def _randrange(seed=None):
-    """Return a randrange generator. ``seed`` can be
-        o None - return randomly seeded generator
-        o int - return a generator seeded with the int
-        o list - the values to be returned will be taken from the list
-          in the order given; the provided list is not modified.
-
-    Examples
-    ========
-
-    >>> from sympy.testing.randtest import _randrange
-    >>> rr = _randrange()
-    >>> rr(1000) # doctest: +SKIP
-    999
-    >>> rr = _randrange(3)
-    >>> rr(1000) # doctest: +SKIP
-    238
-    >>> rr = _randrange([0, 5, 1, 3, 4])
-    >>> rr(3), rr(3)
-    (0, 1)
-    """
-    if seed is None:
-        return randrange
-    elif isinstance(seed, int):
-        return Random(seed).randrange
-    elif is_sequence(seed):
-        seed = list(seed)  # make a copy
-        seed.reverse()
-
-        def give(a, b=None, seq=seed):
-            if b is None:
-                a, b = 0, a
-            a, b = as_int(a), as_int(b)
-            w = b - a
-            if w < 1:
-                raise ValueError('_randrange got empty range')
-            try:
-                x = seq.pop()
-            except IndexError:
-                raise ValueError('_randrange sequence was too short')
-            if a <= x < b:
-                return x
-            else:
-                return give(a, b, seq)
-        return give
-    else:
-        raise ValueError('_randrange got an unexpected seed')
-
-
-def _randint(seed=None):
-    """Return a randint generator. ``seed`` can be
-        o None - return randomly seeded generator
-        o int - return a generator seeded with the int
-        o list - the values to be returned will be taken from the list
-          in the order given; the provided list is not modified.
-
-    Examples
-    ========
-
-    >>> from sympy.testing.randtest import _randint
-    >>> ri = _randint()
-    >>> ri(1, 1000) # doctest: +SKIP
-    999
-    >>> ri = _randint(3)
-    >>> ri(1, 1000) # doctest: +SKIP
-    238
-    >>> ri = _randint([0, 5, 1, 2, 4])
-    >>> ri(1, 3), ri(1, 3)
-    (1, 2)
-    """
-    if seed is None:
-        return randint
-    elif isinstance(seed, int):
-        return Random(seed).randint
-    elif is_sequence(seed):
-        seed = list(seed)  # make a copy
-        seed.reverse()
-
-        def give(a, b, seq=seed):
-            a, b = as_int(a), as_int(b)
-            w = b - a
-            if w < 0:
-                raise ValueError('_randint got empty range')
-            try:
-                x = seq.pop()
-            except IndexError:
-                raise ValueError('_randint sequence was too short')
-            if a <= x <= b:
-                return x
-            else:
-                return give(a, b, seq)
-        return give
-    else:
-        raise ValueError('_randint got an unexpected seed')
+from sympy.core.random import (
+    random_complex_number as A,
+    verify_numerically as B,
+    test_derivative_numerically as C,
+    _randrange as D,
+    _randint as E)
+from sympy.utilities.decorator import deprecated
+
+random_complex_number = deprecated(useinstead="sympy.core.random.random_complex_number",
+    deprecated_since_version="1.10", issue=22433)(A)
+
+
+verify_numerically = deprecated(useinstead="sympy.core.random.verify_numerically",
+    deprecated_since_version="1.10", issue=22433)(B)
+
+
+test_derivative_numerically = deprecated(useinstead="sympy.core.random.test_derivative_numerically",
+    deprecated_since_version="1.10", issue=22433)(C)
+
+
+_randrange = deprecated(useinstead="sympy.core.random._randint",
+    deprecated_since_version="1.10", issue=22433)(D)
+
+
+_randint = deprecated(useinstead="sympy.core.random._randint",
+    deprecated_since_version="1.10", issue=22433)(E)
diff --git a/sympy/utilities/iterables.py b/sympy/utilities/iterables.py
index 111d3f31bd12..cac0d2e1e081 100644
--- a/sympy/utilities/iterables.py
+++ b/sympy/utilities/iterables.py
@@ -2313,7 +2313,6 @@ def random_derangement(t, choice=None, strict=True):
     ========
 
     >>> from sympy.utilities.iterables import random_derangement
-    >>> from random import Random
     >>> t = 'SymPy: a CAS in pure Python'
     >>> d = random_derangement(t)
     >>> all(i != j for i, j in zip(d, t))
@@ -2322,6 +2321,8 @@ def random_derangement(t, choice=None, strict=True):
     A predictable result can be obtained by using a pseudorandom
     generator for the choice:
 
+    >>> from sympy.core.random import random, Random
+    >>> random.seed(1)
     >>> c = Random(1).choice
     >>> d = [''.join(random_derangement(t, c)) for i in range(5)]
     >>> assert len(set(d)) != 1  # we got different values
diff --git a/sympy/utilities/randtest.py b/sympy/utilities/randtest.py
index d9fb1ec4237f..567fd919b9f2 100644
--- a/sympy/utilities/randtest.py
+++ b/sympy/utilities/randtest.py
@@ -2,8 +2,8 @@
 
 SymPyDeprecationWarning(
     feature="Import sympy.utilities.randtest",
-    useinstead="Import from sympy.testing.randtest",
+    useinstead="Import from sympy.core.random",
     issue=18095,
     deprecated_since_version="1.6").warn()
 
-from sympy.testing.randtest import *  # noqa:F401
+from sympy.core.random import *  # noqa:F401


Code review comments:
diff_hunk: @@ -0,0 +1,209 @@
+"""When you need to use random, import it from here so there is only
+one generator working for SymPy. Imports from here should behave the
+same as if they were being imported from Python's random module. But
+only the routines currently used in SymPy are included here. To use
+others import `rng` and access the method directly. For example, to
+capture the current state of the generator use `rng.getstate()`.
+
+There is intentionally no Random to import from here. If you want
+to control the state of the generator, import `seed` and call it
+with or without an argument to set the state.
+
+EXAMPLES
+========
+
+>>> from sympy.core.random import random, seed
+>>> assert random() < 1
+>>> seed(1); a = random()
+>>> b = random()
+>>> seed(1); c = random()
+>>> assert a == c
+>>> assert a != b  # remote possibility this will fail
+"""
+
+import random as _random
+rng = _random.Random()
+
+seed = lambda s=None: rng.seed(s)
+
+choice = rng.choice
+random = rng.random
+randint = rng.randint
+randrange = rng.randrange
+shuffle = rng.shuffle
path: sympy/core/random.py
text: @author:
```suggestion
choice = rng.choice
random = rng.random
randint = rng.randint
randrange = rng.randrange
seed = rng.seed
shuffle = rng.shuffle
```