Repository: scipy/scipy
Pull request: MAINT: Begin translation of hyp2f1 for complex numbers into Cython
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
https://docs.scipy.org/doc/numpy/dev/development_workflow.html#writing-the-commit-message

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
<!--Example: Closes gh-WXYZ.-->
Closes gh-7340
#### What does this implement/fix?
<!--Please explain your changes.-->
This re-implements some special cases of the Gaussian hypergeometric function for complex numbers into Cython, falling back to the original Fortran implementation for unhandled cases. It is the second in a series that also includes #14287, which added extensive testing for hyp2f1 for complex numbers. These two PRs contain portions of the work from #14256, a closed PR which is in the process of being broken up into smaller pieces that are less burdensome to review.

The cases that have been translated to Cython are

1. z = 0
2. Any of a, b, or c a non-positive integer.

Although these cases appear simple, there are subtleties here. In some cases, the choice of value to return is up to convention.
gh-7340 pointed out that some of the choices made in the original implementation may be less than ideal. The behavior implemented in this PR follows mpmath. We give some justification of the choices below. For reference the hypergeometric series is as follows.

<p align="center">
  <img src="https://user-images.githubusercontent.com/1953382/126650139-b5cfacee-8ee8-4dc8-bf98-470241fd112a.png" />
</p>

The cases where any of a, b, or c are non-positive integers are considered degenerate. For now let's exclude the case z = 0. If c is a non-positive integer and neither of a or b are non-positive integers, then the result is unambiguously infinite. However if c is a non-positive integer and one of a or b is a non-positive integer with `c <= a <= 0` or `c <= b <= 0` then although the series technically has undefined value, the polynomial given by terminating at `n = |a| - 1` or `n = |b| - 1` respectively is a solution to the corresponding (in this case degenerate) hypergeometric differential equation


<p align="center">
  <img src="https://user-images.githubusercontent.com/1953382/126651534-c7ab42a0-bbf4-4350-96e7-285317428c29.png" />
</p>

and it's therefore reasonable to return the value of this polynomial evaluated at z. 

For the case z = 0. hyp2f1(a, b, c, 0) = 1, for all choices of a, b, c, except perhaps for c a non-positive integer. As above, when c is a non-positive integer, the choice of value to return is up to convention. The behavior in this PR is to follow mpmath in these cases. That is, hyp2f1(a, b, c, 0) = 1, except for in the case where c is 0. If c = 0 and one of a or b is 0, then hyp2f1(a, b, 0, 0) = 1. If c = 0 and neither a or b is 0, then hyp2f1(a, b, 0, 0) is undefined and the function returns nan.

I'm not sure what the justification is for returning 1 in all cases where c is a negative integer and z = 0, but I think it's reasonable to follow mpmath here as well.

#### Additional information
<!--Any additional information you think is important.-->


Problem statement:
SciPy Hypergeometric function hyp2f1 producing infinities
<<Please describe the issue in detail here, and for bug reports fill in the fields below.>>
The SciPy hypergeometric function special.hyp2f1 seems to be producing infinities for instances where third argument is negative, and fourth argument is complex. The corresponding implementation in MATLAB gives reasonable finite outputs. See for instance,

### Reproducing code example:
```
<<A short code example that reproduces the problem. It should be self-contained, i.e., possible to run as-is via 'python myproblem.py'>>
```
In [2]: special.hyp2f1(-2,-3,-19,3.4+0.4j)
Out[2]: (inf+0j)

### Error message:
```
<<Full error message, if any (starting from line Traceback: ...)>>
```

### Scipy/Numpy/Python version information:
```
scipy version: 0.19.0

<<Output from 'import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)'>>
```

Hints from issue discussion:
This issue is caused by code in both `specfun_wrappers.c` and `specfun.f` that immediately return `inf` whenever c is a negative integer *without* checking if a or b is a negative integer smaller in magnitude than c (in which case the hypergeometric function is a well defined polynomial).

Since I'm working on hypergeometric function at the moment, I can do a quick PR that fixes this issue. 
A related bug is that hyp2f1 will attempt to calculate a value when c=0 and will return nan+nan*1.j instead of inf...

Commit message: DOC: Improve variable name and clarify comments

Patch to review:
diff --git a/scipy/special/_hyp2f1.pxd b/scipy/special/_hyp2f1.pxd
new file mode 100644
index 000000000000..34ca03a6e6cf
--- /dev/null
+++ b/scipy/special/_hyp2f1.pxd
@@ -0,0 +1,195 @@
+"""
+Implementation of Gauss's hypergeometric function for complex values. This is
+an effort to incrementally translate the Fortran implementation from specfun.f
+into Cython so that it's easier to maintain and it's easier to correct defects
+in the original implementation.  Computation of Gauss's hypergeometric function
+involves handling a patchwork of special cases. The goal is to translate the
+cases into Cython little by little, falling back to the Fortran implementation
+for the cases that are yet to be handled. By default the Fortran original is
+followed as closely as possible except for situations when an improvement is
+obvious. Attempts are made to document the why behind certain decisions made by
+the original implementation, with references to the NIST Digital Library of
+Mathematical Functions [1] added where they are appropriate. The review paper
+by Pearson et al [2] is an excellent resource for best practices for numerical
+computation of hypergeometric functions, and the intent is to follow this paper
+when making improvements to and correcting defects in the original
+implementation.
+
+Author: Albert Steppi
+
+Distributed under the same license as Scipy.
+
+References
+----------
+[1] NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/,
+    Release 1.1.1 of 2021-03-15. F. W. J. Olver, A. B. Olde Daalhuis,
+    D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller,
+    B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
+[2] Pearson, J.W., Olver, S. & Porter, M.A.
+    "Numerical methods for the computation of the confluent and Gauss
+    hypergeometric functions."
+    Numer Algor 74, 821-866 (2017). https://doi.org/10.1007/s11075-016-0173-0
+[3] Raimundas Vidunas, "Degenerate Gauss Hypergeometric Functions",
+    Kyushu Journal of Mathematics, 2007, Volume 61, Issue 1, Pages 109-135,
+"""
+
+cimport cython
+from numpy cimport npy_cdouble
+from libc.math cimport fabs, trunc
+
+from . cimport sf_error
+from ._complexstuff cimport (
+    double_complex_from_npy_cdouble,
+    npy_cdouble_from_double_complex,
+    zabs,
+    zpack,
+)
+
+
+cdef extern from "numpy/npy_math.h":
+    double NPY_NAN
+    double NPY_INFINITY
+
+
+cdef extern from 'specfun_wrappers.h':
+    npy_cdouble chyp2f1_wrap(
+        double,
+        double,
+        double,
+        npy_cdouble,
+    ) nogil
+
+
+# Small value from the Fortran original.
+DEF EPS = 1e-15
+
+
+@cython.cdivision(True)
+cdef inline double complex hyp2f1_complex(
+        double a, double b, double c, double complex z
+) nogil:
+    cdef:
+        int num_terms
+        double modulus_z
+        bint a_neg_int, b_neg_int, c_non_pos_int
+    modulus_z = zabs(z)
+    a_neg_int = a == trunc(a) and a < 0
+    b_neg_int = b == trunc(b) and b < 0
+    c_non_pos_int = c == trunc(c) and c <= 0
+    # Special Cases
+    # -------------------------------------------------------------------------
+    # Takes constant value 1 when a = 0 or b = 0, even if c is a non-positive
+    # integer. This follows mpmath.
+    if a == 0 or b == 0:
+        return 1.0 + 0.0j
+    # Equals 1 when z is 0, unless c is 0.
+    if modulus_z == 0:
+        if c != 0:
+            return 1.0 + 0.0j
+        else:
+            # Returning real part NAN and imaginary part 0 follows mpmath.
+            return zpack(NPY_NAN, 0)
+    # Diverges when c is a non-positive integer unless a is an integer with
+    # c <= a <= 0 or b is an integer with c <= b <= 0, (or z equals 0 with
+    # c != 0) Cases z = 0, a = 0, or b = 0 have already been handled. We follow
+    # mpmath in handling the degenerate cases where any of a, b, c are
+    # non-positive integers. See [3] for a treatment of degenerate cases.
+    if c_non_pos_int and not (
+            a_neg_int and c <= a < 0 or b_neg_int and c <= b < 0
+    ):
+        return NPY_INFINITY + 0.0j
+    # Diverges as real(z) -> 1 when c < a + b.
+    if fabs(1 - z.real) < EPS and z.imag == 0 and c - a - b < 0:
+        return NPY_INFINITY + 0.0j
+    # Gauss's Summation Theorem for z = 1; c - a - b > 0 (DLMF 15.4.20).
+    # Fallback to Fortran original. To be translated to Cython later.
+    if z == 1.0 and c - a - b > 0:
+        return double_complex_from_npy_cdouble(
+            chyp2f1_wrap(a, b, c, npy_cdouble_from_double_complex(z))
+        )
+    # Kummer's Theorem for z = -1; c = 1 + a - b (DLMF 15.4.26).
+    # Fall back to Fortran original. To be translated to Cython later.
+    if zabs(z + 1) < EPS and fabs(1 + a - b - c) < EPS:
+        return double_complex_from_npy_cdouble(
+            chyp2f1_wrap(a, b, c, npy_cdouble_from_double_complex(z))
+        )
+    # Reduces to a polynomial when a or b is a negative integer.
+    # If a and b are both negative integers, we take care to terminate
+    # the series at a or b of smaller magnitude. This is to ensure proper
+    # handling of situations like a < c < b <= 0, a, b, c all non-positive
+    # integers, where terminating at a would lead to a term of the form 0 / 0.
+    if a_neg_int and b_neg_int:
+        max_degree = <int> fabs(a) - 1 if a > b else <int> fabs(b) - 1
+        return hyp2f1_series(
+            a, b, c, z, max_degree, False, 0
+        )
+    elif a_neg_int:
+        return hyp2f1_series(
+            a, b, c, z, <int> fabs(a) - 1, False, 0
+        )
+    elif b_neg_int:
+        return hyp2f1_series(
+            a, b, c, z, <int> fabs(b) - 1, False, 0
+        )
+    # Fall through to original Fortran implementation.
+    # -------------------------------------------------------------------------
+    return double_complex_from_npy_cdouble(
+        chyp2f1_wrap(a, b, c, npy_cdouble_from_double_complex(z))
+    )
+
+
+@cython.cdivision(True)
+cdef inline double complex hyp2f1_series(
+        double a,
+        double b,
+        double c,
+        double complex z,
+        int max_degree,
+        bint early_stop,
+        double rtol,
+) nogil:
+    """Return Truncated Maclaurin series for hyp2f1.
+
+    Series is convergent for |z| < 1 but is only practical for numerical
+    computation when |z| < 0.9.
+
+    Parameters
+    ----------
+    a : double
+    b : double
+    c : double
+    z : double complex
+    max_degree : int
+        Maximum degree of terms before truncating.
+    early_stop : bint
+    rtol : double
+        If early_stop is True, truncate when
+            |current_sum - previous_sum| <= rtol * |current_sum|
+        If max_degree is reached before the stopping condition is satisfied,
+        return nan. If early_stop is False, sum all terms up until and
+        including that with degree max_degree.
+
+    Returns
+    -------
+    double complex
+    """
+    cdef:
+        int k
+        double complex term = 1 + 0j
+        double complex previous = 0 + 0j
+        double complex result = 1 + 0j
+    for k in range(max_degree + 1):
+        # Follows the Fortran original exactly. Using *= degrades precision.
+        # Todo: Check generated assembly to see what's going on.
+        term = term * (a + k) * (b + k) / ((k + 1) * (c + k)) * z
+        previous = result
+        result += term
+        if early_stop and zabs(result - previous) < zabs(result) * rtol:
+            break
+    else:
+        # Return nan if max_degree is exceeded without convergence and
+        # early_stop has been set to True.
+        if early_stop:
+            sf_error.error("hyp2f1", sf_error.NO_RESULT, NULL)
+            result = zpack(NPY_NAN, NPY_NAN)
+    return result
diff --git a/scipy/special/functions.json b/scipy/special/functions.json
index 9a66c30a58a6..d9f992bb7afe 100644
--- a/scipy/special/functions.json
+++ b/scipy/special/functions.json
@@ -696,8 +696,8 @@
         "cephes.h": {
             "hyp2f1": "dddd->d"
         },
-        "specfun_wrappers.h": {
-            "chyp2f1_wrap": "dddD->D"
+        "_hyp2f1.pxd": {
+            "hyp2f1_complex": "dddD->D"
         }
     },
     "hyperu": {
diff --git a/scipy/special/tests/test_hyp2f1.py b/scipy/special/tests/test_hyp2f1.py
index c197eca811a2..121e923ffb0d 100644
--- a/scipy/special/tests/test_hyp2f1.py
+++ b/scipy/special/tests/test_hyp2f1.py
@@ -93,7 +93,6 @@ class TestHyp2f1:
                     expected=1 + 0j,
                     rtol=0
                 ),
-                marks=pytest.mark.xfail(reason="gh-7340")
             ),
             pytest.param(
                 Hyp2f1TestCase(
@@ -104,7 +103,6 @@ class TestHyp2f1:
                     expected=1 + 0j,
                     rtol=0
                 ),
-                marks=pytest.mark.xfail(reason="gh-7340"),
             ),
             pytest.param(
                 Hyp2f1TestCase(
@@ -125,7 +123,6 @@ class TestHyp2f1:
                     expected=np.inf + 0j,
                     rtol=0,
                 ),
-                marks=pytest.mark.xfail(reason="gh-7340"),
             ),
             pytest.param(
                 Hyp2f1TestCase(
@@ -136,7 +133,6 @@ class TestHyp2f1:
                     expected=np.nan + 0j,
                     rtol=0,
                 ),
-                marks=pytest.mark.xfail(reason="gh-7340"),
             ),
             pytest.param(
                 Hyp2f1TestCase(
@@ -147,7 +143,6 @@ class TestHyp2f1:
                     expected=(1.0495404166666666+0.05708208333333334j),
                     rtol=1e-15,
                 ),
-                marks=pytest.mark.xfail(reason="gh-7340"),
             ),
             pytest.param(
                 Hyp2f1TestCase(
@@ -158,7 +153,6 @@ class TestHyp2f1:
                     expected=(1.092966013125+0.13455014673750001j),
                     rtol=1e-15,
                 ),
-                marks=pytest.mark.xfail(reason="gh-7340"),
             ),
             pytest.param(
                 Hyp2f1TestCase(
@@ -169,7 +163,6 @@ class TestHyp2f1:
                     expected=(-0.07712512000000005+0.12752814080000005j),
                     rtol=1e-13,
                 ),
-                marks=pytest.mark.xfail(reason="gh-7340"),
             ),
             pytest.param(
                 Hyp2f1TestCase(
@@ -180,7 +173,6 @@ class TestHyp2f1:
                     expected=(1.6400000000000001+0.6400000000000001j),
                     rtol=1e-13,
                 ),
-                marks=pytest.mark.xfail(reason="gh-7340"),
             ),
         ]
     )
@@ -338,7 +330,7 @@ def test_special_case_z_near_minus_1(self, hyp2f1_test_case):
                     c=4.0013768449590685,
                     z=(0.9473684210526314-0.10526315789473695j),
                     expected=(-0.0003054674127221263-9.261359291755414e-05j),
-                    rtol=9e-11,
+                    rtol=1e-10,
                 ),
             ),
             pytest.param(
@@ -348,7 +340,7 @@ def test_special_case_z_near_minus_1(self, hyp2f1_test_case):
                     c=-3.9316537064827854,
                     z=(1.1578947368421053-0.3157894736842106j),
                     expected=(-0.0020809502580892937-0.0041877333232365095j),
-                    rtol=2e-13,
+                    rtol=5e-12,
                 ),
             ),
             pytest.param(


Code review comments:
diff_hunk: @@ -0,0 +1,195 @@
+"""
+Implementation of Gauss's hypergeometric function for complex values. This is
+an effort to incrementally translate the Fortran implementation from specfun.f
+into Cython so that it's easier to maintain and it's easier to correct defects
+in the original implementation.  Computation of Gauss's hypergeometric function
+involves handling a patchwork of special cases. The goal is to translate the
+cases into Cython little by little, falling back to the Fortran implementation
+for the cases that are yet to be handled. By default the Fortran original is
+followed as closely as possible except for situations when an improvement is
+obvious. Attempts are made to document the why behind certain decisions made by
+the original implementation, with references to the NIST Digital Library of
+Mathematical Functions [1] added where they are appropriate. The review paper
+by Pearson et al [2] is an excellent resource for best practices for numerical
+computation of hypergeometric functions, and the intent is to follow this paper
+when making improvements to and correcting defects in the original
+implementation.
+
+Author: Albert Steppi
+
+Distributed under the same license as Scipy.
+
+References
+----------
+[1] NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/,
+    Release 1.1.1 of 2021-03-15. F. W. J. Olver, A. B. Olde Daalhuis,
+    D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller,
+    B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
+[2] Pearson, J.W., Olver, S. & Porter, M.A.
+    "Numerical methods for the computation of the confluent and Gauss
+    hypergeometric functions."
+    Numer Algor 74, 821-866 (2017). https://doi.org/10.1007/s11075-016-0173-0
+[3] Raimundas Vidunas, "Degenerate Gauss Hypergeometric Functions",
+    Kyushu Journal of Mathematics, 2007, Volume 61, Issue 1, Pages 109-135,
+"""
+
+cimport cython
+from numpy cimport npy_cdouble
+from libc.math cimport fabs, trunc
+
+from . cimport sf_error
+from ._complexstuff cimport (
+    double_complex_from_npy_cdouble,
+    npy_cdouble_from_double_complex,
+    zabs,
+    zpack,
+)
+
+
+cdef extern from "numpy/npy_math.h":
+    double NPY_NAN
+    double NPY_INFINITY
+
+
+cdef extern from 'specfun_wrappers.h':
+    npy_cdouble chyp2f1_wrap(
+        double,
+        double,
+        double,
+        npy_cdouble,
+    ) nogil
+
+
+# Small value from the Fortran original.
+DEF EPS = 1e-15
+
+
+@cython.cdivision(True)
+cdef inline double complex hyp2f1_complex(
+        double a, double b, double c, double complex z
+) nogil:
+    cdef:
+        int num_terms
+        double modulus_z
+        bint a_neg_int, b_neg_int, c_non_pos_int
+    modulus_z = zabs(z)
+    a_neg_int = a == trunc(a) and a < 0
+    b_neg_int = b == trunc(b) and b < 0
+    c_non_pos_int = c == trunc(c) and c <= 0
+    # Special Cases
+    # -------------------------------------------------------------------------
+    # Takes constant value 1 when a = 0 or b = 0, even if c is a non-positive
+    # integer. This follows mpmath.
+    if a == 0 or b == 0:
+        return 1.0 + 0.0j
+    # Equals 1 when z is 0, unless c is 0.
+    if modulus_z == 0:
+        if c != 0:
+            return 1.0 + 0.0j
+        else:
+            # Returning real part NAN and imaginary part 0 follows mpmath.
+            return zpack(NPY_NAN, 0)
+    # Diverges when c is a non-positive integer unless a is an integer with
+    # c <= a <= 0 or b is an integer with c <= b <= 0, (or z equals 0 with
+    # c != 0) Cases z = 0, a = 0, or b = 0 have already been handled. We follow
+    # mpmath in handling the degenerate cases where any of a, b, c are
+    # non-positive integers. See [3] for a treatment of degenerate cases.
+    if c_non_pos_int and not (
+            a_neg_int and c <= a < 0 or b_neg_int and c <= b < 0
+    ):
+        return NPY_INFINITY + 0.0j
+    # Diverges as real(z) -> 1 when c < a + b.
+    if fabs(1 - z.real) < EPS and z.imag == 0 and c - a - b < 0:
+        return NPY_INFINITY + 0.0j
+    # Gauss's Summation Theorem for z = 1; c - a - b > 0 (DLMF 15.4.20).
+    # Fallback to Fortran original. To be translated to Cython later.
+    if z == 1.0 and c - a - b > 0:
+        return double_complex_from_npy_cdouble(
+            chyp2f1_wrap(a, b, c, npy_cdouble_from_double_complex(z))
+        )
+    # Kummer's Theorem for z = -1; c = 1 + a - b (DLMF 15.4.26).
+    # Fall back to Fortran original. To be translated to Cython later.
+    if zabs(z + 1) < EPS and fabs(1 + a - b - c) < EPS:
+        return double_complex_from_npy_cdouble(
+            chyp2f1_wrap(a, b, c, npy_cdouble_from_double_complex(z))
+        )
+    # Reduces to a polynomial when a or b is a negative integer.
+    # If a and b are both negative integers, we take care to terminate
+    # the series at a or b of smaller magnitude. This is to ensure proper
+    # handling of situations like a < c < b <= 0, a, b, c all non-positive
+    # integers, where terminating at a would lead to a term of the form 0 / 0.
+    if a_neg_int and b_neg_int:
+        max_degree = <int> fabs(a) - 1 if a > b else <int> fabs(b) - 1
path: scipy/special/_hyp2f1.pxd
text: @user1:
Casting doubles to ints is problematic for large enough doubles. You usually want to keep it as a double instead and rely on the fact that integer arithmetic on integer doubles will be exact.

@author:
Agreed. This was a result of translating the original Fortran code very strictly. I'll have to think about this though. If `|a|` is such a large double that the cast is problematic, then it's also problematic to even attempt to calculate the result in this way. I think there should be a maximum value of `|a|` or `|b|` where we're allowed to handle the degenerate case `a` a negative int or `b` a negative int by truncating the series at `|a|` or `|b|`.

diff_hunk: @@ -0,0 +1,195 @@
+"""
+Implementation of Gauss's hypergeometric function for complex values. This is
+an effort to incrementally translate the Fortran implementation from specfun.f
+into Cython so that it's easier to maintain and it's easier to correct defects
+in the original implementation.  Computation of Gauss's hypergeometric function
+involves handling a patchwork of special cases. The goal is to translate the
+cases into Cython little by little, falling back to the Fortran implementation
+for the cases that are yet to be handled. By default the Fortran original is
+followed as closely as possible except for situations when an improvement is
+obvious. Attempts are made to document the why behind certain decisions made by
+the original implementation, with references to the NIST Digital Library of
+Mathematical Functions [1] added where they are appropriate. The review paper
+by Pearson et al [2] is an excellent resource for best practices for numerical
+computation of hypergeometric functions, and the intent is to follow this paper
+when making improvements to and correcting defects in the original
+implementation.
+
+Author: Albert Steppi
+
+Distributed under the same license as Scipy.
+
+References
+----------
+[1] NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/,
+    Release 1.1.1 of 2021-03-15. F. W. J. Olver, A. B. Olde Daalhuis,
+    D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller,
+    B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
+[2] Pearson, J.W., Olver, S. & Porter, M.A.
+    "Numerical methods for the computation of the confluent and Gauss
+    hypergeometric functions."
+    Numer Algor 74, 821-866 (2017). https://doi.org/10.1007/s11075-016-0173-0
+[3] Raimundas Vidunas, "Degenerate Gauss Hypergeometric Functions",
+    Kyushu Journal of Mathematics, 2007, Volume 61, Issue 1, Pages 109-135,
+"""
+
+cimport cython
+from numpy cimport npy_cdouble
+from libc.math cimport fabs, trunc
+
+from . cimport sf_error
+from ._complexstuff cimport (
+    double_complex_from_npy_cdouble,
+    npy_cdouble_from_double_complex,
+    zabs,
+    zpack,
+)
+
+
+cdef extern from "numpy/npy_math.h":
+    double NPY_NAN
+    double NPY_INFINITY
+
+
+cdef extern from 'specfun_wrappers.h':
+    npy_cdouble chyp2f1_wrap(
+        double,
+        double,
+        double,
+        npy_cdouble,
+    ) nogil
+
+
+# Small value from the Fortran original.
+DEF EPS = 1e-15
+
+
+@cython.cdivision(True)
+cdef inline double complex hyp2f1_complex(
+        double a, double b, double c, double complex z
+) nogil:
+    cdef:
+        int num_terms
+        double modulus_z
+        bint a_neg_int, b_neg_int, c_non_pos_int
+    modulus_z = zabs(z)
+    a_neg_int = a == trunc(a) and a < 0
+    b_neg_int = b == trunc(b) and b < 0
+    c_non_pos_int = c == trunc(c) and c <= 0
+    # Special Cases
+    # -------------------------------------------------------------------------
+    # Takes constant value 1 when a = 0 or b = 0, even if c is a non-positive
+    # integer. This follows mpmath.
+    if a == 0 or b == 0:
+        return 1.0 + 0.0j
+    # Equals 1 when z is 0, unless c is 0.
+    if modulus_z == 0:
+        if c != 0:
+            return 1.0 + 0.0j
+        else:
+            # Returning real part NAN and imaginary part 0 follows mpmath.
+            return zpack(NPY_NAN, 0)
+    # Diverges when c is a non-positive integer unless a is an integer with
+    # c <= a <= 0 or b is an integer with c <= b <= 0, (or z equals 0 with
+    # c != 0) Cases z = 0, a = 0, or b = 0 have already been handled. We follow
+    # mpmath in handling the degenerate cases where any of a, b, c are
+    # non-positive integers. See [3] for a treatment of degenerate cases.
+    if c_non_pos_int and not (
+            a_neg_int and c <= a < 0 or b_neg_int and c <= b < 0
+    ):
+        return NPY_INFINITY + 0.0j
+    # Diverges as real(z) -> 1 when c < a + b.
+    if fabs(1 - z.real) < EPS and z.imag == 0 and c - a - b < 0:
+        return NPY_INFINITY + 0.0j
+    # Gauss's Summation Theorem for z = 1; c - a - b > 0 (DLMF 15.4.20).
+    # Fallback to Fortran original. To be translated to Cython later.
+    if z == 1.0 and c - a - b > 0:
+        return double_complex_from_npy_cdouble(
+            chyp2f1_wrap(a, b, c, npy_cdouble_from_double_complex(z))
+        )
+    # Kummer's Theorem for z = -1; c = 1 + a - b (DLMF 15.4.26).
+    # Fall back to Fortran original. To be translated to Cython later.
+    if zabs(z + 1) < EPS and fabs(1 + a - b - c) < EPS:
+        return double_complex_from_npy_cdouble(
+            chyp2f1_wrap(a, b, c, npy_cdouble_from_double_complex(z))
+        )
+    # Reduces to a polynomial when a or b is a negative integer.
+    # If a and b are both negative integers, we take care to terminate
+    # the series at a or b of smaller magnitude. This is to ensure proper
+    # handling of situations like a < c < b <= 0, a, b, c all non-positive
+    # integers, where terminating at a would lead to a term of the form 0 / 0.
+    if a_neg_int and b_neg_int:
path: scipy/special/_hyp2f1.pxd
text: @user1:
In these code paths it looks like `max_degree` is unbounded-that's risky. We typically want to bail out at some max number of terms (and return `NaN`) to not freeze someones computer; see what we do for `hyp1f1` for example:

https://github.com/scipy/scipy/blob/master/scipy/special/_hypergeometric.pxd#L97

@author:
I agree.  I actually had things more properly at first but tortured the code to make the translation from fortran more exact. I'll fix this.