scipy/special/_precompute/gammainc_data.py      | 13 +++++++++--
 scipy/special/_precompute/hyp2f1_data.py        | 30 ++++++++++++++++++++++++-
 scipy/special/_precompute/wright_bessel_data.py | 12 +++++++++-
 3 files changed, 51 insertions(+), 4 deletions(-)

diff --git a/scipy/special/_precompute/gammainc_data.py b/scipy/special/_precompute/gammainc_data.py
index b3f23cf4cb..2e8c18ebef 100644
--- a/scipy/special/_precompute/gammainc_data.py
+++ b/scipy/special/_precompute/gammainc_data.py
@@ -22,14 +22,23 @@ from time import time
 import numpy as np
 from numpy import pi
 
-from scipy.special._mptestutils import mpf2float
-
 try:
     import mpmath as mp
 except ImportError:
     pass
 
 
+def mpf2float(x):
+    """
+    Convert an mpf to the nearest floating point number. Just using
+    float directly doesn't work because of results like this:
+
+    with mp.workdps(50):
+        float(mpf("0.99999999999999999")) = 0.9999999999999999
+
+    """
+    return float(mpmath.nstr(x, 17, min_fixed=0, max_fixed=0))
+
 def gammainc(a, x, dps=50, maxterms=10**8):
     """Compute gammainc exactly like mpmath does but allow for more
     summands in hypercomb. See
diff --git a/scipy/special/_precompute/hyp2f1_data.py b/scipy/special/_precompute/hyp2f1_data.py
index cfbe2c7f3c..4462d2f37e 100644
--- a/scipy/special/_precompute/hyp2f1_data.py
+++ b/scipy/special/_precompute/hyp2f1_data.py
@@ -81,8 +81,36 @@ from multiprocessing import Pool
 
 
 from scipy.special import hyp2f1
-from scipy.special.tests.test_hyp2f1 import mp_hyp2f1
 
+try:
+    import mpmath
+except ImportError:
+    mpmath = MissingModule("mpmath")
+
+
+def mp_hyp2f1(a, b, c, z):
+    """Return mpmath hyp2f1 calculated on same branch as scipy hyp2f1.
+
+    For most values of a,b,c mpmath returns the x - 0j branch of hyp2f1 on the
+    branch cut x=(1,inf) whereas scipy's hyp2f1 calculates the x + 0j branch.
+    Thus, to generate the right comparison values on the branch cut, we
+    evaluate mpmath.hyp2f1 at x + 1e-15*j.
+
+    The exception to this occurs when c-a=-m in which case both mpmath and
+    scipy calculate the x + 0j branch on the branch cut. When this happens
+    mpmath.hyp2f1 will be evaluated at the original z point.
+    """
+    on_branch_cut = z.real > 1.0 and abs(z.imag) < 1.0e-15
+    cond1 = abs(c - a - round(c - a)) < 1.0e-15 and round(c - a) <= 0
+    cond2 = abs(c - b - round(c - b)) < 1.0e-15 and round(c - b) <= 0
+    # Make sure imaginary part is *exactly* zero
+    if on_branch_cut:
+        z = z.real + 0.0j
+    if on_branch_cut and not (cond1 or cond2):
+        z_mpmath = z.real + 1.0e-15j
+    else:
+        z_mpmath = z
+    return complex(mpmath.hyp2f1(a, b, c, z_mpmath))
 
 def get_region(z):
     """Assign numbers for regions where hyp2f1 must be handled differently."""
diff --git a/scipy/special/_precompute/wright_bessel_data.py b/scipy/special/_precompute/wright_bessel_data.py
index 434874a617..2e6a9bc98f 100644
--- a/scipy/special/_precompute/wright_bessel_data.py
+++ b/scipy/special/_precompute/wright_bessel_data.py
@@ -9,7 +9,6 @@ import os
 from time import time
 
 import numpy as np
-from scipy.special._mptestutils import mpf2float
 
 try:
     import mpmath as mp
@@ -27,6 +26,17 @@ def rgamma_cached(x, dps):
         return mp.rgamma(x)
 
 
+def mpf2float(x):
+    """
+    Convert an mpf to the nearest floating point number. Just using
+    float directly doesn't work because of results like this:
+
+    with mp.workdps(50):
+        float(mpf("0.99999999999999999")) = 0.9999999999999999
+
+    """
+    return float(mpmath.nstr(x, 17, min_fixed=0, max_fixed=0))
+
 def mp_wright_bessel(a, b, x, dps=50, maxterms=2000):
     """Compute Wright's generalized Bessel function as Series with mpmath.
     """