Added Float16 fix in bessel.jl (#144)

Author: Sumegh-git

src/bessel.jl

@@ -131,8 +131,9 @@ for afn in (:airyai, :airyaiprime, :airybi, :airybiprime,
             :airyaix, :airyaiprimex, :airybix, :airybiprimex)
     @eval begin
         $afn(z::Complex) = $afn(float(z))
-        $afn(z::Complex{<:AbstractFloat}) = throw(MethodError($afn,(z,)))
+        $afn(z::Complex{Float16}) = Complex{Float16}($afn(Complex{Float32}(z)))
         $afn(z::Complex{Float32}) = Complex{Float32}($afn(Complex{Float64}(z)))
+        $afn(z::Complex{<:AbstractFloat}) = throw(MethodError($afn,(z,)))
     end
     if afn in (:airyaix, :airyaiprimex)
         @eval $afn(x::Real) = x < 0 ? throw(DomainError(x, "`x` must be nonnegative.")) : real($afn(complex(float(x))))
@@ -158,11 +159,13 @@ for jy in ("j","y"), nu in (0,1)
         @eval begin
             $bjynu(x::Float64) = nan_dom_err(ccall(($jynu,libm),  Float64, (Float64,), x), x)
             $bjynu(x::Float32) = nan_dom_err(ccall(($jynuf,libm), Float32, (Float32,), x), x)
+            $bjynu(x::Float16) = Float16($bjynu(Float32(x)))
         end
     else
         @eval begin
             $bjynu(x::Float64) = ccall(($jynu,libm),  Float64, (Float64,), x)
             $bjynu(x::Float32) = ccall(($jynuf,libm), Float32, (Float32,), x)
+            $bjynu(x::Float16) = Float16($bjynu(Float32(x)))
         end
     end
     @eval begin
@@ -499,8 +502,9 @@ for f in ("i", "ix", "j", "jx", "k", "kx", "y", "yx")
             Tf = promote_type(float(typeof(nu)),float(typeof(real(z))))
             $bfn(Tf(nu), Complex{Tf}(z))
         end
-        $bfn(k::T, z::Complex{T}) where {T<:AbstractFloat} = throw(MethodError($bfn,(k,z)))
+        $bfn(nu::Float16, x::Complex{Float16}) = Complex{Float16}($bfn(Float32(nu), Complex{Float32}(x)))        
         $bfn(nu::Float32, x::Complex{Float32}) = Complex{Float32}($bfn(Float64(nu), Complex{Float64}(x)))
+        $bfn(k::T, z::Complex{T}) where {T<:AbstractFloat} = throw(MethodError($bfn,(k,z)))
     end
 end
 
@@ -515,9 +519,9 @@ for bfn in (:besselh, :besselhx)
             Tf = promote_type(float(typeof(nu)),float(typeof(real(z))))
             $bfn(Tf(nu), k, Complex{Tf}(z))
         end
-
-        $bfn(nu::T, k::Integer, z::Complex{T}) where {T<:AbstractFloat} = throw(MethodError($bfn,(nu,k,z)))
+        $bfn(nu::Float16, k::Integer, x::Complex{Float16}) = Complex{Float16}($bfn(Float32(nu), k, Complex{Float32}(x)))
         $bfn(nu::Float32, k::Integer, x::Complex{Float32}) = Complex{Float32}($bfn(Float64(nu), k, Complex{Float64}(x)))
+        $bfn(nu::T, k::Integer, z::Complex{T}) where {T<:AbstractFloat} = throw(MethodError($bfn,(nu,k,z)))
     end
 end
 

test/runtests.jl

@@ -93,7 +93,7 @@ end
 	@test_throws MethodError ellipe(BigFloat(-1))
 	@test_throws DomainError ellipe(Float16(2.0))
 	@test_throws DomainError ellipe(Float32(2.5))
-end    
+end 
 @testset "sine and cosine integrals" begin
     # Computed via wolframalpha.com: SinIntegral[SetPrecision[Table[x,{x, 1,20,1}],20]] and CosIntegral[SetPrecision[Table[x,{x, 1,20,1}],20]]
     sinintvals = [0.9460830703671830149, 1.605412976802694849, 1.848652527999468256, 1.75820313894905306, 1.54993124494467414, 1.4246875512805065, 1.4545966142480936, 1.5741868217069421, 1.665040075829602, 1.658347594218874, 1.578306806945727416, 1.504971241526373371, 1.499361722862824564, 1.556211050077665054, 1.618194443708368739, 1.631302268270032886, 1.590136415870701122, 1.536608096861185462, 1.518630031769363932, 1.548241701043439840]
@@ -135,13 +135,13 @@ end
     @test_throws AmosException airyai(200im)
     @test_throws AmosException airybi(200)
 
-    for T in [Float32, Float64, Complex{Float32},Complex{Float64}]
+    for T in [Float16, Float32, Float64,Complex{Float16}, Complex{Float32},Complex{Float64}]
         @test airyai(T(1.8)) ≈ 0.0470362168668458052247
         @test airyaiprime(T(1.8)) ≈ -0.0685247801186109345638
         @test airybi(T(1.8)) ≈ 2.595869356743906290060
         @test airybiprime(T(1.8)) ≈ 2.98554005084659907283
     end
-    for T in [Complex{Float32}, Complex{Float64}]
+    for T in [Complex{Float16}, Complex{Float32}, Complex{Float64}]
         z = convert(T,1.8 + 1.0im)
         @test airyaix(z) ≈ airyai(z) * exp(2/3 * z * sqrt(z))
         @test airyaiprimex(z) ≈ airyaiprime(z) * exp(2/3 * z * sqrt(z))
@@ -172,6 +172,8 @@ end
     @testset "$z, $o" for (z, o, f) in bessel_funcs
         @test z(Float32(2.0)) ≈ z(Float64(2.0))
         @test o(Float32(2.0)) ≈ o(Float64(2.0))
+        @test z(Float16(2.0)) ≈ z(Float64(2.0))
+        @test o(Float16(2.0)) ≈ o(Float64(2.0))
         @test z(2) ≈ z(2.0)
         @test o(2) ≈ o(2.0)
         @test z(2.0 + im) ≈ f(0, 2.0 + im)
@@ -187,6 +189,8 @@ end
         true_h133 = 0.30906272225525164362 - 0.53854161610503161800im
         @test besselh(3,1,3) ≈ true_h133
         @test besselh(-3,1,3) ≈ -true_h133
+        @test besselh(Float32(3),1,Float32(3)) ≈ true_h133	
+        @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)
         @testset "Error throwing" begin
@@ -204,12 +208,15 @@ end
         @test besseli(3,-3) ≈ -true_i33
         @test besseli(-3,-3) ≈ -true_i33
         @test besseli(Float32(-3),Complex{Float32}(-3,0)) ≈ -true_i33
+        @test besseli(Float16(-3),Complex{Float16}(-3,0)) ≈ -true_i33
         true_im3p1_3 = 0.84371226532586351965
         @test besseli(-3.1,3) ≈ true_im3p1_3
         for i in [-5 -3 -1 1 3 5]
             @test besseli(i,0) == 0.0
             @test besseli(i,Float32(0)) == 0
             @test besseli(i,Complex{Float32}(0)) == 0
+            @test besseli(i,Float16(0)) == 0
+            @test besseli(i,Complex{Float16}(0)) == 0
         end
         @testset "Error throwing" begin
             @test_throws AmosException besseli(1,1000)
@@ -225,7 +232,7 @@ end
         for i in [-5 -3 -1 1 3 5]
             @test besselj(i,0) == 0
             @test besselj(i,Float32(0)) == 0
-            @test besselj(i,Complex{Float32}(0)) == 0
+            @test besselj(i,Complex{Float32}(0)) == 0.0
         end
 
         j33 = besselj(3,3.)
@@ -257,6 +264,7 @@ end
 
         true_jm3p1_3 = -0.45024252862270713882
         @test besselj(-3.1,3) ≈ true_jm3p1_3
+        @test besselj(Float16(-3.1),Float16(3)) ≈ true_jm3p1_3
 
         @testset "Error throwing" begin
             @test_throws DomainError    besselj(0.1, -0.4)
@@ -268,6 +276,8 @@ end
     @testset "besselk" begin
         true_k33 = 0.12217037575718356792
         @test besselk(3,3) ≈ true_k33
+        @test besselk(Float32(3),Float32(3)) ≈ true_k33
+        @test besselk(Float16(3),Float16(3)) ≈ true_k33
         @test besselk(-3,3) ≈ true_k33
         true_k3m3 = -0.1221703757571835679 - 3.0151549516807985776im
         @test besselk(3,complex(-3)) ≈ true_k3m3
@@ -311,7 +321,7 @@ end
     end
 
     @testset "besselhx" begin
-        for elty in [Complex{Float32},Complex{Float64}]
+        for elty in [Complex{Float16},Complex{Float32},Complex{Float64}]
             z = convert(elty, 1.0 + 1.9im)
             @test besselhx(1.0, 1, z) ≈ convert(elty,-0.5949634147786144 - 0.18451272807835967im)
             @test besselhx(Float32(1.0), 1, z) ≈ convert(elty,-0.5949634147786144 - 0.18451272807835967im)
@@ -341,6 +351,10 @@ end
             @test besselix(i,Float32(0)) == 0
             @test besseljx(i,Complex{Float32}(0)) == 0
             @test besselix(i,Complex{Float32}(0)) == 0
+            @test besseljx(i,Float16(0)) == 0
+            @test besselix(i,Float16(0)) == 0
+            @test besseljx(i,Complex{Float16}(0)) == 0
+            @test besselix(i,Complex{Float16}(0)) == 0
         end
         @testset "Error throwing" begin
             @test_throws AmosException hankelh1x(1, 0)
@@ -358,6 +372,7 @@ end
         @testset "$f" for f in (besselj,bessely,besseli,besselk,hankelh1,hankelh2)
             @test f(0,1) ≈ f(0,Complex{Float64}(1))
             @test f(0,1) ≈ f(0,Complex{Float32}(1))
+            @test f(0,1) ≈ f(0,Complex{Float16}(1))
         end
     end
 end
@@ -596,7 +611,6 @@ end
 end
 
 @test sprint(showerror, AmosException(1)) == "AmosException with id 1: input error."
-
 # Used to check method existence below
 struct NotAFloat <: AbstractFloat
 end
@@ -664,7 +678,6 @@ end
         @test lgamma(-Inf*im) === -Inf - Inf*im
         @test lgamma(Inf + Inf*im) === lgamma(NaN + 0im) === lgamma(NaN*im) === NaN + NaN*im
     end
-
     @testset "Other float types" begin
         let x = one(Float16)
             @test gamma(x) ≈ one(Float16)