Skip to content

Commit 0e7a1f8

Browse files
committed
add evalpoly code for complex numbers; address copy issue with poly + scalar
1 parent 07da326 commit 0e7a1f8

File tree

1 file changed

+26
-2
lines changed

1 file changed

+26
-2
lines changed

src/polynomials/Polynomial.jl

Lines changed: 26 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -68,11 +68,34 @@ function (p::Polynomial{T})(x::S) where {T,S}
6868
length(p) == 0 && return zero(T) * oS
6969
b = p[end] * oS
7070
@inbounds for i in (lastindex(p) - 1):-1:0
71-
b = p[i]*oS .+ x * b
71+
b = p[i]*oS .+ x * b # not muladd(x,b,p[i]), unless we want to add methods for matrices, ...
7272
end
7373
return b
7474
end
7575

76+
## From base/math.jl from Julia 1.4
77+
function (p::Polynomial{T})(z::S) where {T,S <: Complex}
78+
d = degree(p)
79+
d == -1 && zero(z)
80+
d == 0 && return p[0]
81+
N = d + 1
82+
a = p[end]
83+
b = p[end-1]
84+
85+
x = real(z)
86+
y = imag(z)
87+
r = 2x
88+
s = muladd(x, x, y*y)
89+
for i in d-2:-1:0
90+
ai = a
91+
a = muladd(r, ai, b)
92+
b = p[i] - s * ai
93+
end
94+
ai = a
95+
muladd(ai, z, b)
96+
end
97+
98+
7699
function fromroots(P::Type{<:Polynomial}, r::AbstractVector{T}; var::SymbolLike = :x) where {T <: Number}
77100
n = length(r)
78101
c = zeros(T, n + 1)
@@ -173,7 +196,8 @@ end
173196

174197
function Base.:+(p::Polynomial{T}, c::S) where {T,S<:Number}
175198
U = promote_type(T, S)
176-
p2 = U == S ? copy(p) : convert(Polynomial{U}, p)
199+
q = copy(p)
200+
p2 = U == S ? q : convert(Polynomial{U}, q)
177201
p2[0] += c
178202
return p2
179203
end

0 commit comments

Comments
 (0)