Repository: scipy/scipy
Pull request: MAINT: special.hankel2: fix edge case
#### Reference issue
Closes gh-4517

#### What does this implement/fix?
gh-4517 reported that the output of `hankel2` was inaccurate at 0, 0. This PR fixes that issue (assuming the definition in terms of $J_n$ and $Y_n$ is appropriate here) and only that issue.

#### Additional information
There are undoubtedly other edge cases in related functions (and through `special`) that are not perfect, but since this one was deemed special enough to be open as its own issue for some time, I thought it might as well be fixed. Rather than wait for all of them to be fixed to close this issue, perhaps we need a meta-issue to find and fix edge cases if there is not one already, or maybe it is OK to wait to fix them as they are reported given that there are other issues.

Problem statement:
MAINT: special.hankel2: `(0, 0)` delivers (nan+nan*j) instead of (1+nan*j)
When evaulating the zeroth-order Hankel function of second kind at zero

```
In[121]: scipy.special.hankel2(0,0)
Out[121]: (nan+nan*j)
```

the real value is nan instead of 1. Bessel and Neumann function deliver the correct results

```
In[126]: scipy.special.jv(0,0)
Out[126]: 1.0

In[127]: scipy.special.yn(0,0)
Out[127]: -inf
```

The combination of both according to the definition of the Hankel function again the wrong result

```
In[125]: scipy.special.jv(0,0) - 1j*scipy.special.yn(0,0)
Out[125]: (nan+inf*j)
```

Hints from issue discussion:
The second result is correct if you follow C99 complex arithmetic (which is what Python and C complex numbers follow): `nan+infj` is one representation for the complex infinity.

The first result is apparently given not correctly by the AMOS codes, which do not calculate non-finite results. Overflow needs to be hardcoded likely.

Checking against two other sources, mpmath produces the same results as
scipy and matlab gives simply nan.

On Wed, Feb 11, 2015 at 10:33 AM, Pauli Virtanen notifications@github.com
wrote:

> The second result is correct if you follow C99 complex arithmetic:
> `nan+infj` is one representation for the complex infinity.
> 
> —
> Reply to this email directly or view it on GitHub
> https://github.com/scipy/scipy/issues/4517#issuecomment-73900758.

Commit message: Apply suggestions from code review

Co-authored-by: Albert Steppi <albert.steppi@gmail.com>

Patch to review:
diff --git a/scipy/special/tests/test_basic.py b/scipy/special/tests/test_basic.py
index efb34b2c5238..32bf2fab2f76 100644
--- a/scipy/special/tests/test_basic.py
+++ b/scipy/special/tests/test_basic.py
@@ -2642,6 +2642,12 @@ def test_hankl2e(self):
         hankrl2e = special.hankel2e(1,.1)
         assert_almost_equal(hank2e,hankrl2e,8)
 
+    def test_hankel2_gh4517(self):
+        # Test edge case reported in https://github.com/scipy/scipy/issues/4517
+        res = special.hankel2(0, 0)
+        assert res.real == 1
+        assert np.isposinf(res.imag)
+
 
 class TestHyper:
     def test_h1vp(self):
diff --git a/scipy/special/xsf/bessel.h b/scipy/special/xsf/bessel.h
index ca88c1ae4661..28dd831b4e16 100644
--- a/scipy/special/xsf/bessel.h
+++ b/scipy/special/xsf/bessel.h
@@ -1172,6 +1172,11 @@ inline std::complex<double> cyl_hankel_2(double v, std::complex<double> z) {
     if (isnan(v) || isnan(z.real()) || isnan(z.imag())) {
         return cy;
     }
+    if (v == 0 && z == 0.0) {
+        cy.real(NAN);
+        cy.imag(INFINITY);
+        return cy;
+    }
     if (v < 0) {
         v = -v;
         sign = -1;


Code review comments:
diff_hunk: @@ -2642,6 +2642,12 @@ def test_hankl2e(self):
         hankrl2e = special.hankel2e(1,.1)
         assert_almost_equal(hank2e,hankrl2e,8)
 
+    def test_hankel2_gh4517(self):
+        # Test edge case reported in https://github.com/scipy/scipy/issues/4517
+        res = special.hankel2(0, 0)
+        assert res.real == 1
path: scipy/special/tests/test_basic.py
text: @author:
```suggestion
        assert np.isnan(res.real)
```