Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

export @evalpoly macro (née @chorner) #7146

Merged
merged 2 commits into from
Jun 9, 2014
Merged

export @evalpoly macro (née @chorner) #7146

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
3d36d38
export @evalpoly macro (née @chorner)
stevengj Jun 5, 2014
9d9176a
fix doc, add NEWS for @evalpoly
stevengj Jun 5, 2014
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Next Next commit
export @evalpoly macro (née @chorner)
  • Loading branch information
@stevengj
stevengj committed Jun 5, 2014
commit 3d36d38d3d0584d85467af9b3233d4abb29cac07
1 change: 1 addition & 0 deletions base/exports.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -268,6 +268,7 @@ export
At_rdiv_Bt,

# scalar math
@evalpoly,
abs,
abs2,
acos,
Expand Down
23 changes: 12 additions & 11 deletions base/math.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 = {}
Expand All @@ -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)
Expand Down
4 changes: 2 additions & 2 deletions base/special/gamma.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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}))
Expand All @@ -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
Expand Down
8 changes: 8 additions & 0 deletions doc/stdlib/base.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
------------

Expand Down