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.

     """