From 055bbe4aaeb930cfb19394b1848fd7300f1b0a86 Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@mit.edu>
Date: Thu, 5 Jun 2014 10:28:10 -0400
Subject: [PATCH] fixed type-stability of chorner and polygamma (workaround for
 #7060)

---
 base/math.jl          | 2 +-
 base/special/gamma.jl | 9 ++++++---
 test/math.jl          | 3 +++
 3 files changed, 10 insertions(+), 4 deletions(-)

diff --git a/base/math.jl b/base/math.jl
index 180bd9cc10f2f..61aa80f7000fd 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -70,7 +70,7 @@ macro chorner(z, p...)
              :(r = x + x),
              :(s = x*x + y*y),
              as...,
-             :(Complex($ai * x + $b, $ai * y)))
+             :($ai * $(esc(z)) + $b))
     R = Expr(:macrocall, symbol("@horner"), esc(z), p...)
     :(isa($(esc(z)), Real) ? $R : $C)
 end
diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 3775cc189305d..93019511d6382 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -224,6 +224,9 @@ function inv_oftype(x::Complex, y::Real)
 end
 inv_oftype(x::Real, y::Real) = oftype(x, inv(y))
 
+zeta_returntype{T}(s::Integer, z::T) = T
+zeta_returntype(s, z) = Complex128
+
 # Hurwitz zeta function, which is related to polygamma 
 # (at least for integer m > 0 and real(z) > 0) by:
 #    polygamma(m, z) = (-1)^(m+1) * gamma(m+1) * zeta(m+1, z).
@@ -237,7 +240,7 @@ inv_oftype(x::Real, y::Real) = oftype(x, inv(y))
 # (which Julia already exports).
 function zeta(s::Union(Int,Float64,Complex{Float64}),
               z::Union(Float64,Complex{Float64}))
-    ζ = isa(s, Int) ? zero(z) : complex(zero(z))
+    ζ = zero(zeta_returntype(s, z))
     z == 1 && return oftype(ζ, zeta(s))
     s == 2 && return oftype(ζ, trigamma(z))
 
@@ -245,13 +248,13 @@ function zeta(s::Union(Int,Float64,Complex{Float64}),
 
     # annoying s = Inf or NaN cases:
     if !isfinite(s)
-        (isnan(s) || isnan(z)) && return Complex(NaN, NaN)
+        (isnan(s) || isnan(z)) && return (s*z)^2 # type-stable NaN+Nan*im
         if real(s) == Inf
             z == 1 && return one(ζ)
             if x > 1 || (x >= 0.5 ? abs(z) > 1 : abs(z - round(x)) > 1)
                 return zero(ζ) # distance to poles is > 1
             end
-            x > 0 && imag(z) == 0 && imag(s) == 0 && return Complex(Inf, 0.0)
+            x > 0 && imag(z) == 0 && imag(s) == 0 && return oftype(ζ, Inf)
         end
         throw(DomainError()) # nothing clever to return
     end
diff --git a/test/math.jl b/test/math.jl
index fdfb4cd12a0e8..1b15d72544d1d 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -274,3 +274,6 @@ end
 @test polygamma(4, -0.0) == Inf == -polygamma(4, +0.0)
 @test zeta(4, +0.0) == Inf == zeta(4, -0.0)
 @test zeta(5, +0.0) == Inf == -zeta(5, -0.0)
+@test isa([digamma(x) for x in [1.0]], Vector{Float64})
+@test isa([trigamma(x) for x in [1.0]], Vector{Float64})
+@test isa([polygamma(3,x) for x in [1.0]], Vector{Float64})