Add evalpoly function (JuliaLang#32753)

NEWS.md

@@ -36,7 +36,7 @@ New library functions
 * `readdir` output is now guaranteed to be sorted. The `sort` keyword allows opting out of sorting to get names in OS-native order ([#33542]).
 * The new `only(x)` function returns the one-and-only element of a collection `x`, and throws an `ArgumentError` if `x` contains zero or multiple elements. ([#33129])
 * `takewhile` and `dropwhile` have been added to the Iterators submodule ([#33437]).
-
+* There is a now an `evalpoly` (generated) function meant to take the role of the `@evalpoly` macro. The function is just as efficient as the macro while giving added flexibility, so it should be preferred over `@evalpoly`. `evalpoly` takes a list of coefficients as a tuple, so where one might write `@evalpoly(x, p1, p2, p3)` one would instead write `evalpoly(x, (p1, p2, p3))`.
 
 Standard library changes
 ------------------------

base/exports.jl

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

base/math.jl

@@ -13,7 +13,7 @@ export sin, cos, sincos, tan, sinh, cosh, tanh, asin, acos, atan,
  cbrt, sqrt, significand,
  hypot, max, min, minmax, ldexp, frexp,
  clamp, clamp!, modf, ^, mod2pi, rem2pi,
- @evalpoly
+ @evalpoly, evalpoly
 
 import .Base: log, exp, sin, cos, tan, sinh, cosh, tanh, asin,
  acos, atan, asinh, acosh, atanh, sqrt, log2, log10,
@@ -81,18 +81,143 @@ function clamp!(x::AbstractArray, lo, hi)
  x
 end
 
+
+"""
+ evalpoly(x, p::Tuple)
+Evaluate the polynomial ``\\sum_k p[k] x^{k-1}`` for the coefficients `p[1]`, `p[2]`, ...;
+that is, the coefficients are given in ascending order by power of `x`. This function
+generates efficient code using Horner's method with loops unrolled at compile time.
+# Example
+```jldoctest
+julia> evalpoly(2, (1, 2, 3))
+17
+```
+"""
+function evalpoly(x, p::Tuple)
+ if @generated
+ N = length(p.parameters)
+ ex = :(p[end])
+ for i in N-1:-1:1
+ ex = :(muladd(x, $ex, p[$i]))
+ end
+ ex
+ else
+ _evalpoly(x, p)
+ end
+end
+
+"""
+ evalpoly(x, p::AbstractVector)
+Evaluate the polynomial ``\\sum_k p[k] x^{k-1}`` for the coefficients `p[1]`, `p[2]`, ...;
+that is, the coefficients are given in ascending order by power of `x`. This function
+uses Horner's method *without* loops unrolled at compile time. Use this method when the
+number of coefficients is not known at compile time.
+# Example
+```jldoctest
+julia> evalpoly(2, [1, 2, 3])
+17
+```
+"""
+evalpoly(x, p::AbstractVector) = _evalpoly(x, p)
+
+function _evalpoly(x, p)
+ N = length(p)
+ ex = p[end]
+ for i in N-1:-1:1
+ ex = muladd(x, ex, p[i])
+ end
+ ex
+end
+
+"""
+ evalpoly(z::Complex, p::Tuple)
+Evaluate the polynomial ``\\sum_k p[k] z^{k-1}`` for the coefficients `p[1]`, `p[2]`, ...;
+that is, the coefficients are given in ascending order by power of `z`. This function
+generates efficient code using a Goertzel-like algorithm specialized for complex arguments
+with loops unrolled at compile time.
+The Goertzel-like algorthim is described in Knuth's Art of Computer Programming,
+Volume 2: Seminumerical Algorithms, Sec. 4.6.4.
+# Example
+```jldoctest
+julia> evalpoly(2 + im, (1, 2, 3))
+14 + 14im
+```
+"""
+function evalpoly(z::Complex, p::Tuple)
+ if @generated
+ N = length(p.parameters)
+ a = :(p[end])
+ b = :(p[end-1])
+ as = []
+ for i in N-2:-1:1
+ ai = Symbol("a", i)
+ push!(as, :($ai = $a))
+ a = :(muladd(r, $ai, $b))
+ b = :(p[$i] - s * $ai)
+ end
+ ai = :a0
+ push!(as, :($ai = $a))
+ C = Expr(:block,
+ :(x = real(z)),
+ :(y = imag(z)),
+ :(r = x + x),
+ :(s = muladd(x, x, y*y)),
+ as...,
+ :(muladd($ai, z, $b)))
+ else
+ _evalpoly(z, p)
+ end
+end
+evalpoly(z::Complex, p::Tuple{<:Any}) = p[1]
+
+
+"""
+ evalpoly(z::Complex, p::AbstractVector)
+Evaluate the polynomial ``\\sum_k p[k] z^{k-1}`` for the coefficients `p[1]`, `p[2]`, ...;
+that is, the coefficients are given in ascending order by power of `z`. This function
+generates efficient code using a Goertzel-like algorithm specialized for complex arguments
+,*without* loops unrolled at compile time. Use this method when the number of coefficients
+is not known at compile time.
+The Goertzel-like algorthim is described in Knuth's Art of Computer Programming,
+Volume 2: Seminumerical Algorithms, Sec. 4.6.4.
+# Example
+```jldoctest
+julia> evalpoly(2 + im, [1, 2, 3])
+14 + 14im
+
+julia> evalpoly(2 + im, 1:3)
+14 + 14im
+```
+"""
+evalpoly(z::Complex, p::AbstractVector) = _evalpoly(z, p)
+
+function _evalpoly(z::Complex, p)
+ length(p) == 1 && return p[1]
+ N = length(p)
+ a = p[end]
+ b = p[end-1]
+
+ x = real(z)
+ y = imag(z)
+ r = 2x
+ s = muladd(x, x, y*y)
+ for i in N-2:-1:1
+ ai = a
+ a = muladd(r, ai, b)
+ b = p[i] - s * ai
+ end
+ ai = a
+ muladd(ai, z, b)
+end
+
 """
  @horner(x, p...)
 
 Evaluate `p[1] + x * (p[2] + x * (....))`, i.e. a polynomial via Horner's rule.
 """
 macro horner(x, p...)
- ex = esc(p[end])
- for i = length(p)-1:-1:1
- ex = :(muladd(t, $ex, $(esc(p[i]))))
- end
- ex = quote local r = $ex end # structure this to add exactly one line number node for the macro
- return Expr(:block, :(local t = $(esc(x))), ex, :r)
+ xesc, pesc = esc(x), esc.(p)
+ :(invoke(evalpoly, Tuple{Any, Tuple}, $xesc, ($(pesc...),)))
 end
 
 # Evaluate p[1] + z*p[2] + z^2*p[3] + ... + z^(n-1)*p[n]. This uses
@@ -121,28 +246,8 @@ julia> @evalpoly(2, 1, 1, 1)
 ```
 """
 macro evalpoly(z, p...)
- a = :($(esc(p[end])))
- b = :($(esc(p[end-1])))
- as = []
- for i = length(p)-2:-1:1
- ai = Symbol("a", i)
- push!(as, :($ai = $a))
- a = :(muladd(r, $ai, $b))
- b = :($(esc(p[i])) - s * $ai) # see issue #15985 on fused mul-subtract
- end
- ai = :a0
- push!(as, :($ai = $a))
- C = Expr(:block,
- :(x = real(tt)),
- :(y = imag(tt)),
- :(r = x + x),
- :(s = muladd(x, x, y*y)),
- as...,
- :(muladd($ai, tt, $b)))
- R = Expr(:macrocall, Symbol("@horner"), (), :tt, map(esc, p)...)
- :(let tt = $(esc(z))
- isa(tt, Complex) ? $C : $R
- end)
+ zesc, pesc = esc(z), esc.(p)
+ :(evalpoly($zesc, ($(pesc...),)))
 end
 
 """

test/math.jl

@@ -487,17 +487,30 @@ end
 
 @testset "evalpoly" begin
  @test @evalpoly(2,3,4,5,6) == 3+2*(4+2*(5+2*6)) == @evalpoly(2+0im,3,4,5,6)
- @test let evalcounts=0
- @evalpoly(begin
- evalcounts += 1
- 4
- end, 1,2,3,4,5)
- evalcounts
- end == 1
  a0 = 1
  a1 = 2
  c = 3
  @test @evalpoly(c, a0, a1) == 7
+ @test @evalpoly(1, 2) == 2
+end
+
+@testset "evalpoly real" begin
+ for x in -1.0:2.0, p1 in -3.0:3.0, p2 in -3.0:3.0, p3 in -3.0:3.0
+ evpm = @evalpoly(x, p1, p2, p3)
+ @test evalpoly(x, (p1, p2, p3)) == evpm
+ @test evalpoly(x, [p1, p2, p3]) == evpm
+ end
+end
+
+@testset "evalpoly complex" begin
+ for x in -1.0:2.0, y in -1.0:2.0, p1 in -3.0:3.0, p2 in -3.0:3.0, p3 in -3.0:3.0
+ z = x + im * y
+ evpm = @evalpoly(z, p1, p2, p3)
+ @test evalpoly(z, (p1, p2, p3)) == evpm
+ @test evalpoly(z, [p1, p2, p3]) == evpm
+ end
+ @test evalpoly(1+im, (2,)) == 2
+ @test evalpoly(1+im, [2,]) == 2
 end
 
 @testset "cis" begin