fix ellipk and ellipe (#499)

* fix ellipk and ellipe

* use ===

* Apply suggestion from @inkydragon

---------

Co-authored-by: Chengyu Han <cyhan.dev@outlook.com>

Author: ogauthe

src/ellip.jl

@@ -36,6 +36,9 @@ As suggested in this paper, the domain is restricted to ``(-\infty,1]``.
 ellipk(m::Real) = _ellipk(float(m))
 
 function _ellipk(m::Float64)
+    isnan(m) && return NaN
+    (m === -Inf) && return 0.0
+
     flag_is_m_neg = m < 0.0
     if flag_is_m_neg
         x               = m / (m-1)         #dealing with negative args
@@ -47,7 +50,7 @@ function _ellipk(m::Float64)
         return Float64(halfπ)
 
     elseif x == 1.0
-        return Inf
+        return flag_is_m_neg ? 0.0 : Inf
 
     elseif x > 1.0
         throw(DomainError(m, "`m` must lie between -Inf and 1 ---- Domain: (-Inf,1.0]"))
@@ -206,6 +209,9 @@ As suggested in this paper, the domain is restricted to ``(-\infty,1]``.
 ellipe(m::Real) = _ellipe(float(m))
 
 function _ellipe(m::Float64)
+    isnan(m) && return NaN
+    (m === -Inf) && return Inf
+
     flag_is_m_neg = m < 0.0
     if flag_is_m_neg
         x               = m / (m-1)         #dealing with negative args
@@ -216,7 +222,7 @@ function _ellipe(m::Float64)
     if x == 0.0
         return Float64(halfπ)
     elseif x == 1.0
-        return 1.0
+        return flag_is_m_neg ? Inf : 1.0
 
     elseif x > 1.0
         throw(DomainError(m,"`m` must lie between -Inf and 1 ---- Domain : (-inf,1.0]"))

test/ellip.jl

@@ -17,6 +17,9 @@
         @test ellipk(1.0) == Inf
         @test ellipk(Float16(0.92)) ≈ 2.683551406315229344  rtol=2*eps(Float16)
         @test ellipk(Float32(0.92)) ≈ 2.683551406315229344  rtol=2*eps(Float32)
+        @test isnan(ellipk(NaN))
+        @test ellipk(-Inf) == 0.0
+        @test ellipk(-1e30) ≈ 0.0                          atol=1e-13
         @test_throws MethodError ellipk(BigFloat(0.5))
         @test_throws DomainError ellipk(1.1)
     end
@@ -38,6 +41,9 @@
         @test ellipe(1.0)   ≈ 1.00                          rtol=2*eps()
         @test ellipe(Float16(0.865)) ≈ 1.1322436887003925   rtol=2*eps(Float16)
         @test ellipe(Float32(0.865)) ≈ 1.1322436887003925   rtol=2*eps(Float32)
+        @test isnan(ellipe(NaN))
+        @test ellipe(-Inf) == Inf
+        @test ellipe(-1e16) > ellipe(-1e15)
         @test_throws MethodError ellipe(BigFloat(-1))
         @test_throws DomainError ellipe(1.2)
     end