Skip to content

Commit ad58cd1

Browse files
authored
closes #338 (#339)
* add RationalFunction(p) constructor * change broadcast behaviour * close issue #338 * adjust evaluation * version bump * fix printing * bug in plot recipe
1 parent 168f7e8 commit ad58cd1

File tree

5 files changed

+48
-19
lines changed

5 files changed

+48
-19
lines changed

Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@ name = "Polynomials"
22
uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
33
license = "MIT"
44
author = "JuliaMath"
5-
version = "2.0.10"
5+
version = "2.0.11"
66

77
[deps]
88
Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"

src/rational-functions/common.jl

Lines changed: 28 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -19,11 +19,11 @@ abstract type AbstractRationalFunction{T,X,P} end
1919

2020

2121
function Base.show(io::IO, pq::AbstractRationalFunction)
22-
p,q = pq
22+
p,q = pqs(pq)
2323
print(io,"(")
24-
print(io, p)
24+
printpoly(io, p)
2525
print(io, ") // (")
26-
print(io, q)
26+
printpoly(io, q)
2727
print(io, ")")
2828
end
2929

@@ -40,7 +40,8 @@ function Base.convert(::Type{PQ}, pq′::PQ′) where {T,X,P,PQ <: AbstractRatio
4040
T′,X′,P′,PQ′<:AbstractRationalFunction{T′,X′,P′} }
4141
!isconstant(pq′) && assert_same_variable(X,X′)
4242
p′,q′=pqs(pq′)
43-
p,q = convert(P, p′), convert(P, q′)
43+
𝑷 = isconstant(pq′) ? P : promote_type(P, P′)
44+
p,q = convert(𝑷, p′), convert(𝑷, q′)
4445
rational_function(PQ, p, q)
4546
end
4647

@@ -81,7 +82,11 @@ function Base.promote_rule(::Type{PQ}, ::Type{PQ′}) where {T,X,P,PQ <: Abstrac
8182
𝑷𝑸{𝑻,X,𝑷{𝑻,X}}
8283
end
8384
Base.promote_rule(::Type{PQ}, ::Type{P}) where {PQ <: AbstractRationalFunction, P<:AbstractPolynomial} = PQ
84-
Base.promote_rule(::Type{PQ}, ::Type{P}) where {PQ <: AbstractRationalFunction, P<:Number} = PQ
85+
function Base.promote_rule(::Type{PQ}, ::Type{S}) where {T,X, P<:AbstractPolynomial{T,X}, PQ <: AbstractRationalFunction{T,X,P}, S<:Number}
86+
R = promote_type(S,T)
87+
P′ = constructorof(P){R,X}
88+
constructorof(PQ){R,X,P′}
89+
end
8590

8691

8792
## Look like rational numbers
@@ -101,9 +106,13 @@ function Base.://(p::PQ, q::Union{Number,AbstractPolynomial}) where {PQ <: Abstr
101106
p0, p1 = p
102107
rational_function(PQ, p0, p1*q)
103108
end
109+
110+
Base.://(p::AbstractPolynomial,q::Number) = p // (q*one(p))
111+
Base.://(p::Number, q::AbstractPolynomial) = (p*one(q)) // q
112+
104113

105114
function Base.copy(pq::PQ) where {PQ <: AbstractRationalFunction}
106-
p,q = pq
115+
p,q = pqs(pq)
107116
rational_function(PQ, p, q)
108117
end
109118

@@ -122,7 +131,7 @@ function Base.iterate(pq::AbstractRationalFunction, state=nothing)
122131
nothing
123132
end
124133
Base.collect(pq::AbstractRationalFunction{T,X,P}) where {T,X,P} = collect(P, pq)
125-
134+
Base.broadcastable(pq::AbstractRationalFunction) = Ref(pq)
126135

127136
Base.eltype(pq::Type{<:AbstractRationalFunction{T,X,P}}) where {T,X,P} = P
128137
Base.eltype(pq::Type{<:AbstractRationalFunction{T,X}}) where {T,X} = Polynomial{T,X}
@@ -134,7 +143,6 @@ Base.eltype(pq::Type{<:AbstractRationalFunction}) = Polynomial{Float64,:x}
134143
pqs(pq)
135144
136145
Return `(p,q)`, where `pq=p/q`, as polynomials.
137-
Alternative to simply `p,q=pq` in case `pq` is not stored as two polynomials.
138146
"""
139147
pqs(pq::AbstractRationalFunction) = (numerator(pq), denominator(pq))
140148

@@ -171,7 +179,10 @@ function indeterminate(::Type{PQ}, var=:x) where {PQ<:AbstractRationalFunction}
171179
X
172180
end
173181

174-
isconstant(pq::AbstractRationalFunction; kwargs...) = all(isconstant.(lowest_terms(pq;kwargs...)))
182+
function isconstant(pq::AbstractRationalFunction; kwargs...)
183+
p,q = pqs(lowest_terms(pq, kwargs...))
184+
isconstant(p) && isconstant(q)
185+
end
175186
isconstant(::Number) = true
176187

177188
function constantterm(pq::AbstractRationalFunction; kwargs...)
@@ -202,14 +213,16 @@ end
202213
# use degree as largest degree of p,q after reduction
203214
function degree(pq::AbstractRationalFunction)
204215
pq′ = lowest_terms(pq)
205-
maximum(degree.(pq′))
216+
maximum(degree.(pqs(pq′)))
206217
end
207218

208219
# Evaluation
209-
function eval_rationalfunction(x, pq::AbstractRationalFunction)
210-
md = minimum(degree.(pq))
220+
function eval_rationalfunction(x, pq::AbstractRationalFunction{T}) where {T}
211221
num, den = pqs(pq)
222+
dn, dd = degree(num), degree(den)
223+
md = min(dn, dd)
212224
result = num(x)/den(x)
225+
md < 0 && return result
213226
while md >= 0
214227
!isnan(result) && return result
215228
num,den = derivative(num), derivative(den)
@@ -240,8 +253,8 @@ end
240253

241254

242255
function Base.isapprox(pq₁::PQ₁, pq₂::PQ₂,
243-
rtol::Real = sqrt(eps(real(promote_type(T,S)))),
244-
atol::Real = zero(real(promote_type(T,S)))) where {T,X,P,PQ₁<:AbstractRationalFunction{T,X,P},
256+
rtol::Real = sqrt(eps(float(real(promote_type(T,S))))),
257+
atol::Real = zero(float(real(promote_type(T,S))))) where {T,X,P,PQ₁<:AbstractRationalFunction{T,X,P},
245258
S,Y,Q,PQ₂<:AbstractRationalFunction{S,Y,Q}}
246259

247260
p₁,q₁ = pqs(pq₁)
@@ -253,7 +266,7 @@ end
253266

254267
# Arithmetic
255268
function Base.:-(pq::PQ) where {PQ <: AbstractRationalFunction}
256-
p, q = copy.(pq)
269+
p, q = copy(pq)
257270
rational_function(PQ, -p, q)
258271
end
259272

src/rational-functions/plot-recipes.jl

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -55,7 +55,8 @@ function rational_function_trim(pq, a, b, xlims, ylims)
5555
n = 601
5656
xs = range(a,stop=b, length=n)
5757
ys = pq.(xs)
58-
M = max(5, 3*maximum(abs, pq.(cps)), 1.25*maximum(abs, pq.((a,b))))
58+
Mcps = isempty(cps) ? 5 : 3*maximum(abs, pq.(cps))
59+
M = max(5, Mcps, 1.25*maximum(abs, pq.((a,b))))
5960

6061
lo = ylims[1] == nothing ? -M : ylims[1]
6162
hi = ylims[2] == nothing ? M : ylims[2]

src/rational-functions/rational-function.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -55,6 +55,7 @@ end
5555

5656
RationalFunction(p,q) = RationalFunction(promote(p,q)...)
5757
RationalFunction(p::ImmutablePolynomial,q::ImmutablePolynomial) = throw(ArgumentError("Sorry, immutable #polynomials are not a valid polynomial type for RationalFunction"))
58+
RationalFunction(p::AbstractPolynomial) = RationalFunction(p,one(p))
5859

5960
# evaluation
6061
(pq::RationalFunction)(x) = eval_rationalfunction(x, pq)

test/rational-functions.jl

Lines changed: 16 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@ using LinearAlgebra
1212
@test p // q isa RationalFunction
1313
@test p // r isa RationalFunction
1414
@test_throws ArgumentError r // s
15+
@test RationalFunction(p) == p // one(p)
1516

1617
# We expect p::P // q::P (same type polynomial).
1718
# As Immutable Polynomials have N as type parameter, we disallow
@@ -26,19 +27,25 @@ using LinearAlgebra
2627
@test eltype(p//q) == typeof(pp)
2728
u = gcd(p//q...)
2829
@test u/u[end] fromroots(Polynomial, [2,3])
29-
@test degree.(p // q) == [degree(p), degree(q)]
30+
@test degree.(p//q) == degree(p//q) # no broadcast over rational functions
3031

3132
# evaluation
3233
pq = p//q
3334
@test pq(2.5) p(2.5) / q(2.5)
3435
@test pq(2) fromroots([1,3])(2) / fromroots([3,4])(2)
36+
@test zero(pq)(10) == 0
37+
@test one(pq)(10) == 1
38+
3539

3640
# arithmetic
3741
rs = r // (r-1)
3842
x = 1.2
3943
@test (pq + rs)(x) (pq(x) + rs(x))
4044
@test (pq * rs)(x) (pq(x) * rs(x))
4145
@test (-pq)(x) -p(x)/q(x)
46+
@test pq .* (pq, pq) == (pq*pq, pq*pq)
47+
@test pq .* [pq pq] == [pq*pq pq*pq]
48+
4249

4350
# derivative
4451
pq = p // one(p)
@@ -56,9 +63,15 @@ using LinearAlgebra
5663

5764
# lowest terms
5865
pq = p // q
59-
pp, qq = lowest_terms(pq)
66+
pp, qq = Polynomials.pqs(lowest_terms(pq))
6067
@test all(abs.(pp.(roots(qq))) .> 1/2)
6168

69+
#
70+
@test (2p)//(2one(p)) p//one(p)
71+
@test p//one(p) p//one(p) + eps()
72+
x = variable(p)
73+
@test (x*p)//x p // one(p)
74+
6275

6376
end
6477

@@ -102,6 +115,7 @@ end
102115

103116
## T, Polynomial{T} promotes
104117
@test eltype([1, p, pp]) == PP
118+
@test eltype(eltype(eltype([im, p, pp]))) == Complex{Int}
105119

106120
## test mixed types promote polynomial type
107121
@test eltype([pp rr p r]) == PP

0 commit comments

Comments
 (0)