Added faster method for besselh with real input (#200)

* Added faster method for besselh with real input
* Added tests for issue #29

Author: eprovst

src/bessel.jl

@@ -314,6 +314,18 @@ function besselh(nu::Float64, k::Integer, z::Complex{Float64})
     return _besselh(nu,Int32(k),z,Int32(1))
 end
 
+function besselh(nu::Float64, k::Integer, x::Float64)
+    # Given that x is real, Jnu(x) and Ynu(x) are also real.
+    if k == 1
+        return complex(besselj(nu, x), bessely(nu, x))
+    elseif k == 2
+        return complex(besselj(nu, x), -bessely(nu, x))
+    else
+        # We emulate ZBESH's behaviour
+        throw(AmosException(1))
+    end
+end
+
 """
     besselhx(nu, [k=1,] z)
 
@@ -583,9 +595,8 @@ end
 for bfn in (:besselh, :besselhx)
     @eval begin
         $bfn(nu, z) = $bfn(nu, 1, z)
-        $bfn(nu::Real, k::Integer, x::Real) = $bfn(nu, k, float(x))
-        $bfn(nu::Real, k::Integer, x::AbstractFloat) = $bfn(float(nu), k, complex(x))
-
+        $bfn(nu::Real, k::Integer, x::Real) = $bfn(float(nu), k, float(x))
+        $bfn(nu::AbstractFloat, k::Integer, x::AbstractFloat) = $bfn(float(nu), k, complex(x))
         function $bfn(nu::Real, k::Integer, z::Complex)
             Tf = promote_type(float(typeof(nu)),float(typeof(real(z))))
             $bfn(Tf(nu), k, Complex{Tf}(z))
@@ -596,6 +607,9 @@ for bfn in (:besselh, :besselhx)
     end
 end
 
+besselh(nu::Float16, k::Integer, x::Float16) = Complex{Float16}(besselh(Float32(nu), k, Float32(x)))
+besselh(nu::Float32, k::Integer, x::Float32) = Complex{Float32}(besselh(Float64(nu), k, Float64(x)))
+
 """
     besselj0(x)
 

test/runtests.jl

@@ -260,8 +260,10 @@ end
         @test besselh(Float16(3),1,Float16(3)) ≈ true_h133
         @test besselh(3,2,3) ≈ conj(true_h133)
         @test besselh(-3,2,3) ≈ -conj(true_h133)
+        @test besselh(1,1,0) ≈ complex(0,-Inf)
+        @test besselh(1,2,0) ≈ complex(0,Inf)
         @testset "Error throwing" begin
-            @test_throws AmosException besselh(1,0)
+            @test_throws AmosException besselh(1,5,0)
             @test_throws MethodError besselh(1,big(1.0))
             @test_throws MethodError besselh(1,complex(big(1.0)))
             @test_throws MethodError besselhx(1,big(1.0))