From a736cb99b12785b91b7cdde49f977acfa2076f7a Mon Sep 17 00:00:00 2001
From: Jeff Bezanson <jeff.bezanson@gmail.com>
Date: Mon, 7 Oct 2013 11:30:56 -0400
Subject: [PATCH] use jn, jnf, yn, ynf for bessel functions of integer order
 and real argument also add missing domain errors in bessely fixes #4187

---
 base/math.jl | 52 ++++++++++++++++++++++++++++++++++++----------------
 base/mpfr.jl |  9 +++++++++
 2 files changed, 45 insertions(+), 16 deletions(-)

diff --git a/base/math.jl b/base/math.jl
index 7a8fcde320400..db090b64160da 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -396,13 +396,22 @@ modf(x) = rem(x,one(x)), trunc(x)
 
 # special functions
 
-for (jy,nu) in (("j",0), ("j",1), ("y",0), ("y",1))
+for jy in ("j","y"), nu in (0,1)
     jynu = Expr(:quote, symbol(string(jy,nu)))
     jynuf = Expr(:quote, symbol(string(jy,nu,"f")))
     bjynu = symbol(string("bessel",jy,nu))
+    if jy == "y"
+        @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)
+        end
+    else
+        @eval begin
+            $bjynu(x::Float64) = ccall(($jynu,libm),  Float64, (Float64,), x)
+            $bjynu(x::Float32) = ccall(($jynuf,libm), Float32, (Float32,), x)
+        end
+    end
     @eval begin
-        $bjynu(x::Float64) = ccall(($jynu,libm),  Float64, (Float64,), x)
-        $bjynu(x::Float32) = ccall(($jynuf,libm), Float32, (Float32,), x)
         $bjynu(x::Real) = $bjynu(float(x))
         $bjynu(x::Complex) = $(symbol(string("bessel",jy)))($nu,x)
         @vectorize_1arg Number $bjynu
@@ -539,18 +548,10 @@ function besselj(nu::Float64, z::Complex128)
 end
 
 function besselj(nu::Integer, x::FloatingPoint)
-    if x == 0
-        return (nu == 0) ? one(x) : zero(x)
-    end
-    if nu < 0
-        nu = -nu
-        x = -x
-    end
-    ans = _besselj(float64(nu), complex128(abs(x)))
-    if (x < 0) && isodd(nu)
-        ans = -ans
-    end
-    oftype(x, real(ans))
+    return oftype(x, ccall((:jn, libm), Float64, (Cint, Float64), nu, x))
+end
+function besselj(nu::Integer, x::Float32)
+    return ccall((:jnf, libm), Float32, (Cint, Float32), nu, x)
 end
 
 besselk(nu::Float64, z::Complex128) = _besselk(abs(nu), z)
@@ -581,7 +582,11 @@ end
 @vectorize_2arg Number besseli
 
 function besselj(nu::FloatingPoint, x::FloatingPoint)
-    if x < 0 && !isinteger(nu)
+    if isinteger(nu)
+        if typemin(Int32) <= nu <= typemax(Int32)
+            return besselj(int(nu), x)
+        end
+    elseif x < 0
         throw(DomainError())
     end
     oftype(x, real(besselj(float64(nu), complex128(x))))
@@ -610,8 +615,23 @@ function bessely(nu::Real, x::FloatingPoint)
     if x < 0
         throw(DomainError())
     end
+    if isinteger(nu) && typemin(Int32) <= nu <= typemax(Int32)
+        return bessely(int(nu), x)
+    end
     oftype(x, real(bessely(float64(nu), complex128(x))))
 end
+function bessely(nu::Integer, x::FloatingPoint)
+    if x < 0
+        throw(DomainError())
+    end
+    return oftype(x, ccall((:yn, libm), Float64, (Cint, Float64), nu, x))
+end
+function bessely(nu::Integer, x::Float32)
+    if x < 0
+        throw(DomainError())
+    end
+    return ccall((:ynf, libm), Float32, (Cint, Float32), nu, x)
+end
 @vectorize_2arg Number bessely
 
 hankelh1(nu, z) = besselh(nu, 1, z)
diff --git a/base/mpfr.jl b/base/mpfr.jl
index 51fc6a3bc3fd7..b84aab1f4ebe3 100644
--- a/base/mpfr.jl
+++ b/base/mpfr.jl
@@ -443,18 +443,27 @@ function besselj(n::Integer, x::BigFloat)
 end
 
 function bessely0(x::BigFloat)
+    if x < 0
+        throw(DomainError())
+    end
     z = BigFloat()
     ccall((:mpfr_y0, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
     return z
 end
 
 function bessely1(x::BigFloat)
+    if x < 0
+        throw(DomainError())
+    end
     z = BigFloat()
     ccall((:mpfr_y1, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
     return z
 end
 
 function bessely(n::Integer, x::BigFloat)
+    if x < 0
+        throw(DomainError())
+    end
     z = BigFloat()
     ccall((:mpfr_yn, :libmpfr), Int32, (Ptr{BigFloat}, Clong, Ptr{BigFloat}, Int32), &z, n, &x, ROUNDING_MODE[end])
     return z