fix evalpoly type instability and 0-length case

Author: stevengj

base/math.jl

@@ -23,7 +23,7 @@ import .Base: log, exp, sin, cos, tan, sinh, cosh, tanh, asin,
 using .Base: sign_mask, exponent_mask, exponent_one,
             exponent_half, uinttype, significand_mask,
             significand_bits, exponent_bits, exponent_bias,
-            exponent_max, exponent_raw_max, clamp, clamp!
+            exponent_max, exponent_raw_max, clamp, clamp!, reduce_empty_iter
 
 using Core.Intrinsics: sqrt_llvm
 
@@ -103,17 +103,20 @@ function evalpoly(x, p::Tuple)
         _evalpoly(x, p)
     end
 end
+evalpoly(x, ::Tuple{}) = zero(one(x)) # dimensionless zero, i.e. 0 * x^0
 
 evalpoly(x, p::AbstractVector) = _evalpoly(x, p)
 
 function _evalpoly(x, p)
     Base.require_one_based_indexing(p)
     N = length(p)
-    ex = p[end]
+    @inbounds p0 = N == 0 ? reduce_empty_iter(+, p) : p[N]
+    s = oftype(one(x) * p0, p0)
+    N <= 1 && return s
     for i in N-1:-1:1
-        ex = muladd(x, ex, p[i])
+        @inbounds s = muladd(x, s, p[i])
     end
-    ex
+    return s
 end
 
 function evalpoly(z::Complex, p::Tuple)
@@ -142,16 +145,17 @@ function evalpoly(z::Complex, p::Tuple)
     end
 end
 evalpoly(z::Complex, p::Tuple{<:Any}) = p[1]
-
+evalpoly(z::Complex, ::Tuple{}) = zero(one(z)) # dimensionless zero, i.e. 0 * z^0
 
 evalpoly(z::Complex, p::AbstractVector) = _evalpoly(z, p)
 
 function _evalpoly(z::Complex, p)
     Base.require_one_based_indexing(p)
-    length(p) == 1 && return p[1]
     N = length(p)
-    a = p[end]
-    b = p[end-1]
+    @inbounds p0 = N == 0 ? reduce_empty_iter(+, p) : p[N]
+    N <= 1 && return oftype(one(z) * p0, p0)
+    a = p0
+    @inbounds b = p[N-1]
 
     x = real(z)
     y = imag(z)
@@ -160,7 +164,7 @@ function _evalpoly(z::Complex, p)
     for i in N-2:-1:1
         ai = a
         a = muladd(r, ai, b)
-        b = muladd(-s, ai, p[i])
+        @inbounds b = muladd(-s, ai, p[i])
     end
     ai = a
     muladd(ai, z, b)

test/math.jl

@@ -715,6 +715,8 @@ end
         @test evalpoly(x, (p1, p2, p3)) == evpm
         @test evalpoly(x, [p1, p2, p3]) == evpm
     end
+    @test evalpoly(1.0f0, ()) === 0.0f0 # issue #56699
+    @test @inferred(evalpoly(1.0f0, [2])) === 2.0f0 # type-stability
 end
 
 @testset "evalpoly complex" begin
@@ -726,6 +728,8 @@ end
     end
     @test evalpoly(1+im, (2,)) == 2
     @test evalpoly(1+im, [2,]) == 2
+    @test evalpoly(1.0f0+im, ()) === 0.0f0+0im # issue #56699
+    @test @inferred(evalpoly(1.0f0+im, [2])) === 2.0f0+0im # type-stability
 end
 
 @testset "cis" begin