Repository: mpmath/mpmath

Issue:
Hankel functions incorrect for large imaginary parts and standard precision
For example:
```
>>> hankel1(1,1+30j)
mpc(real='0.0', imag='0.0')
>>> hankel2(1,1-30j)
mpc(real='0.0', imag='0.0')
```
whereas Mathematica gives
```
In[1]:= HankelH1[1, 1. + 30 I]
Out[1]= -7.25495*10^-15 - 1.17346*10^-14 I
In[2]:= HankelH2[1, 1. - 30 I]
Out[2]= -7.25495*10^-15 + 1.17346*10^-14 I
```
Correct values can be found by increasing the precision, but this should be automatic as clearly Mathematica is doing.

The problem arises because the standard definitions of the functions involve two terms with many digits in common, and there need to be sufficient digits available when they are combined.

Discussion:
@fredrik-johansson: Yeah, a simple fix is to use ``sum_accurately`` in the implementations of both Hankel functions. A better fix is to implement the appropriate asymptotic expansions of these functions.

Solution patch:
diff --git a/mpmath/functions/bessel.py b/mpmath/functions/bessel.py
index 836bed6..7c69afa 100644
--- a/mpmath/functions/bessel.py
+++ b/mpmath/functions/bessel.py
@@ -180,11 +180,17 @@ def h(n):
 
 @defun_wrapped
 def hankel1(ctx,n,x,**kwargs):
-    return ctx.besselj(n,x,**kwargs) + ctx.j*ctx.bessely(n,x,**kwargs)
+    def terms():
+        return [ctx.besselj(n,x,**kwargs),
+                ctx.j*ctx.bessely(n,x,**kwargs)]
+    return ctx.sum_accurately(terms)
 
 @defun_wrapped
 def hankel2(ctx,n,x,**kwargs):
-    return ctx.besselj(n,x,**kwargs) - ctx.j*ctx.bessely(n,x,**kwargs)
+    def terms():
+        return [ctx.besselj(n,x,**kwargs),
+                -ctx.j*ctx.bessely(n,x,**kwargs)]
+    return ctx.sum_accurately(terms)
 
 @defun_wrapped
 def whitm(ctx,k,m,z,**kwargs):


Test patch:
diff --git a/mpmath/tests/test_functions2.py b/mpmath/tests/test_functions2.py
index 188327f..610f7ca 100644
--- a/mpmath/tests/test_functions2.py
+++ b/mpmath/tests/test_functions2.py
@@ -2433,3 +2433,7 @@ def test_issue_634():
 
 def test_issue_908():
     assert mp.besselj(-10+0j, 0+0j) == 0
+
+def test_issue_637():
+    assert hankel1(1, 1 + 30j).ae(-7.25495e-15 - 1.17346e-14j)
+    assert hankel2(1, 1 - 30j).ae(-7.25495e-15 + 1.17346e-14j)