Make pochhammer more robust

src/specialfunctions.jl

@@ -49,8 +49,8 @@ pochhammer(x::AbstractArray{T,1},n::Integer) where {T<:Number} = [pochhammer(x[i
 pochhammer(x::AbstractArray{T,2},n::Integer) where {T<:Number} = [pochhammer(x[i,j],n) for i=1:size(x,1),j=1:size(x,2)]
 pochhammer(x::AbstractArray{T},n::Integer) where {T<:Number} = reshape([ pochhammer(x[i],n) for i in eachindex(x) ], size(x))
 
-pochhammer(x::Number,n::Number) = gamma(x+n)/gamma(x)
-pochhammer(x::AbstractArray{T},n::Number) where {T<:Number} = gamma(x+n)./gamma(x)
+pochhammer(x::Number,n::Number) = isinteger(n) ? pochhammer(x,Int(n)) : ogamma(x)/ogamma(x+n)
+pochhammer(x::AbstractArray{T},n::Number) where {T<:Number} = isinteger(n) ? pochhammer(x,Int(n)) : ogamma.(x)./ogamma.(x.+n)
 
 function pochhammer(x::Number,n::UnitRange{T}) where T<:Real
     ret = Vector{promote_type(typeof(x),T)}(length(n))
@@ -61,6 +61,14 @@ function pochhammer(x::Number,n::UnitRange{T}) where T<:Real
     ret
 end
 
+function ogamma(x::Number)
+    if isinteger(x) && x<0
+        zero(float(x))
+    else
+        inv(gamma(x))
+    end
+end
+
 """
 Stirling's asymptotic series for ``\\Gamma(z)``.
 """

test/runtests.jl

@@ -28,3 +28,5 @@ include("butterflytests.jl")
 include("sphericalharmonics/sphericalharmonictests.jl")
 
 include("toeplitztests.jl")
+
+include("specialfunctionstests.jl")

test/specialfunctionstests.jl

@@ -0,0 +1,19 @@
+using FastTransforms
+using Compat.Test
+import FastTransforms.pochhammer
+
+const sqrtpi = 1.772453850905516027298
+
+@testset "Pochhammer" begin
+    @test pochhammer(2,3) == 24
+    @test pochhammer(0.5,3) == 0.5*1.5*2.5
+    @test pochhammer(0.5,0.5) == 1/sqrtpi
+    @test pochhammer(0,1) == 0
+    @test pochhammer(-1,2) == 0
+    @test pochhammer(-5,3) == -60
+    @test pochhammer(-1,-0.5) == 0
+    @test 1.0/pochhammer(-0.5,-0.5) == 0
+    @test pochhammer(-1+0im,-1) == -2
+    @test pochhammer([2,-5,-1],3) == [24,-60,0]
+    @test pochhammer([-1,-0.5],-0.5) == [0,-Inf]
+end