From 3d36d38d3d0584d85467af9b3233d4abb29cac07 Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@mit.edu>
Date: Thu, 5 Jun 2014 15:52:41 -0400
Subject: [PATCH 1/2] =?UTF-8?q?export=20@evalpoly=20macro=20(n=C3=A9e=20@c?=
 =?UTF-8?q?horner)?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 base/exports.jl       |  1 +
 base/math.jl          | 23 ++++++++++++-----------
 base/special/gamma.jl |  4 ++--
 doc/stdlib/base.rst   |  8 ++++++++
 4 files changed, 23 insertions(+), 13 deletions(-)

diff --git a/base/exports.jl b/base/exports.jl
index 7f16e442d7059..ab762defb5d74 100644
--- a/base/exports.jl
+++ b/base/exports.jl
@@ -268,6 +268,7 @@ export
     At_rdiv_Bt,
 
 # scalar math
+    @evalpoly,
     abs,
     abs2,
     acos,
diff --git a/base/math.jl b/base/math.jl
index 61aa80f7000fd..e347a28f6abb6 100644
--- a/base/math.jl
+++ b/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)
diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 08770004039f5..3c24b4b05b962 100644
--- a/base/special/gamma.jl
+++ b/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
diff --git a/doc/stdlib/base.rst b/doc/stdlib/base.rst
index a388ce71a6845..f4b9efdd61386 100644
--- a/doc/stdlib/base.rst
+++ b/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 algorithm due to Knuth.
+
 Data Formats
 ------------
 

From 9d9176a31b642f8cfbf53601c913eddda9642ea7 Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@mit.edu>
Date: Thu, 5 Jun 2014 17:29:45 -0400
Subject: [PATCH 2/2] fix doc, add NEWS for @evalpoly

---
 NEWS.md             | 2 ++
 doc/stdlib/base.rst | 2 +-
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/NEWS.md b/NEWS.md
index 661366014f116..03eb323389e7f 100644
--- a/NEWS.md
+++ b/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
 ------------------
diff --git a/doc/stdlib/base.rst b/doc/stdlib/base.rst
index f4b9efdd61386..7cf50a0a6f7e5 100644
--- a/doc/stdlib/base.rst
+++ b/doc/stdlib/base.rst
@@ -3227,7 +3227,7 @@ Mathematical Functions
    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 algorithm due to Knuth.
+   complex ``z``, a more efficient Goertzel-like algorithm.
 
 Data Formats
 ------------