Skip to content

Commit

Permalink
Fixes in _pow for Interval coeffs, and tests
Browse files Browse the repository at this point in the history
  • Loading branch information
lbenet committed Aug 24, 2024
1 parent 2b5f19e commit 8acc4ff
Show file tree
Hide file tree
Showing 2 changed files with 48 additions and 52 deletions.
73 changes: 28 additions & 45 deletions ext/TaylorSeriesIAExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,67 +8,50 @@ import TaylorSeries: evaluate, _evaluate, normalize_taylor, square

isdefined(Base, :get_extension) ? (using IntervalArithmetic) : (using ..IntervalArithmetic)

# Method used for Taylor1{Interval{T}}^n

# TS._pow
for T in (:Taylor1, :TaylorN)
@eval begin
function ^(a::$T{Interval{S}}, n::Integer) where {S<:Real}
n == 0 && return one(a)
n == 1 && return copy(a)
n == 2 && return square(a)
n < 0 && return a^float(n)
return power_by_squaring(a, n)
# Use power_by_squaring!
function TS._pow(a::$T{Interval{S}}, n::Integer) where {S<:Real}
aux = one(constant_term(a))^n
c = $T( zero(aux), a.order)
caux = zero(c)
TS.power_by_squaring!(c, a, caux, n)
return c
end
^(a::$T{Interval{S}}, r::Rational) where {S<:Real} = a^float(r)
end
end

function ^(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real}
function TS._pow(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real}
isinteger(r) && r >= 0 && return a^Int(r)
a0 = constant_term(a) Interval(zero(T), T(Inf))
aux = one(a0)^r

iszero(r) && return Taylor1(aux, a.order)
aa = one(aux) * a
aa[0] = one(aux) * a0
r == 1 && return aa
r == 2 && return square(aa)
r == 1/2 && return sqrt(aa)

a[0] = aux * a0
#
l0 = findfirst(a)
lnull = trunc(Int, r*l0 )
if (a.order-lnull < 0) || (lnull > a.order)
return Taylor1( zero(aux), a.order)
end
c_order = l0 == 0 ? a.order : min(a.order, trunc(Int,r*a.order))
(lnull > a.order) && return Taylor1( zero(aux), a.order)
c_order = l0 == 0 ? a.order : min(a.order, trunc(Int, r*a.order))
c = Taylor1(zero(aux), c_order)
for k = 0:c_order
TS.pow!(c, aa, c, r, k)
aux0 = deepcopy(c)
for k in eachindex(c)
TS.pow!(c, a, aux0, r, k)
end

return c
end
function ^(a::TaylorN{Interval{T}}, r::S) where {T<:Real, S<:Real}

function TS._pow(a::TaylorN{Interval{T}}, r::S) where {T<:Real, S<:Real}
isinteger(r) && r >= 0 && return a^Int(r)
a0 = constant_term(a) Interval(zero(T), T(Inf))
a0r = a0^r
aux = one(a0r)

iszero(r) && return TaylorN(aux, a.order)
aa = aux * a
aa[0] = aux * a0
r == 1 && return aa
r == 2 && return square(aa)
r == 1/2 && return sqrt(aa)
isinteger(r) && return aa^round(Int,r)

# @assert !iszero(a0)
iszero(a0) && throw(DomainError(a,
"""The 0-th order TaylorN coefficient must be non-zero
in order to expand `^` around 0."""))

c = TaylorN( a0r, a.order)
for ord in 1:a.order
TS.pow!(c, aa, c, r, ord)
aux = one(a0^r)
a[0] = aux * a0
#
c = TaylorN( zero(aux), a.order)
aux0 = zero(c)
for ord in eachindex(c)
TS.pow!(c, a, aux0, r, ord)
end

return c
end

Expand Down
27 changes: 20 additions & 7 deletions test/intervals.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ using Test
@testset "Tests Taylor1 and TaylorN expansions over Intervals" begin
a = 1..2
b = -1 .. 1
p3(x, a) = x^3 + 3*a*x^2 + 3*a^2*x + a^3
p4(x, a) = x^4 + 4*a*x^3 + 6*a^2*x^2 + 4*a^3*x + a^4
p5(x, a) = x^5 + 5*a*x^4 + 10*a^2*x^3 + 10*a^3*x^2 + 5*a^4*x + a^5

Expand All @@ -24,22 +25,34 @@ using Test
@test normalize_taylor(ti) == ti
@test normalize_taylor(x) == x

@test p4(ti,-a) == (ti-a)^4
@test p5(ti,-a) == (ti-a)^5
@test p4(ti,-b) == (ti-b)^4
@test p3(ti,-a) == (ti-a)^3 == (ti-a)^3.0
@test p4(ti,-a) == (ti-a)^4 == (ti-a)^4.0
@test p5(ti,-a) == (ti-a)^5 == (ti-a)^5.0
@test (ti-b)^3 == (ti-b)^3.0
@test all((p3(ti,-b)).coeffs .⊆ ((ti-b)^3).coeffs)
@test p4(ti,-b) == (ti-b)^4 == (ti-b)^4.0
@test (ti-b)^5 == (ti-b)^5.0
@test all((p5(ti,-b)).coeffs .⊆ ((ti-b)^5).coeffs)


@test p4(x,-y) == (x-y)^4
@test p5(x,-y) == (x-y)^5
@test p4(x,-a) == (x-a)^4
@test p5(x,-a) == (x-a)^5
@test p3(x,-y) == (x-y)^3 == (x-y)^3.0
@test p4(x,-y) == (x-y)^4 == (x-y)^4.0
@test p5(x,-y) == (x-y)^5 == (x-y)^5.0
@test p3(x,-a) == (x-a)^3 == (x-a)^3.0
@test p4(x,-a) == (x-a)^4 == (x-a)^4.0
@test p5(x,-a) == (x-a)^5 == (x-a)^5.0
@test (x-b)^3 == (x-b)^3.0
for ind in eachindex(p3(x,-b))
@test all((p3(x,-b)[ind]).coeffs .⊆ (((x-b)^3)[ind]).coeffs)
end
@test p4(x,-b) == (x-b)^4
@test (x-b)^5 == (x-b)^5.0
for ind in eachindex(p5(x,-b))
@test all((p5(x,-b)[ind]).coeffs .⊆ (((x-b)^5)[ind]).coeffs)
end

# Tests `evaluate`
@test evaluate(p3(x,y), IntervalBox(a,-b)) == p3(a, -b)
@test evaluate(p4(x,y), IntervalBox(a,-b)) == p4(a, -b)
@test (p5(x,y))(IntervalBox(a,b)) == p5(a, b)
@test (a-b)^4 ((x-y)^4)(a × b)
Expand Down

0 comments on commit 8acc4ff

Please sign in to comment.