Merge pull request #7125 from stevengj/polygamma

RFC: add complex polygamma and Hurwitz zeta functions

NEWS.md

@@ -161,6 +161,16 @@ Library improvements
 
     * Broadcasting `.//` is now included ([#7094]).
 
+    * `prevfloat` and `nextfloat` now saturate at -Inf and Inf,
+      respectively, and have otherwise been fixed to follow the IEEE-754
+      standard functions `nextDown` and `nextUp` ([#5025]).
+
+    * New function `widen` for widening numeric types and values, and `widemul`
+      for multiplying to a larger type ([#6169]).
+
+    * `polygamma`, `digamma`, and `trigamma` now accept complex
+      arguments, and `zeta(s, z)` now provides the Hurwitz zeta ([#7125]).
+
   * `String` improvements
 
     * Triple-quoted regex strings, `r"""..."""` ([#4934]).
@@ -259,10 +269,6 @@ Library improvements
   * New function `deleteat!` deletes a specified index or indices and
     returns the updated collection
 
-  * `prevfloat` and `nextfloat` now saturate at -Inf and Inf, respectively, and
-    have otherwise been fixed to follow the IEEE-754 standard functions `nextDown`
-    and `nextUp` ([#5025]).
-
   * The `setenv` function for external processes now accepts a `dir` keyword
     argument for specifying the directory to start the child process in ([#4888]).
 
@@ -272,9 +278,6 @@ Library improvements
   * Ranges and arrays with the same elements are now unequal. This allows hashing
     and comparing ranges to be faster. ([#5778])
 
-  * New function `widen` for widening numeric types and values, and `widemul`
-    for multiplying to a larger type ([#6169])
-
   * Broadcasting now works on arbitrary `AbstractArrays` ([#5387])
 
   * Reduction functions that accept a pre-allocated output array, including

base/math.jl

@@ -48,6 +48,33 @@ macro horner(x, p...)
     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...)
+    a = :($(esc(p[end])))
+    b = :($(esc(p[end-1])))
+    as = {}
+    for i = length(p)-2:-1:1
+        ai = symbol(string("a", i))
+        push!(as, :($ai = $a))
+        a = :($b + r*$ai)
+        b = :($(esc(p[i])) - s * $ai)
+    end
+    ai = :a0
+    push!(as, :($ai = $a))
+    C = Expr(:block,
+             :(x = real($(esc(z)))),
+             :(y = imag($(esc(z)))),
+             :(r = x + x),
+             :(s = x*x + y*y),
+             as...,
+             :($ai * $(esc(z)) + $b))
+    R = Expr(:macrocall, symbol("@horner"), esc(z), p...)
+    :(isa($(esc(z)), Real) ? $R : $C)
+end
+
 rad2deg(z::Real) = oftype(z, 57.29577951308232*z)
 deg2rad(z::Real) = oftype(z, 0.017453292519943295*z)
 rad2deg(z::Integer) = rad2deg(float(z))