export @evalpoly macro (née @chorner)

NEWS.md

@@ -306,6 +306,8 @@ Library improvements
 
   * New functions `randsubseq` and `randsubseq!` to create a random subsequence of an array ([#6726])
 
+  * New macro `@evalpoly` for efficient inline evaluation of polynomials ([#7146]).
+
 
 Build improvements
 ------------------

base/exports.jl

@@ -268,6 +268,7 @@ export
     At_rdiv_Bt,
 
 # scalar math
+    @evalpoly,
     abs,
     abs2,
     acos,

base/math.jl

@@ -43,16 +43,16 @@ clamp{T}(x::AbstractArray{T}, lo, hi) =
 macro horner(x, p...)
     ex = esc(p[end])
     for i = length(p)-1:-1:1
-        ex = :($(esc(p[i])) + $(esc(x)) * $ex)
+        ex = :($(esc(p[i])) + t * $ex)
     end
-    ex
+    Expr(:block, :(t = $(esc(x))), ex)
 end
 
-# Evaluate p[1] + z*p[2] + z^2*p[3] + ... + z^(n-1)*p[n], assuming
-# the p[k] are *real* coefficients.   This uses Horner's method if z
-# is real, but for complex z it uses a more efficient algorithm described
-# in Knuth, TAOCP vol. 2, section 4.6.4, equation (3).
-macro chorner(z, p...)
+# Evaluate p[1] + z*p[2] + z^2*p[3] + ... + z^(n-1)*p[n].  This uses
+# Horner's method if z is real, but for complex z it uses a more
+# efficient algorithm described in Knuth, TAOCP vol. 2, section 4.6.4,
+# equation (3).
+macro evalpoly(z, p...)
     a = :($(esc(p[end])))
     b = :($(esc(p[end-1])))
     as = {}
@@ -65,14 +65,15 @@ macro chorner(z, p...)
     ai = :a0
     push!(as, :($ai = $a))
     C = Expr(:block,
-             :(x = real($(esc(z)))),
-             :(y = imag($(esc(z)))),
+             :(t = $(esc(z))),
+             :(x = real(t)),
+             :(y = imag(t)),
              :(r = x + x),
              :(s = x*x + y*y),
              as...,
-             :($ai * $(esc(z)) + $b))
+             :($ai * t + $b))
     R = Expr(:macrocall, symbol("@horner"), esc(z), p...)
-    :(isa($(esc(z)), Real) ? $R : $C)
+    :(isa($(esc(z)), Complex) ? $C : $R)
 end
 
 rad2deg(z::Real) = oftype(z, 57.29577951308232*z)

base/special/gamma.jl

@@ -91,7 +91,7 @@ function digamma(z::Union(Float64,Complex{Float64}))
     ψ += log(z) - 0.5*t
     t *= t # 1/z^2
     # the coefficients here are float64(bernoulli[2:9] .// (2*(1:8)))
-    ψ -= t * @chorner(t,0.08333333333333333,-0.008333333333333333,0.003968253968253968,-0.004166666666666667,0.007575757575757576,-0.021092796092796094,0.08333333333333333,-0.4432598039215686)
+    ψ -= t * @evalpoly(t,0.08333333333333333,-0.008333333333333333,0.003968253968253968,-0.004166666666666667,0.007575757575757576,-0.021092796092796094,0.08333333333333333,-0.4432598039215686)
 end
 
 function trigamma(z::Union(Float64,Complex{Float64}))
@@ -114,7 +114,7 @@ function trigamma(z::Union(Float64,Complex{Float64}))
     w = t * t # 1/z^2
     ψ += t + 0.5*w
     # the coefficients here are float64(bernoulli[2:9])
-    ψ += t*w * @chorner(w,0.16666666666666666,-0.03333333333333333,0.023809523809523808,-0.03333333333333333,0.07575757575757576,-0.2531135531135531,1.1666666666666667,-7.092156862745098)
+    ψ += t*w * @evalpoly(w,0.16666666666666666,-0.03333333333333333,0.023809523809523808,-0.03333333333333333,0.07575757575757576,-0.2531135531135531,1.1666666666666667,-7.092156862745098)
 end
 
 signflip(m::Number, z) = (-1+0im)^m * z

doc/stdlib/base.rst

@@ -3221,6 +3221,14 @@ Mathematical Functions
 
    Multiply ``x`` and ``y``, giving the result as a larger type.
 
+.. function:: @evalpoly(z, c...)
+
+   Evaluate the polynomial :math:`\sum_k c[k] z^{k-1}` for the
+   coefficients ``c[1]``, ``c[2]``, ...; that is, the coefficients are
+   given in ascending order by power of ``z``.  This macro expands to
+   efficient inline code that uses either Horner's method or, for
+   complex ``z``, a more efficient Goertzel-like algorithm.
+
 Data Formats
 ------------