Merge pull request #97 from yha/hypergeom-nan-checks

Test hypergeometric arguments for NaN values

src/manual_wrappers.jl

@@ -231,6 +231,9 @@ end
 #(c) 2013 Jiahao Chen <[email protected]>
 export hypergeom, hypergeom_e
 
+# NaN values not handled correctly by GSL, see GSL.jl issue #96
+@inline _hypergeom_any_nan(a,b,x) = any(isnan,a) || any(isnan,b) || isnan(x)
+
 """
     hypergeom(a, b, x::Float64) -> Float64
 
@@ -243,6 +246,7 @@ and length-0 `a` and/or `b` may be input as simply `[]`.
 Supported values of ``(p, q)`` are ``(0, 0)``, ``(0, 1)``, ``(1, 1)``, ``(2, 0)`` and ``(2, 1)``.
 """
 function hypergeom(a, b, x)
+    _hypergeom_any_nan(a,b,x) && return NaN
     n = length(a), length(b)
     if n == (0, 0)
         exp(x)
@@ -271,6 +275,7 @@ and length-0 `a` and/or `b` may be input as simply `[]`.
 Supported values of ``(p, q)`` are ``(0, 0)``, ``(0, 1)``, ``(1, 1)``, ``(2, 0)`` and ``(2, 1)``.
 """
 function hypergeom_e(a, b, x)
+    _hypergeom_any_nan(a,b,x) && return gsl_sf_result(NaN, NaN)
     n = length(a), length(b)
     if n == (0, 0)
         sf_exp_err_e(x,0.0)

test/hypergeom.jl

@@ -34,5 +34,15 @@ using SpecialFunctions
     #Kummer's theorem
     res = hypergeom_e([a, b], 1+a-b, -1.0)
     @test isapprox(res.val, (gamma(1+a-b)*gamma(1+a/2)/(gamma(1+a)*gamma(1+a/2-b))), atol=res.err)
-
+
+    #NaN handling
+    for (h, condition) in ((hypergeom, isnan), (hypergeom_e, r->isnan(r.val)))
+        @test all(condition,
+                (h((),(),NaN),
+                 h((),NaN,x), h((),c,NaN),
+                 h(NaN,c,x), h(a,NaN,x), h(a,c,NaN),
+                 h((NaN,b),(),x), h((a,NaN),(),x), h((a,b),(),NaN),
+                 h((NaN,b),c,x), h((a,NaN),c,x), h((a,b),NaN,x), h((a,b),c,NaN)
+                ))
+    end
 end