From e7b85e6a0bea4e0a56f7dedeb22458e8454f262a Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Sun, 3 Dec 2023 13:02:14 -0600 Subject: [PATCH 01/13] Adapt changes in IA --- ext/TaylorSeriesIAExt.jl | 101 ++++++++++++++++------------ test/intervals.jl | 141 ++++++++++++++++++++------------------- 2 files changed, 131 insertions(+), 111 deletions(-) diff --git a/ext/TaylorSeriesIAExt.jl b/ext/TaylorSeriesIAExt.jl index d526bd16..80ca2a5c 100644 --- a/ext/TaylorSeriesIAExt.jl +++ b/ext/TaylorSeriesIAExt.jl @@ -10,20 +10,17 @@ isdefined(Base, :get_extension) ? (using IntervalArithmetic) : (using ..Interval # Method used for Taylor1{Interval{T}}^n 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) - end - ^(a::$T{Interval{S}}, r::Rational) where {S<:Real} = a^float(r) + @eval function ^(a::$T{Interval{T}}, n::S) where {T<:Real, S<:Integer} + 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) end end function ^(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real} - a0 = constant_term(a) ∩ Interval(zero(T), T(Inf)) + a0 = constant_term(a) ∩ interval(zero(T), T(Inf)) aux = one(a0)^r iszero(r) && return Taylor1(aux, a.order) @@ -47,7 +44,7 @@ function ^(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real} return c end function ^(a::TaylorN{Interval{T}}, r::S) where {T<:Real, S<:Real} - a0 = constant_term(a) ∩ Interval(zero(T), T(Inf)) + a0 = constant_term(a) ∩ interval(zero(T), T(Inf)) a0r = a0^r aux = one(a0r) @@ -77,7 +74,7 @@ for T in (:Taylor1, :TaylorN) @eval function log(a::$T{Interval{S}}) where {S<:Real} iszero(constant_term(a)) && throw(DomainError(a, """The 0-th order coefficient must be non-zero in order to expand `log` around 0.""")) - a0 = constant_term(a) ∩ Interval(S(0.0), S(Inf)) + a0 = constant_term(a) ∩ interval(zero(S), S(Inf)) order = a.order aux = log(a0) aa = one(aux) * a @@ -90,7 +87,7 @@ for T in (:Taylor1, :TaylorN) end @eval function asin(a::$T{Interval{S}}) where {S<:Real} - a0 = constant_term(a) ∩ Interval(-one(S), one(S)) + a0 = constant_term(a) ∩ interval(-one(S), one(S)) a0^2 == one(a0) && throw(DomainError(a, """Series expansion of asin(x) diverges at x = ±1.""")) @@ -107,7 +104,7 @@ for T in (:Taylor1, :TaylorN) end @eval function acos(a::$T{Interval{S}}) where {S<:Real} - a0 = constant_term(a) ∩ Interval(-one(S), one(S)) + a0 = constant_term(a) ∩ interval(-one(S), one(S)) a0^2 == one(a0) && throw(DomainError(a, """Series expansion of asin(x) diverges at x = ±1.""")) @@ -124,7 +121,7 @@ for T in (:Taylor1, :TaylorN) end @eval function acosh(a::$T{Interval{S}}) where {S<:Real} - a0 = constant_term(a) ∩ Interval(one(S), S(Inf)) + a0 = constant_term(a) ∩ interval(one(S), S(Inf)) a0^2 == one(a0) && throw(DomainError(a, """Series expansion of acosh(x) diverges at x = ±1.""")) @@ -142,7 +139,7 @@ for T in (:Taylor1, :TaylorN) @eval function atanh(a::$T{Interval{S}}) where {S<:Real} order = a.order - a0 = constant_term(a) ∩ Interval(-one(S), one(S)) + a0 = constant_term(a) ∩ interval(-one(S), one(S)) aux = atanh(a0) aa = one(aux) * a aa[0] = one(aux) * a0 @@ -179,32 +176,53 @@ function evaluate(a::Taylor1, dx::Interval{S}) where {S<:Real} return sum_even + sum_odd*dx end -function evaluate(a::TaylorN, dx::IntervalBox{N,T}) where {T<:Real,N} - @assert N == get_numvars() - a_length = length(a) - suma = zero(constant_term(a)) + Interval{T}(0, 0) +function evaluate(a::TaylorN, dx::AbstractVector{Interval{T}}) where {T<:Real} + @assert length(dx) == get_numvars() + suma = zero(constant_term(a)) + interval(zero(T)) @inbounds for homPol in reverse(eachindex(a)) - suma += evaluate(a[homPol], dx) + suma += evaluate(a.coeffs[homPol+1], dx) end return suma end -function evaluate(a::HomogeneousPolynomial, dx::IntervalBox{N,T}) where {T<:Real,N} - @assert N == get_numvars() - dx == IntervalBox(-1..1, Val(N)) && return _evaluate(a, dx, Val(true)) - dx == IntervalBox( 0..1, Val(N)) && return _evaluate(a, dx, Val(false)) +function evaluate(a::Taylor1{TaylorN{T}}, dx::Interval{S}) where {T<:Real, S<:Real} + order = a.order + uno = one(dx) + dx2 = dx^2 + if iseven(order) + kend = order-2 + @inbounds sum_even = a[end]*uno + @inbounds sum_odd = a[end-1]*zero(dx) + else + kend = order-3 + @inbounds sum_odd = a[end]*uno + @inbounds sum_even = a[end-1]*uno + end + @inbounds for k in kend:-2:0 + sum_odd = sum_odd*dx2 + a[k+1] + sum_even = sum_even*dx2 + a[k] + end + return sum_even + sum_odd*dx +end + + +function evaluate(a::HomogeneousPolynomial, dx::AbstractVector{Interval{T}}) where {T<:Real} + @assert length(dx) == get_numvars() + all(dx .== interval(-one(T), one(T))) && return _evaluate(a, dx, Val(true)) + all(dx .== interval(zero(T), one(T))) && return _evaluate(a, dx, Val(false)) return evaluate(a, dx...) end -function _evaluate(a::HomogeneousPolynomial, dx::IntervalBox{N,T}, ::Val{true} ) where {T<:Real,N} - a.order == 0 && return a[1] + Interval{T}(0, 0) +function _evaluate(a::HomogeneousPolynomial, dx::AbstractVector{Interval{T}}, + ::Val{true} ) where {T<:Real} + a.order == 0 && return a[1] + interval(zero(T)) ct = TS.coeff_table[a.order+1] - @inbounds suma = a[1]*Interval{T}(0,0) + @inbounds suma = a[1]*interval(zero(T)) - Ieven = Interval{T}(0,1) + Ieven = interval(zero(T), one(T)) for (i, a_coeff) in enumerate(a.coeffs) iszero(a_coeff) && continue if isodd(sum(ct[i])) @@ -221,8 +239,9 @@ function _evaluate(a::HomogeneousPolynomial, dx::IntervalBox{N,T}, ::Val{true} ) return suma end -function _evaluate(a::HomogeneousPolynomial, dx::IntervalBox{N,T}, ::Val{false} ) where {T<:Real,N} - a.order == 0 && return a[1] + Interval{T}(0, 0) +function _evaluate(a::HomogeneousPolynomial, dx::AbstractVector{Interval{T}}, + ::Val{false} ) where {T<:Real} + a.order == 0 && return a[1] + interval(zero(T)) @inbounds suma = zero(a[1])*dx[1] @inbounds for homPol in a.coeffs @@ -243,14 +262,14 @@ normalize_taylor(a::Taylor1, I::Interval{T}, symI::Bool=true) where {T<:Real} = _normalize(a, I, Val(symI)) """ - normalize_taylor(a::TaylorN, I::IntervalBox, symI::Bool=true) + normalize_taylor(a::TaylorN, I::AbstractVector{Interval{T}}, symI::Bool=true) -Normalize `a::TaylorN` such that the intervals in `I::IntervalBox` +Normalize `a::TaylorN` such that the intervals in `I::AbstractVector{Interval{T}}` are mapped by an affine transformation to the intervals `-1..1` (`symI=true`) or to `0..1` (`symI=false`). """ -normalize_taylor(a::TaylorN, I::IntervalBox{N,T}, symI::Bool=true) where {T<:Real,N} = - _normalize(a, I, Val(symI)) +normalize_taylor(a::TaylorN, I::AbstractVector{Interval{T}}, + symI::Bool=true) where {T<:Real} = _normalize(a, I, Val(symI)) # I -> -1..1 function _normalize(a::Taylor1, I::Interval{T}, ::Val{true}) where {T<:Real} @@ -269,20 +288,20 @@ function _normalize(a::Taylor1, I::Interval{T}, ::Val{false}) where {T<:Real} end # I -> IntervalBox(-1..1, Val(N)) -function _normalize(a::TaylorN, I::IntervalBox{N,T}, ::Val{true}) where {T<:Real,N} +function _normalize(a::TaylorN, I::AbstractVector{Interval{T}}, ::Val{true}) where {T<:Real} order = get_order(a) - x = Vector{typeof(a)}(undef, N) - for ind in eachindex(x) + x = Vector{typeof(a)}(undef, length(I)) + @inbounds for ind in eachindex(x) x[ind] = mid(I[ind]) + TaylorN(ind, order=order)*radius(I[ind]) end return a(x) end # I -> IntervalBox(0..1, Val(N)) -function _normalize(a::TaylorN, I::IntervalBox{N,T}, ::Val{false}) where {T<:Real,N} +function _normalize(a::TaylorN, I::AbstractVector{Interval{T}}, ::Val{false}) where {T<:Real} order = get_order(a) - x = Vector{typeof(a)}(undef, N) - for ind in eachindex(x) + x = Vector{typeof(a)}(undef, length(I)) + @inbounds for ind in eachindex(x) x[ind] = inf(I[ind]) + TaylorN(ind, order=order)*diam(I[ind]) end return a(x) diff --git a/test/intervals.jl b/test/intervals.jl index d163acfa..6c69163c 100644 --- a/test/intervals.jl +++ b/test/intervals.jl @@ -7,8 +7,9 @@ using Test # eeuler = Base.MathConstants.e @testset "Tests Taylor1 and TaylorN expansions over Intervals" begin - a = 1..2 - b = -1 .. 1 + a = interval(1, 2) + b = interval(-1, 1) + c = interval(0, 1) 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 @@ -40,65 +41,65 @@ using Test end # Tests `evaluate` - @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) - @test (((x-y)^4)[4])(a × b) == -39 .. 81 - - p4n = normalize_taylor(p4(x,y), a × b, true) - @test (0..16) ⊆ p4n((-1..1)×(-1..1)) - p5n = normalize_taylor(p5(x,y), a × b, true) - @test (-32 .. 32) ⊆ p5n((-1..1)×(-1..1)) - - p4n = normalize_taylor(p4(x,y), a × b, false) - @test (0..16) ⊆ p4n((0..1)×(0..1)) - p5n = normalize_taylor(p5(x,y), a × b, false) - @test (0..32) ⊆ p5n((0..1)×(0..1)) - - @test evaluate(x*y^3, (-1..1)×(-1..1)) == (-1..1) - @test evaluate(x*y^2, (-1..1)×(-1..1)) == (-1..1) - @test evaluate(x^2*y^2, (-1..1)×(-1..1)) == (0..1) - - ii = -1..1 + @test evaluate(p4(x,y), [a, -b]) == p4(a, -b) + @test (p5(x,y))([a, b]) == p5(a, b) + @test (a-b)^4 ⊆ ((x-y)^4)([a, b]) + @test (((x-y)^4)[4])([a, b]) == interval(-39, 81) + + p4n = normalize_taylor(p4(x,y), [a, b], true) + @test interval(0, 16) ⊆ p4n([b, b]) + p5n = normalize_taylor(p5(x,y), [a, b], true) + @test interval(-32, 32) ⊆ p5n([b, b]) + + p4n = normalize_taylor(p4(x,y), [a, b], false) + @test interval(0, 16) ⊆ p4n([c, c]) + p5n = normalize_taylor(p5(x,y), [a, b], false) + @test interval(0, 32) ⊆ p5n([c, c]) + + @test evaluate(x*y^3, [b, b]) == b + @test evaluate(x*y^2, [b, b]) == b + @test evaluate(x^2*y^2, [b, b]) == c + + ii = b t = Taylor1(1) - @test 0..2 ⊆ (1+t)(ii) + @test interval(0, 2) ⊆ (1+t)(ii) t = Taylor1(2) - @test 0..4 ⊆ ((1+t)^2)(ii) + @test interval(0, 4) ⊆ ((1+t)^2)(ii) - ii = 0..6 + ii = interval(0, 6) t = Taylor1(4) f(x) = 0.1 * x^3 - 0.5*x^2 + 1 ft = f(t) f1 = normalize_taylor(ft, ii, true) f2 = normalize_taylor(ft, ii, false) - @test Interval(-23/27, f(6)) ⊆ f(ii) - @test Interval(-23/27, f(6)) ⊆ ft(ii) - @test Interval(-23/27, f(6)) ⊆ f1(-1..1) - @test Interval(-23/27, f(6)) ⊆ f2(0..1) - @test f1(-1..1) ⊆ f(ii) - @test diam(f1(-1..1)) < diam(f2(0..1)) + @test interval(-23/27, f(6)) ⊆ f(ii) + @test interval(-23/27, f(6)) ⊆ ft(ii) + @test interval(-23/27, f(6)) ⊆ f1b + @test interval(-23/27, f(6)) ⊆ f2(c) + @test f1b ⊆ f(ii) + @test diam(f1b) < diam(f2(c)) # An example from Makino's thesis - ii = 0..1 + ii = c t = Taylor1(5) g(x) = 1 - x^4 + x^5 gt = g(t) - g1 = normalize_taylor(gt, 0..1, true) - @test Interval(g(4/5),1) ⊆ g(ii) - @test Interval(g(4/5),1) ⊆ gt(ii) - @test Interval(g(4/5),1) ⊆ g1(-1..1) - @test g1(-1..1) ⊂ g(ii) - @test diam(g1(-1..1)) < diam(gt(ii)) + g1 = normalize_taylor(gt, c, true) + @test interval(g(4/5), 1) ⊆ g(ii) + @test interval(g(4/5), 1) ⊆ gt(ii) + @test interval(g(4/5), 1) ⊆ g1b + @test g1b ⊂ g(ii) + @test diam(g1b) < diam(gt(ii)) # Test display for Taylor1{Complex{Interval{T}}} - vc = [complex(1.5 .. 2, 0..0 ), complex(-2 .. -1, -1 .. 1 ), - complex( -1 .. 1.5, -1 .. 1.5), complex( 0..0, -1 .. 1.5)] + vc = [complex(interval(1.5, 2), interval(0, 0)), complex(interval(-2, -1), interval(-1, 1 )), + complex(interval(-1, 1.5), interval(-1, 1.5)), complex( interval(0,0), interval(-1, 1.5))] displayBigO(false) @test string(Taylor1(vc, 5)) == - " ( [1.5, 2] + [0, 0]im ) - ( [1, 2] + [-1, 1]im ) t + ( [-1, 1.5] + [-1, 1.5]im ) t² + ( [0, 0] + [-1, 1.5]im ) t³ " + " ( Interval(1.5, 2.0) + Interval(0.0, 0.0)im ) - ( Interval(1.0, 2.0) + Interval(-1.0, 1.0)im ) t + ( Interval(-1.0, 1.5) + Interval(-1.0, 1.5)im ) t² + ( Interval(0.0, 0.0) + Interval(-1.0, 1.5)im ) t³ " displayBigO(true) @test string(Taylor1(vc, 5)) == - " ( [1.5, 2] + [0, 0]im ) - ( [1, 2] + [-1, 1]im ) t + ( [-1, 1.5] + [-1, 1.5]im ) t² + ( [0, 0] + [-1, 1.5]im ) t³ + 𝒪(t⁶)" + " ( Interval(1.5, 2.0) + Interval(0.0, 0.0)im ) - ( Interval(1.0, 2.0) + Interval(-1.0, 1.0)im ) t + ( Interval(-1.0, 1.5) + Interval(-1.0, 1.5)im ) t² + ( Interval(0.0, 0.0) + Interval(-1.0, 1.5)im ) t³ + 𝒪(t⁶)" # Iss 351 (inspired by a test in ReachabilityAnalysis) p1 = Taylor1([0 .. 0, (0 .. 0.1) + (0 .. 0.01) * y], 4) @@ -109,38 +110,38 @@ using Test # Tests related to Iss #311 # `sqrt` and `pow` defined on Interval(0,Inf) @test_throws DomainError sqrt(ti) - @test sqrt(Interval(0.0, 1.e-15) + ti) == sqrt(Interval(-1.e-15, 1.e-15) + ti) - aa = sqrt(sqrt(Interval(0.0, 1.e-15) + ti)) - @test aa == sqrt(sqrt(Interval(-1.e-15, 1.e-15) + ti)) - bb = (Interval(0.0, 1.e-15) + ti)^(1/4) - @test bb == (Interval(-1.e-15, 1.e-15) + ti)^(1/4) + @test sqrt(interval(0.0, 1.e-15) + ti) == sqrt(interval(-1.e-15, 1.e-15) + ti) + aa = sqrt(sqrt(interval(0.0, 1.e-15) + ti)) + @test aa == sqrt(sqrt(interval(-1.e-15, 1.e-15) + ti)) + bb = (interval(0.0, 1.e-15) + ti)^(1/4) + @test bb == (interval(-1.e-15, 1.e-15) + ti)^(1/4) @test all(aa.coeffs[2:end] .⊂ bb.coeffs[2:end]) @test_throws DomainError sqrt(x) - @test sqrt(Interval(-1,1)+x) == sqrt(Interval(0,1)+x) - @test (Interval(-1,1)+x)^(1/4) == (Interval(0,1)+x)^(1/4) + @test sqrt(interval(-1,1)+x) == sqrt(interval(0,1)+x) + @test (interval(-1,1)+x)^(1/4) == (interval(0,1)+x)^(1/4) # `log` defined on Interval(0,Inf) @test_throws DomainError log(ti) - @test log(Interval(0.0, 1.e-15) + ti) == log(Interval(-1.e-15, 1.e-15) + ti) + @test log(interval(0.0, 1.e-15) + ti) == log(interval(-1.e-15, 1.e-15) + ti) @test_throws DomainError log(y) - @test log(Interval(0.0, 1.e-15) + y) == log(Interval(-1.e-15, 1.e-15) + y) - # `asin` and `acos` defined on Interval(-1,1) - @test_throws DomainError asin(Interval(1.0 .. 2.0) + ti) - @test asin(Interval(-2.0 .. 0.0) + ti) == asin(Interval(-1,0) + ti) - @test_throws DomainError acos(Interval(1.0 .. 2.0) + ti) - @test acos(Interval(-2.0 .. 0.0) + ti) == acos(Interval(-1,0) + ti) - @test_throws DomainError asin(Interval(1.0 .. 2.0) + x) - @test asin(Interval(-2.0 .. 0.0) + x) == asin(Interval(-1,0) + x) - @test_throws DomainError acos(Interval(1.0 .. 2.0) + x) - @test acos(Interval(-2.0 .. 0.0) + x) == acos(Interval(-1,0) + x) - # acosh defined on Interval(1,Inf) - @test_throws DomainError acosh(Interval(0.0 .. 1.0) + ti) - @test acosh(Interval(0.0 .. 2.0) + ti) == acosh(Interval(1.0 .. 2.0) + ti) - @test_throws DomainError acosh(Interval(0.0 .. 1.0) + x) - @test acosh(Interval(0.0 .. 2.0) + x) == acosh(Interval(1.0 .. 2.0) + x) - # atanh defined on Interval(-1,1) - @test_throws DomainError atanh(Interval(1.0 .. 1.0) + ti) - @test atanh(Interval(-2.0 .. 0.0) + ti) == atanh(Interval(-1.0 .. 0.0) + ti) - @test_throws DomainError atanh(Interval(1.0 .. 1.0) + y) - @test atanh(Interval(-2.0 .. 0.0) + y) == atanh(Interval(-1.0 .. 0.0) + y) + @test log(interval(0.0, 1.e-15) + y) == log(interval(-1.e-15, 1.e-15) + y) + # `asin` and `acos` defined on interval(-1,1) + @test_throws DomainError asin(interval(1.0, 2.0) + ti) + @test asin(interval(-2.0, 0.0) + ti) == asin(interval(-1,0) + ti) + @test_throws DomainError acos(interval(1.0, 2.0) + ti) + @test acos(interval(-2.0, 0.0) + ti) == acos(interval(-1,0) + ti) + @test_throws DomainError asin(interval(1.0, 2.0) + x) + @test asin(interval(-2.0, 0.0) + x) == asin(interval(-1,0) + x) + @test_throws DomainError acos(interval(1.0, 2.0) + x) + @test acos(interval(-2.0, 0.0) + x) == acos(interval(-1,0) + x) + # acosh defined on interval(1,Inf) + @test_throws DomainError acosh(interval(0.0, 1.0) + ti) + @test acosh(interval(0.0, 2.0) + ti) == acosh(interval(1.0, 2.0) + ti) + @test_throws DomainError acosh(interval(0.0, 1.0) + x) + @test acosh(interval(0.0, 2.0) + x) == acosh(interval(1.0, 2.0) + x) + # atanh defined on interval(-1,1) + @test_throws DomainError atanh(interval(1.0, 1.0) + ti) + @test atanh(interval(-2.0, 0.0) + ti) == atanh(interval(-1.0, 0.0) + ti) + @test_throws DomainError atanh(interval(1.0, 1.0) + y) + @test atanh(interval(-2.0, 0.0) + y) == atanh(interval(-1.0, 0.0) + y) end From 25090702e70b165d426ff9030c05d1726e885658 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Mon, 11 Dec 2023 19:40:16 -0600 Subject: [PATCH 02/13] Add _isthinzero as alias of iszero --- src/auxiliary.jl | 65 +++++++++++++++++++++++------------------------- src/printing.jl | 20 +++++++-------- 2 files changed, 41 insertions(+), 44 deletions(-) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 534b1054..9ddd4f67 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -290,41 +290,38 @@ for T in (:HomogeneousPolynomial, :TaylorN) end -# Finds the first non zero entry; extended to Taylor1 -function Base.findfirst(a::Taylor1{T}) where {T<:Number} - first = findfirst(x->!iszero(x), a.coeffs) - isnothing(first) && return -1 - return first-1 -end -# Finds the last non-zero entry; extended to Taylor1 -function Base.findlast(a::Taylor1{T}) where {T<:Number} - last = findlast(x->!iszero(x), a.coeffs) - isnothing(last) && return -1 - return last-1 -end +## _isthinzero +""" + _isthinzero(x) -# Finds the first non zero entry; extended to HomogeneousPolynomial -function Base.findfirst(a::HomogeneousPolynomial{T}) where {T<:Number} - first = findfirst(x->!iszero(x), a.coeffs) - isa(first, Nothing) && return -1 - return first -end -function Base.findfirst(a::TaylorN{T}) where {T<:Number} - first = findfirst(x->!iszero(x), a.coeffs) - isa(first, Nothing) && return -1 - return first-1 -end -# Finds the last non-zero entry; extended to HomogeneousPolynomial -function Base.findlast(a::HomogeneousPolynomial{T}) where {T<:Number} - last = findlast(x->!iszero(x), a.coeffs) - isa(last, Nothing) && return -1 - return last -end -# Finds the last non-zero entry; extended to TaylorN -function Base.findlast(a::TaylorN{T}) where {T<:Number} - last = findlast(x->!iszero(x), a.coeffs) - isa(last, Nothing) && return -1 - return last-1 +Generic wrapper to function `iszero`, which allows using the correct +function for `Interval`s +""" +_isthinzero(x) = iszero(x) + + +## findfirst, findlast +for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) + # Finds the first non zero entry + @eval function Base.findfirst(a::$T{T}) where {T<:Number} + first = findfirst(x->!_isthinzero(x), a.coeffs) + isnothing(first) && return -1 + if $T == HomogeneousPolynomial + return first + else + return first-1 + end + end + # Finds the last non-zero entry + @eval function Base.findlast(a::$T{T}) where {T<:Number} + last = findlast(x->!_isthinzero(x), a.coeffs) + isnothing(last) && return -1 + if $T == HomogeneousPolynomial + return last + else + return last-1 + end + end end diff --git a/src/printing.jl b/src/printing.jl index 0ba7c835..c5f0b861 100644 --- a/src/printing.jl +++ b/src/printing.jl @@ -27,7 +27,7 @@ function pretty_print(a::Taylor1) bigO = bigOnotation[end] ? string("+ 𝒪(", var, superscriptify(a.order+1), ")") : string("") - # iszero(a) && return string(space, z, space, bigO) + TS._isthinzero(a) && return string(space, z, space, bigO) strout::String = space ifirst = true for i in eachindex(a) @@ -50,14 +50,14 @@ function pretty_print(a::Taylor1{T}) where {T<:NumberNotSeries} bigO = bigOnotation[end] ? string("+ 𝒪(", var, superscriptify(a.order+1), ")") : string("") - iszero(a) && return string(space, z, space, bigO) + TS._isthinzero(a) && return string(space, z, space, bigO) strout::String = space ifirst = true for i in eachindex(a) monom::String = i==0 ? string("") : i==1 ? string(" ", var) : string(" ", var, superscriptify(i)) @inbounds c = a[i] - iszero(c) && continue + TS._isthinzero(c) && continue cadena = numbr2str(c, ifirst) strout = string(strout, cadena, monom, space) ifirst = false @@ -73,14 +73,14 @@ function pretty_print(a::Taylor1{T} where {T <: AbstractSeries{S}}) where {S<:Nu bigO = bigOnotation[end] ? string("+ 𝒪(", var, superscriptify(a.order+1), ")") : string("") - iszero(a) && return string(space, z, space, bigO) + TS._isthinzero(a) && return string(space, z, space, bigO) strout::String = space ifirst = true for i in eachindex(a) monom::String = i==0 ? string("") : i==1 ? string(" ", var) : string(" ", var, superscriptify(i)) @inbounds c = a[i] - iszero(c) && continue + TS._isthinzero(c) && continue cadena = numbr2str(c, ifirst) ccad::String = i==0 ? cadena : ifirst ? string("(", cadena, ")") : string(cadena[1:2], "(", cadena[3:end], ")") @@ -94,7 +94,7 @@ end function pretty_print(a::HomogeneousPolynomial{T}) where {T<:Number} z = zero(a[1]) space = string(" ") - iszero(a) && return string(space, z) + TS._isthinzero(a) && return string(space, z) strout::String = homogPol2str(a) strout end @@ -105,12 +105,12 @@ function pretty_print(a::TaylorN{T}) where {T<:Number} bigO::String = bigOnotation[end] ? string(" + 𝒪(‖x‖", superscriptify(a.order+1), ")") : string("") - iszero(a) && return string(space, z, space, bigO) + TS._isthinzero(a) && return string(space, z, space, bigO) strout::String = space ifirst = true for ord in eachindex(a) pol = a[ord] - iszero(pol) && continue + TS._isthinzero(pol) && continue cadena::String = homogPol2str( pol ) strsgn = (ifirst || ord == 0 || cadena[2] == '-') ? string("") : string(" +") @@ -141,7 +141,7 @@ function homogPol2str(a::HomogeneousPolynomial{T}) where {T<:Number} end end @inbounds c = a[pos] - iszero(c) && continue + TS._isthinzero(c) && continue cadena = numbr2str(c, ifirst) strout = string(strout, cadena, monom, space) ifirst = false @@ -170,7 +170,7 @@ function homogPol2str(a::HomogeneousPolynomial{Taylor1{T}}) where {T<:Number} end end @inbounds c = a[pos] - iszero(c) && continue + TS._isthinzero(c) && continue cadena = numbr2str(c, ifirst) ccad::String = (pos==1 || ifirst) ? string("(", cadena, ")") : string(cadena[1:2], "(", cadena[3:end], ")") From c9a4931defec02fc644e73f256a2b2be4873ef42 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Tue, 12 Dec 2023 10:58:47 -0600 Subject: [PATCH 03/13] Updates to newer version of IA in TaylorSeriesIAExt --- ext/TaylorSeriesIAExt.jl | 573 ++++++++++++++++++++++++++++++++++++--- test/intervals.jl | 87 +++--- 2 files changed, 582 insertions(+), 78 deletions(-) diff --git a/ext/TaylorSeriesIAExt.jl b/ext/TaylorSeriesIAExt.jl index 80ca2a5c..ff5ff657 100644 --- a/ext/TaylorSeriesIAExt.jl +++ b/ext/TaylorSeriesIAExt.jl @@ -2,32 +2,67 @@ module TaylorSeriesIAExt using TaylorSeries -import Base: ^, log, asin, acos, acosh, atanh, power_by_squaring +import Base: ^, sqrt, log, asin, acos, acosh, atanh, iszero, == -import TaylorSeries: evaluate, _evaluate, normalize_taylor, square +import TaylorSeries: evaluate, _evaluate, normalize_taylor isdefined(Base, :get_extension) ? (using IntervalArithmetic) : (using ..IntervalArithmetic) -# Method used for Taylor1{Interval{T}}^n -for T in (:Taylor1, :TaylorN) - @eval function ^(a::$T{Interval{T}}, n::S) where {T<:Real, S<:Integer} - 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) +# Some functions require special interval functions (isequal_interval, isthinzero) +for I in (:Interval, :ComplexI) + @eval begin + TS._isthinzero(x::$I{T}) where {T<:Real} = isthinzero(x) + + function ==(a::Taylor1{$I{T}}, b::Taylor1{$I{S}}) where {T<:Real, S<:Real} + if a.order != b.order + a, b = fixorder(a, b) + end + return all(isequal_interval.(a.coeffs, b.coeffs)) + end + function ==(a::HomogeneousPolynomial{$I{T}}, + b::HomogeneousPolynomial{$I{S}}) where {T<:Real, S<:Real} + a.order == b.order && return all(isequal_interval.(a.coeffs, b.coeffs)) + return all(TS._isthinzero.(a.coeffs)) && all(TS._isthinzero.(b.coeffs)) + end + + iszero(a::Taylor1{$I{T}}) where {T<:Real} = all(TS._isthinzero.(a.coeffs)) + iszero(a::HomogeneousPolynomial{$I{T}}) where {T<:Real} = + all(TS._isthinzero.(a.coeffs)) + end + +# Methods related to power, sqr, sqrt, ... + for T in (:Taylor1, :TaylorN) + @eval begin + function ^(a::$T{$I{T}}, n::S) where {T<:Real, S<:Integer} + n == 0 && return one(a) + n == 1 && return copy(a) + n == 2 && return TS.square(a) + n < 0 && return a^float(n) + return Base.power_by_squaring(a, n) + end + + function TS.square(a::$T{$I{T}}) where {T<:Real} + c = $T( constant_term(a)^2, a.order) + for k in 1:a.order + TS.sqr!(c, a, k) + end + return c + end + + ^(a::$T{$I{T}}, r::Rational) where {T<:Real} = a^float(r) + end end end function ^(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real} - a0 = constant_term(a) ∩ interval(zero(T), T(Inf)) + a0 = intersect_interval(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 == 2 && return TS.square(aa) r == 1/2 && return sqrt(aa) l0 = findfirst(a) @@ -43,8 +78,34 @@ function ^(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real} return c end +# function ^(a::Taylor1{Complex{Interval{T}}}, r::S) where {T<:Real, S<:Real} +# a0 = intersect_interval(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 TS.square(aa) +# r == 1/2 && return sqrt(aa) + +# 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)) +# c = Taylor1(zero(aux), c_order) +# for k = 0:c_order +# TS.pow!(c, aa, r, k) +# end + +# return c +# end + function ^(a::TaylorN{Interval{T}}, r::S) where {T<:Real, S<:Real} - a0 = constant_term(a) ∩ interval(zero(T), T(Inf)) + isinteger(r) && return a^round(Int,r) + a0 = intersect_interval(constant_term(a), interval(zero(T), T(Inf))) a0r = a0^r aux = one(a0r) @@ -52,12 +113,12 @@ function ^(a::TaylorN{Interval{T}}, r::S) where {T<:Real, S<:Real} aa = aux * a aa[0] = aux * a0 r == 1 && return aa - r == 2 && return square(aa) + r == 2 && return TS.square(aa) r == 1/2 && return sqrt(aa) - isinteger(r) && return aa^round(Int,r) + # isinteger(r) && return aa^round(Int,r) # @assert !iszero(a0) - iszero(a0) && throw(DomainError(a, + TS._isthinzero(a0) && throw(DomainError(a, """The 0-th order TaylorN coefficient must be non-zero in order to expand `^` around 0.""")) @@ -68,13 +129,282 @@ function ^(a::TaylorN{Interval{T}}, r::S) where {T<:Real, S<:Real} return c end +# function ^(a::TaylorN{Complex{Interval{T}}}, r::S) where {T<:Real, S<:Real} +# isinteger(r) && return a^round(Int,r) +# a0 = intersect_interval(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 TS.square(aa) +# r == 1/2 && return sqrt(aa) +# # isinteger(r) && return aa^round(Int,r) + +# # @assert !iszero(a0) +# TS._isthinzero(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, r, ord) +# end + +# return c +# end + + +# sqr! +for T = (:Taylor1, :TaylorN) + @eval begin + @inline function TS.sqr!(c::$T{Interval{T}}, a::$T{Interval{T}}, + k::Int) where {T<:Real} + if k == 0 + @inbounds c[0] = constant_term(a)^2 + return nothing + end + + kodd = k%2 + kend = div(k - 2 + kodd, 2) + @inbounds for i = 0:kend + if $T == Taylor1 + c[k] += a[i] * a[k-i] + else + TS.mul!(c[k], a[i], a[k-i]) + end + end + @inbounds c[k] = interval(2) * c[k] + kodd == 1 && return nothing + + if $T == Taylor1 + @inbounds c[k] += a[div(k,2)]^2 + else + TS.sqr!(c[k], a[div(k,2)]) + end + + return nothing + end + # @inline function TS.sqr!(c::$T{Complex{Interval{T}}}, + # a::$T{Complex{Interval{T}}}, k::Int) where {T<:Real} + # if k == 0 + # @inbounds c[0] = constant_term(a)^2 + # return nothing + # end + + # kodd = k%2 + # kend = div(k - 2 + kodd, 2) + # @inbounds for i = 0:kend + # if $T == Taylor1 + # c[k] += a[i] * a[k-i] + # else + # TS.mul!(c[k], a[i], a[k-i]) + # end + # end + # @inbounds c[k] = interval(2) * c[k] + # kodd == 1 && return nothing + + # if $T == Taylor1 + # @inbounds c[k] += a[div(k,2)]^2 + # else + # TS.sqr!(c[k], a[div(k,2)]) + # end + + # return nothing + # end + end +end +@inline function TS.sqr!(c::HomogeneousPolynomial{Interval{T}}, + a::HomogeneousPolynomial{Interval{T}}) where {T<:Real} + iszero(a) && return nothing + + @inbounds num_coeffs_a = TS.size_table[a.order+1] + + @inbounds posTb = TS.pos_table[c.order+1] + @inbounds idxTb = TS.index_table[a.order+1] + + @inbounds for na = 1:num_coeffs_a + ca = a[na] + TS._isthinzero(ca) && continue + inda = idxTb[na] + pos = posTb[2*inda] + c[pos] += ca^2 + @inbounds for nb = na+1:num_coeffs_a + cb = a[nb] + TS._isthinzero(cb) && continue + indb = idxTb[nb] + pos = posTb[inda+indb] + c[pos] += interval(2) * ca * cb + end + end + + return nothing +end +# @inline function TS.sqr!(c::HomogeneousPolynomial{Complex{Interval{T}}}, +# a::HomogeneousPolynomial{Complex{Interval{T}}}) where {T<:Real} +# iszero(a) && return nothing + +# @inbounds num_coeffs_a = TS.size_table[a.order+1] + +# @inbounds posTb = TS.pos_table[c.order+1] +# @inbounds idxTb = TS.index_table[a.order+1] + +# @inbounds for na = 1:num_coeffs_a +# ca = a[na] +# isthinzero(ca) && continue +# inda = idxTb[na] +# pos = posTb[2*inda] +# c[pos] += ca^2 +# @inbounds for nb = na+1:num_coeffs_a +# cb = a[nb] +# isthinzero(cb) && continue +# indb = idxTb[nb] +# pos = posTb[inda+indb] +# c[pos] += interval(2) * ca * cb +# end +# end +# return nothing +# end + + +function sqrt(a::Taylor1{Interval{T}}) where {T<:Real} + # First non-zero coefficient + l0nz = findfirst(a) + aux = sqrt(zero( constant_term(a) )) + if l0nz < 0 + return Taylor1(aux, a.order) + elseif l0nz%2 == 1 # l0nz must be pair + throw(DomainError(a, + """First non-vanishing Taylor1 coefficient must correspond + to an **even power** in order to expand `sqrt` around 0.""")) + end + + # The last l0nz coefficients are set to zero. + lnull = l0nz >> 1 # integer division by 2 + c_order = l0nz == 0 ? a.order : a.order >> 1 + c = Taylor1( aux, c_order ) + @inbounds c[lnull] = sqrt( a[l0nz] ) + aa = one(aux) * a + for k = lnull+1:c_order + TS.sqrt!(c, aa, k, lnull) + end + + return c +end +# function sqrt(a::Taylor1{Complex{Interval{T}}}) where {T<:Real} +# # First non-zero coefficient +# l0nz = findfirst(a) +# aux = sqrt(zero( constant_term(a) )) +# if l0nz < 0 +# return Taylor1(aux, a.order) +# elseif l0nz%2 == 1 # l0nz must be pair +# throw(DomainError(a, +# """First non-vanishing Taylor1 coefficient must correspond +# to an **even power** in order to expand `sqrt` around 0.""")) +# end + +# # The last l0nz coefficients are set to zero. +# lnull = l0nz >> 1 # integer division by 2 +# c_order = l0nz == 0 ? a.order : a.order >> 1 +# c = Taylor1( aux, c_order ) +# @inbounds c[lnull] = sqrt( a[l0nz] ) +# aa = one(aux) * a +# for k = lnull+1:c_order +# TS.sqrt!(c, aa, k, lnull) +# end + +# return c +# end + +function sqrt(a::TaylorN{Interval{T}}) where {T<:Real} + @inbounds p0 = sqrt( constant_term(a) ) + if TS._isthinzero(p0) + throw(DomainError(a, + """The 0-th order TaylorN coefficient must be non-zero + in order to expand `sqrt` around 0.""")) + end + + c = TaylorN( p0, a.order) + aa = one(p0)*a + for k in 1:a.order + TS.sqrt!(c, aa, k) + end + return c +end +# function sqrt(a::TaylorN{Complex{Interval{T}}}) where {T<:Real} +# @inbounds p0 = sqrt( constant_term(a) ) +# if isthinzero(p0) +# throw(DomainError(a, +# """The 0-th order TaylorN coefficient must be non-zero +# in order to expand `sqrt` around 0.""")) +# end + +# c = TaylorN( p0, a.order) +# aa = one(p0)*a +# for k in 1:a.order +# TS.sqrt!(c, aa, k) +# end + +# return c +# end + +@inline function sqrt!(c::Taylor1{Interval{T}}, a::Taylor1{Interval{T}}, + k::Int, k0::Int=0) where {T<:Real} + if k == k0 + @inbounds c[k] = sqrt(a[2*k0]) + return nothing + end + kodd = (k - k0)%2 + kend = div(k - k0 - 2 + kodd, 2) + imax = min(k0+kend, a.order) + imin = max(k0+1, k+k0-a.order) + imin ≤ imax && ( @inbounds c[k] = c[imin] * c[k+k0-imin] ) + @inbounds for i = imin+1:imax + c[k] += c[i] * c[k+k0-i] + end + if k+k0 ≤ a.order + @inbounds aux = a[k+k0] - interval(2) * c[k] + else + @inbounds aux = - interval(2) * c[k] + end + if kodd == 0 + @inbounds aux = aux - c[kend+k0+1]^2 + end + @inbounds c[k] = aux / (interval(2) * c[k0]) + return nothing +end +# @inline function sqrt!(c::TaylorN{Interval{T}}, a::TaylorN{Interval{T}}, +# k::Int) where {T<:Real} +# if k == 0 +# @inbounds c[0] = sqrt( constant_term(a) ) +# return nothing +# end + +# kodd = k%2 +# kend = div(k - 2 + kodd, 2) +# @inbounds for i = 1:kend +# mul!(c[k], c[i], c[k-i]) +# end +# @inbounds aux = a[k] - interval(2) * c[k] +# if kodd == 0 +# @inbounds aux = aux - c[kend+1]^2 +# end +# @inbounds c[k] = aux / (interval(2) * constant_term(c)) + +# return nothing +# end + + +# several math functions for T in (:Taylor1, :TaylorN) @eval function log(a::$T{Interval{S}}) where {S<:Real} - iszero(constant_term(a)) && throw(DomainError(a, + TS._isthinzero(constant_term(a)) && throw(DomainError(a, """The 0-th order coefficient must be non-zero in order to expand `log` around 0.""")) - a0 = constant_term(a) ∩ interval(zero(S), S(Inf)) + a0 = intersect_interval(constant_term(a), interval(zero(S), S(Inf))) order = a.order aux = log(a0) aa = one(aux) * a @@ -87,16 +417,18 @@ for T in (:Taylor1, :TaylorN) end @eval function asin(a::$T{Interval{S}}) where {S<:Real} - a0 = constant_term(a) ∩ interval(-one(S), one(S)) - a0^2 == one(a0) && throw(DomainError(a, - """Series expansion of asin(x) diverges at x = ±1.""")) + a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) + a0sqr = a0^2 + isequal_interval(a0sqr, one(a0)) && + throw(DomainError(a, + """Series expansion of asin(x) diverges at x = ±1.""")) order = a.order aux = asin(a0) aa = one(aux) * a aa[0] = one(aux) * a0 c = $T( aux, order ) - r = $T( sqrt(1 - a0^2), order ) + r = $T( sqrt(1 - a0sqr), order ) for k in eachindex(a) TS.asin!(c, aa, r, k) end @@ -104,16 +436,18 @@ for T in (:Taylor1, :TaylorN) end @eval function acos(a::$T{Interval{S}}) where {S<:Real} - a0 = constant_term(a) ∩ interval(-one(S), one(S)) - a0^2 == one(a0) && throw(DomainError(a, - """Series expansion of asin(x) diverges at x = ±1.""")) + a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) + a0sqr = a0^2 + isequal_interval(a0sqr, one(a0)) && + throw(DomainError(a, + """Series expansion of acos(x) diverges at x = ±1.""")) order = a.order aux = acos(a0) aa = one(aux) * a aa[0] = one(aux) * a0 c = $T( aux, order ) - r = $T( sqrt(1 - a0^2), order ) + r = $T( sqrt(1 - a0sqr), order ) for k in eachindex(a) TS.acos!(c, aa, r, k) end @@ -121,8 +455,9 @@ for T in (:Taylor1, :TaylorN) end @eval function acosh(a::$T{Interval{S}}) where {S<:Real} - a0 = constant_term(a) ∩ interval(one(S), S(Inf)) - a0^2 == one(a0) && throw(DomainError(a, + a0 = intersect_interval(constant_term(a), interval(one(S), S(Inf))) + a0sqr = a0^2 + isequal_interval(a0sqr, one(a0)) && throw(DomainError(a, """Series expansion of acosh(x) diverges at x = ±1.""")) order = a.order @@ -130,7 +465,7 @@ for T in (:Taylor1, :TaylorN) aa = one(aux) * a aa[0] = one(aux) * a0 c = $T( aux, order ) - r = $T( sqrt(a0^2 - 1), order ) + r = $T( sqrt(a0sqr - 1), order ) for k in eachindex(a) TS.acosh!(c, aa, r, k) end @@ -139,13 +474,13 @@ for T in (:Taylor1, :TaylorN) @eval function atanh(a::$T{Interval{S}}) where {S<:Real} order = a.order - a0 = constant_term(a) ∩ interval(-one(S), one(S)) + a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) aux = atanh(a0) aa = one(aux) * a aa[0] = one(aux) * a0 c = $T( aux, order) r = $T(one(aux) - a0^2, order) - iszero(constant_term(r)) && throw(DomainError(a, + TS._isthinzero(constant_term(r)) && throw(DomainError(a, """Series expansion of atanh(x) diverges at x = ±1.""")) for k in eachindex(a) @@ -182,7 +517,6 @@ function evaluate(a::TaylorN, dx::AbstractVector{Interval{T}}) where {T<:Real} @inbounds for homPol in reverse(eachindex(a)) suma += evaluate(a.coeffs[homPol+1], dx) end - return suma end @@ -209,10 +543,31 @@ end function evaluate(a::HomogeneousPolynomial, dx::AbstractVector{Interval{T}}) where {T<:Real} @assert length(dx) == get_numvars() - all(dx .== interval(-one(T), one(T))) && return _evaluate(a, dx, Val(true)) - all(dx .== interval(zero(T), one(T))) && return _evaluate(a, dx, Val(false)) + all(isequal_interval.(dx, interval(-one(T), one(T)))) && + return _evaluate(a, dx, Val(true)) + all(isequal_interval.(dx, interval(zero(T), one(T)))) && + return _evaluate(a, dx, Val(false)) + return _evaluate(a, dx) +end + +evaluate(a::TaylorN{Taylor1}, vals::AbstractVector{Interval{T}}) where {T<:Real} = + _evaluate(a, (vals...,), Val(false)) + +# _evaluate +function _evaluate(a::HomogeneousPolynomial, dx::AbstractVector{Interval{T}}) where + {T<:Real} + a.order == 0 && return a[1] + interval(zero(T)) + + ct = TS.coeff_table[a.order+1] + @inbounds suma = a[1]*interval(zero(T)) + + for (i, a_coeff) in enumerate(a.coeffs) + TS._isthinzero(a_coeff) && continue + @inbounds tmp = prod(dx .^ ct[i]) + suma += a_coeff * tmp + end - return evaluate(a, dx...) + return suma end function _evaluate(a::HomogeneousPolynomial, dx::AbstractVector{Interval{T}}, @@ -224,15 +579,15 @@ function _evaluate(a::HomogeneousPolynomial, dx::AbstractVector{Interval{T}}, Ieven = interval(zero(T), one(T)) for (i, a_coeff) in enumerate(a.coeffs) - iszero(a_coeff) && continue + TS._isthinzero(a_coeff) && continue if isodd(sum(ct[i])) - tmp = dx[1] - else - tmp = Ieven - for n in eachindex(ct[i]) - iseven(ct[i][n]) && continue - tmp *= dx[1] - end + suma += sum(a_coeff) * dx[1] + continue + end + @inbounds tmp = iseven(ct[i][1]) ? Ieven : dx[1] + for n in 2:length(dx) + @inbounds vv = iseven(ct[i][n]) ? Ieven : dx[1] + tmp *= vv end suma += a_coeff * tmp end @@ -250,6 +605,31 @@ function _evaluate(a::HomogeneousPolynomial, dx::AbstractVector{Interval{T}}, return suma end +function _evaluate(a::TaylorN, vals::NTuple{N,TaylorN{Interval{T}}}) where {N,T<:Real} + @assert get_numvars() == N + R = promote_type(TS.numtype(a), typeof(vals[1])) + a_length = length(a) + suma = zeros(R, a_length) + @inbounds for homPol in 1:a_length + suma[homPol] = _evaluate(a.coeffs[homPol], vals) + end + return suma +end + +function _evaluate(a::HomogeneousPolynomial, + vals::NTuple{N,TaylorN{Interval{T}}}) where {N,T<:Real} + ct = TS.coeff_table[a.order+1] + suma = zero(a[1])*vals[1] + + for (i, a_coeff) in enumerate(a.coeffs) + TS._isthinzero(a_coeff) && continue + @inbounds tmp = prod( vals .^ ct[i] ) + suma += a_coeff * tmp + end + + return suma +end + """ normalize_taylor(a::Taylor1, I::Interval, symI::Bool=true) @@ -288,7 +668,8 @@ function _normalize(a::Taylor1, I::Interval{T}, ::Val{false}) where {T<:Real} end # I -> IntervalBox(-1..1, Val(N)) -function _normalize(a::TaylorN, I::AbstractVector{Interval{T}}, ::Val{true}) where {T<:Real} +function _normalize(a::TaylorN, I::AbstractVector{Interval{T}}, + ::Val{true}) where {T<:Real} order = get_order(a) x = Vector{typeof(a)}(undef, length(I)) @inbounds for ind in eachindex(x) @@ -307,4 +688,108 @@ function _normalize(a::TaylorN, I::AbstractVector{Interval{T}}, ::Val{false}) wh return a(x) end +# Printing-related methods +# function TS.pretty_print(a::Taylor1{Interval{T}}) where {T<:Real} +# z = zero(a[0]) +# var = TS._params_Taylor1_.var_name +# space = string(" ") +# bigO = TS.bigOnotation[end] ? +# string("+ 𝒪(", var, TS.superscriptify(a.order+1), ")") : +# string("") +# TS._isthinzero(a) && return string(space, z, space, bigO) +# strout::String = space +# ifirst = true +# for i in eachindex(a) +# monom::String = i==0 ? string("") : +# i==1 ? string(" ", var) : string(" ", var, TS.superscriptify(i)) +# @inbounds c = a[i] +# TS._isthinzero(c) && continue +# cadena = TS.numbr2str(c, ifirst) +# strout = string(strout, cadena, monom, space) +# ifirst = false +# end +# strout = strout * bigO +# return strout +# end + +# function TS.pretty_print(a::HomogeneousPolynomial{Interval{T}}) where +# {T<:Real} +# z = zero(a[1]) +# space = string(" ") +# TS._isthinzero(a) && return string(space, z) +# strout::String = TS.homogPol2str(a) +# return strout +# end + +# function TS.pretty_print(a::TaylorN{Interval{T}}) where {T<:Real} +# z = zero(a[0]) +# space = string("") +# bigO::String = TS.bigOnotation[end] ? +# string(" + 𝒪(‖x‖", TS.superscriptify(a.order+1), ")") : +# string("") +# TS._isthinzero(a) && return string(space, z, space, bigO) +# strout::String = space +# ifirst = true +# for ord in eachindex(a) +# pol = a[ord] +# TS._isthinzero(pol) && continue +# cadena::String = TS.homogPol2str( pol ) +# strsgn = (ifirst || ord == 0) ? +# string("") : string(" +") +# strout = string( strout, strsgn, cadena) +# ifirst = false +# end +# strout = strout * bigO +# strout +# end + + +# function TS.homogPol2str(a::HomogeneousPolynomial{Interval{T}}) where +# {T<:Real} +# numVars = get_numvars() +# order = a.order +# z = zero(a.coeffs[1]) +# space = string(" ") +# strout::String = space +# ifirst = true +# iIndices = zeros(Int, numVars) +# for pos = 1:TS.size_table[order+1] +# monom::String = string("") +# @inbounds iIndices[:] = TS.coeff_table[order+1][pos] +# for ivar = 1:numVars +# powivar = iIndices[ivar] +# if powivar == 1 +# monom = string(monom, TS.name_taylorNvar(ivar)) +# elseif powivar > 1 +# monom = string(monom, TS.name_taylorNvar(ivar), +# TS.superscriptify(powivar)) +# end +# end +# @inbounds c = a[pos] +# TS._isthinzero(c) && continue +# cadena = TS.numbr2str(c, ifirst) +# strout = string(strout, cadena, monom, space) +# ifirst = false +# end +# return strout[1:prevind(strout, end)] +# end + + +function TS.numbr2str(zz::Interval, ifirst::Bool=false) + TS._isthinzero(zz) && return string( zz ) + plusmin = ifelse( ifirst, string(""), string("+ ") ) + return string(plusmin, zz) +end +function TS.numbr2str(zz::ComplexI, ifirst::Bool=false) + zT = zero(zz.re) + TS._isthinzero(zz) && return string(zT) + if ifirst + cadena = string("( ", zz, " )") + else + cadena = string("+ ( ", zz, " )") + end + return cadena +end + + end \ No newline at end of file diff --git a/test/intervals.jl b/test/intervals.jl index 6c69163c..e853fbb5 100644 --- a/test/intervals.jl +++ b/test/intervals.jl @@ -12,12 +12,16 @@ using Test c = interval(0, 1) 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 + # ifour = interval(4) + # ifive = interval(5) + # isix = interval(6) + # iten = interval(10) + # p4(x, a) = x^4 + ifour*a*x^3 + isix*a^2*x^2 + ifour*a^3*x + a^4 + # p5(x, a) = x^5 + ifive*a*x^4 + iten*a^2*x^3 + iten*a^3*x^2 + ifive*a^4*x + a^5 ti = Taylor1(Interval{Float64}, 10) x, y = set_variables(Interval{Float64}, "x y") - # @test eltype(ti) == Interval{Float64} - # @test eltype(x) == Interval{Float64} @test eltype(ti) == Taylor1{Interval{Float64}} @test eltype(x) == TaylorN{Interval{Float64}} @test TS.numtype(ti) == Interval{Float64} @@ -28,8 +32,7 @@ using Test @test p4(ti,-a) == (ti-a)^4 @test p5(ti,-a) == (ti-a)^5 @test p4(ti,-b) == (ti-b)^4 - @test all((p5(ti,-b)).coeffs .⊆ ((ti-b)^5).coeffs) - + @test all(issubset_interval.((p5(ti,-b)).coeffs, ((ti-b)^5).coeffs)) @test p4(x,-y) == (x-y)^4 @test p5(x,-y) == (x-y)^5 @@ -37,34 +40,35 @@ using Test @test p5(x,-a) == (x-a)^5 @test p4(x,-b) == (x-b)^4 for ind in eachindex(p5(x,-b)) - @test all((p5(x,-b)[ind]).coeffs .⊆ (((x-b)^5)[ind]).coeffs) + @test all(issubset_interval.((p5(x,-b)[ind]).coeffs, (((x-b)^5)[ind]).coeffs)) end + # Tests `evaluate` - @test evaluate(p4(x,y), [a, -b]) == p4(a, -b) - @test (p5(x,y))([a, b]) == p5(a, b) - @test (a-b)^4 ⊆ ((x-y)^4)([a, b]) - @test (((x-y)^4)[4])([a, b]) == interval(-39, 81) + @test isequal_interval(evaluate(p4(x,y), [a, -b]), p4(a, -b)) + @test isequal_interval((p5(x,y))([a, b]), p5(a, b)) + @test issubset_interval((a-b)^4, ((x-y)^4)([a, b])) + @test isequal_interval((((x-y)^4)[4])([a, b]), interval(-39, 81)) p4n = normalize_taylor(p4(x,y), [a, b], true) - @test interval(0, 16) ⊆ p4n([b, b]) + @test issubset_interval(interval(0, 16), p4n([b, b])) p5n = normalize_taylor(p5(x,y), [a, b], true) - @test interval(-32, 32) ⊆ p5n([b, b]) + @test issubset_interval(interval(-32, 32), p5n([b, b])) p4n = normalize_taylor(p4(x,y), [a, b], false) - @test interval(0, 16) ⊆ p4n([c, c]) + @test issubset_interval(interval(0, 16), p4n([c, c])) p5n = normalize_taylor(p5(x,y), [a, b], false) - @test interval(0, 32) ⊆ p5n([c, c]) + @test issubset_interval(interval(0, 32), p5n([c, c])) - @test evaluate(x*y^3, [b, b]) == b - @test evaluate(x*y^2, [b, b]) == b - @test evaluate(x^2*y^2, [b, b]) == c + @test isequal_interval(evaluate(x*y^3, [b, b]), b) + @test isequal_interval(evaluate(x*y^2, [b, b]), b) + @test isequal_interval(evaluate(x^2*y^2, [b, b]), c) ii = b t = Taylor1(1) - @test interval(0, 2) ⊆ (1+t)(ii) + @test issubset_interval(interval(0, 2), (1+t)(ii)) t = Taylor1(2) - @test interval(0, 4) ⊆ ((1+t)^2)(ii) + @test issubset_interval(interval(0, 4), ((1+t)^2)(ii)) ii = interval(0, 6) t = Taylor1(4) @@ -72,12 +76,17 @@ using Test ft = f(t) f1 = normalize_taylor(ft, ii, true) f2 = normalize_taylor(ft, ii, false) - @test interval(-23/27, f(6)) ⊆ f(ii) - @test interval(-23/27, f(6)) ⊆ ft(ii) - @test interval(-23/27, f(6)) ⊆ f1b - @test interval(-23/27, f(6)) ⊆ f2(c) - @test f1b ⊆ f(ii) - @test diam(f1b) < diam(f2(c)) + @test !isguaranteed(f(ii)) + @test !isguaranteed(ft(ii)) + @test !isguaranteed(f1(b)) + @test !isguaranteed(f2(c)) + @test issubset_interval(interval(-23/27, f(6)), f(ii)) + @test issubset_interval(interval(-23/27, f(6)), ft(ii)) + @test issubset_interval(interval(-23/27, f(6)), f1(b)) + @test issubset_interval(interval(-23/27, f(6)), f2(c)) + @test issubset_interval(f1(b), f(ii)) + @test diam(f1(b)) < diam(f2(c)) + # An example from Makino's thesis ii = c @@ -85,21 +94,31 @@ using Test g(x) = 1 - x^4 + x^5 gt = g(t) g1 = normalize_taylor(gt, c, true) - @test interval(g(4/5), 1) ⊆ g(ii) - @test interval(g(4/5), 1) ⊆ gt(ii) - @test interval(g(4/5), 1) ⊆ g1b - @test g1b ⊂ g(ii) - @test diam(g1b) < diam(gt(ii)) + @test issubset_interval(interval(g(4/5), 1), g(ii)) + @test issubset_interval(interval(g(4/5), 1), gt(ii)) + @test issubset_interval(interval(g(4/5), 1), g1(b)) + @test isinterior(g1(b), g(ii)) + @test diam(g1(b)) < diam(gt(ii)) + # Test display for Taylor1{Complex{Interval{T}}} - vc = [complex(interval(1.5, 2), interval(0, 0)), complex(interval(-2, -1), interval(-1, 1 )), - complex(interval(-1, 1.5), interval(-1, 1.5)), complex( interval(0,0), interval(-1, 1.5))] + vc = [complex(interval(1.5, 2), interval(0, 0)), + complex(interval(-2, -1), interval(-1, 1 )), + complex(interval(-1, 1.5), interval(-1, 1.5)), + complex(interval(0, 0), interval(-1, 1.5))] displayBigO(false) @test string(Taylor1(vc, 5)) == - " ( Interval(1.5, 2.0) + Interval(0.0, 0.0)im ) - ( Interval(1.0, 2.0) + Interval(-1.0, 1.0)im ) t + ( Interval(-1.0, 1.5) + Interval(-1.0, 1.5)im ) t² + ( Interval(0.0, 0.0) + Interval(-1.0, 1.5)im ) t³ " + " ( Interval{Float64}(1.5, 2.0, com) + Interval{Float64}(0.0, 0.0, com)im ) +" * + " ( Interval{Float64}(-2.0, -1.0, com) + Interval{Float64}(-1.0, 1.0, com)im ) t" * + " + ( Interval{Float64}(-1.0, 1.5, com) + Interval{Float64}(-1.0, 1.5, com)im ) t²" * + " + ( Interval{Float64}(0.0, 0.0, com) + Interval{Float64}(-1.0, 1.5, com)im ) t³ " displayBigO(true) @test string(Taylor1(vc, 5)) == - " ( Interval(1.5, 2.0) + Interval(0.0, 0.0)im ) - ( Interval(1.0, 2.0) + Interval(-1.0, 1.0)im ) t + ( Interval(-1.0, 1.5) + Interval(-1.0, 1.5)im ) t² + ( Interval(0.0, 0.0) + Interval(-1.0, 1.5)im ) t³ + 𝒪(t⁶)" + " ( Interval{Float64}(1.5, 2.0, com) + Interval{Float64}(0.0, 0.0, com)im ) +" * + " ( Interval{Float64}(-2.0, -1.0, com) + Interval{Float64}(-1.0, 1.0, com)im ) t" * + " + ( Interval{Float64}(-1.0, 1.5, com) + Interval{Float64}(-1.0, 1.5, com)im ) t²" * + " + ( Interval{Float64}(0.0, 0.0, com) + Interval{Float64}(-1.0, 1.5, com)im ) t³ + 𝒪(t⁶)" + # Iss 351 (inspired by a test in ReachabilityAnalysis) p1 = Taylor1([0 .. 0, (0 .. 0.1) + (0 .. 0.01) * y], 4) @@ -115,7 +134,7 @@ using Test @test aa == sqrt(sqrt(interval(-1.e-15, 1.e-15) + ti)) bb = (interval(0.0, 1.e-15) + ti)^(1/4) @test bb == (interval(-1.e-15, 1.e-15) + ti)^(1/4) - @test all(aa.coeffs[2:end] .⊂ bb.coeffs[2:end]) + @test all(isinterior.(aa.coeffs[2:end], bb.coeffs[2:end])) @test_throws DomainError sqrt(x) @test sqrt(interval(-1,1)+x) == sqrt(interval(0,1)+x) @test (interval(-1,1)+x)^(1/4) == (interval(0,1)+x)^(1/4) From 3725c34b34a5641b0c660ca1ce9ea191a871d98a Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Tue, 12 Dec 2023 20:09:38 -0600 Subject: [PATCH 04/13] Further additions and clean up --- ext/TaylorSeriesIAExt.jl | 610 ++++++++++++++++++++++++++++----------- 1 file changed, 446 insertions(+), 164 deletions(-) diff --git a/ext/TaylorSeriesIAExt.jl b/ext/TaylorSeriesIAExt.jl index ff5ff657..a901ddff 100644 --- a/ext/TaylorSeriesIAExt.jl +++ b/ext/TaylorSeriesIAExt.jl @@ -401,92 +401,460 @@ end # several math functions for T in (:Taylor1, :TaylorN) - @eval function log(a::$T{Interval{S}}) where {S<:Real} - TS._isthinzero(constant_term(a)) && throw(DomainError(a, - """The 0-th order coefficient must be non-zero in order to expand `log` around 0.""")) - a0 = intersect_interval(constant_term(a), interval(zero(S), S(Inf))) - order = a.order - aux = log(a0) - aa = one(aux) * a - aa[0] = one(aux) * a0 - c = $T( aux, order ) - for k in eachindex(a) - TS.log!(c, aa, k) + @eval begin + function log(a::$T{Interval{S}}) where {S<:Real} + TS._isthinzero(constant_term(a)) && throw(DomainError(a, + """The 0-th order coefficient must be non-zero in order to expand `log` around 0.""")) + a0 = intersect_interval(constant_term(a), interval(zero(S), S(Inf))) + order = a.order + aux = log(a0) + aa = one(aux) * a + aa[0] = one(aux) * a0 + c = $T( aux, order ) + for k in eachindex(a) + TS.log!(c, aa, k) + end + return c end - return c - end - @eval function asin(a::$T{Interval{S}}) where {S<:Real} - a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) - a0sqr = a0^2 - isequal_interval(a0sqr, one(a0)) && - throw(DomainError(a, - """Series expansion of asin(x) diverges at x = ±1.""")) - - order = a.order - aux = asin(a0) - aa = one(aux) * a - aa[0] = one(aux) * a0 - c = $T( aux, order ) - r = $T( sqrt(1 - a0sqr), order ) - for k in eachindex(a) - TS.asin!(c, aa, r, k) + function asin(a::$T{Interval{S}}) where {S<:Real} + a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) + a0sqr = a0^2 + isequal_interval(a0sqr, one(a0)) && + throw(DomainError(a, + """Series expansion of asin(x) diverges at x = ±1.""")) + + order = a.order + aux = asin(a0) + aa = one(aux) * a + aa[0] = one(aux) * a0 + c = $T( aux, order ) + r = $T( sqrt(1 - a0sqr), order ) + for k in eachindex(a) + TS.asin!(c, aa, r, k) + end + return c end - return c - end - @eval function acos(a::$T{Interval{S}}) where {S<:Real} - a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) - a0sqr = a0^2 - isequal_interval(a0sqr, one(a0)) && - throw(DomainError(a, - """Series expansion of acos(x) diverges at x = ±1.""")) - - order = a.order - aux = acos(a0) - aa = one(aux) * a - aa[0] = one(aux) * a0 - c = $T( aux, order ) - r = $T( sqrt(1 - a0sqr), order ) - for k in eachindex(a) - TS.acos!(c, aa, r, k) + function acos(a::$T{Interval{S}}) where {S<:Real} + a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) + a0sqr = a0^2 + isequal_interval(a0sqr, one(a0)) && + throw(DomainError(a, + """Series expansion of acos(x) diverges at x = ±1.""")) + + order = a.order + aux = acos(a0) + aa = one(aux) * a + aa[0] = one(aux) * a0 + c = $T( aux, order ) + r = $T( sqrt(1 - a0sqr), order ) + for k in eachindex(a) + TS.acos!(c, aa, r, k) + end + return c end - return c - end - @eval function acosh(a::$T{Interval{S}}) where {S<:Real} - a0 = intersect_interval(constant_term(a), interval(one(S), S(Inf))) - a0sqr = a0^2 - isequal_interval(a0sqr, one(a0)) && throw(DomainError(a, - """Series expansion of acosh(x) diverges at x = ±1.""")) - - order = a.order - aux = acosh(a0) - aa = one(aux) * a - aa[0] = one(aux) * a0 - c = $T( aux, order ) - r = $T( sqrt(a0sqr - 1), order ) - for k in eachindex(a) - TS.acosh!(c, aa, r, k) + function acosh(a::$T{Interval{S}}) where {S<:Real} + a0 = intersect_interval(constant_term(a), interval(one(S), S(Inf))) + a0sqr = a0^2 + isequal_interval(a0sqr, one(a0)) && throw(DomainError(a, + """Series expansion of acosh(x) diverges at x = ±1.""")) + + order = a.order + aux = acosh(a0) + aa = one(aux) * a + aa[0] = one(aux) * a0 + c = $T( aux, order ) + r = $T( sqrt(a0sqr - 1), order ) + for k in eachindex(a) + TS.acosh!(c, aa, r, k) + end + return c + end + + function atanh(a::$T{Interval{S}}) where {S<:Real} + order = a.order + a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) + aux = atanh(a0) + aa = one(aux) * a + aa[0] = one(aux) * a0 + c = $T( aux, order) + r = $T(one(aux) - a0^2, order) + TS._isthinzero(constant_term(r)) && throw(DomainError(a, + """Series expansion of atanh(x) diverges at x = ±1.""")) + + for k in eachindex(a) + TS.atanh!(c, aa, r, k) + end + return c + end + + # Some internal functions + @inline function exp!(c::$T{Interval{T}}, a::$T{Interval{T}}, k::Int) where {T<:Real} + if k == 0 + @inbounds c[0] = exp(constant_term(a)) + return nothing + end + if $T == Taylor1 + @inbounds c[k] = interval(k) * a[k] * c[0] + else + @inbounds mul!(c[k], interval(k) * a[k], c[0]) + end + @inbounds for i = 1:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * a[k-i] * c[i] + else + mul!(c[k], interval(k-i) * a[k-i], c[i]) + end + end + @inbounds c[k] = c[k] / interval(k) + return nothing + end + + @inline function expm1!(c::$T{Interval{T}}, a::$T{Interval{T}}, k::Int) where {T<:Real} + if k == 0 + @inbounds c[0] = expm1(constant_term(a)) + return nothing + end + c0 = c[0]+one(c[0]) + if $T == Taylor1 + @inbounds c[k] = interval(k) * a[k] * c0 + else + @inbounds mul!(c[k], interval(k) * a[k], c0) + end + @inbounds for i = 1:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * a[k-i] * c[i] + else + mul!(c[k], interval(k-i) * a[k-i], c[i]) + end + end + @inbounds c[k] = c[k] / interval(k) + return nothing + end + + @inline function log!(c::$T{Interval{T}}, a::$T{T}, k::Int) where {T<:Real} + if k == 0 + @inbounds c[0] = log(constant_term(a)) + return nothing + elseif k == 1 + @inbounds c[1] = a[1] / constant_term(a) + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = interval(k-1) * a[1] * c[k-1] + else + @inbounds mul!(c[k], interval(k-1)*a[1], c[k-1]) + end + @inbounds for i = 2:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * a[i] * c[k-i] + else + mul!(c[k], interval(k-i)*a[i], c[k-i]) + end + end + @inbounds c[k] = (a[k] - c[k]/interval(k)) / constant_term(a) + return nothing + end + + @inline function log1p!(c::$T{Interval{T}}, a::$T{Interval{T}}, k::Int) where {T<:Real} + a0 = constant_term(a) + a0p1 = a0+one(a0) + if k == 0 + @inbounds c[0] = log1p(a0) + return nothing + elseif k == 1 + @inbounds c[1] = a[1] / a0p1 + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = interval(k-1) * a[1] * c[k-1] + else + @inbounds mul!(c[k], interval(k-1)*a[1], c[k-1]) + end + @inbounds for i = 2:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * a[i] * c[k-i] + else + mul!(c[k], interval(k-i)*a[i], c[k-i]) + end + end + @inbounds c[k] = (a[k] - c[k]/interval(k)) / a0p1 + return nothing + end + + @inline function sincos!(s::$T{Interval{T}}, c::$T{Interval{T}}, a::$T{Interval{T}}, + k::Int) where {T<:Real} + if k == 0 + a0 = constant_term(a) + @inbounds s[0], c[0] = sincos( a0 ) + return nothing + end + + x = a[1] + if $T == Taylor1 + @inbounds s[k] = x * c[k-1] + @inbounds c[k] = -x * s[k-1] + else + mul!(s[k], x, c[k-1]) + mul!(c[k], -x, s[k-1]) + end + @inbounds for i = 2:k + x = interval(i) * a[i] + if $T == Taylor1 + s[k] += x * c[k-i] + c[k] -= x * s[k-i] + else + mul!(s[k], x, c[k-i]) + mul!(c[k], -x, s[k-i]) + end + end + + @inbounds s[k] = s[k] / interval(k) + @inbounds c[k] = c[k] / interval(k) + return nothing + end + + @inline function tan!(c::$T{Interval{T}}, a::$T{Interval{T}}, c2::$T{Interval{T}}, + k::Int) where {T<:Real} + if k == 0 + @inbounds aux = tan( constant_term(a) ) + @inbounds c[0] = aux + @inbounds c2[0] = aux^2 + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = interval(k) * a[k] * c2[0] + else + @inbounds mul!(c[k], interval(k) * a[k], c2[0]) + end + @inbounds for i = 1:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * a[k-i] * c2[i] + else + mul!(c[k], interval(k-i) * a[k-i], c2[i]) + end + end + @inbounds c[k] = a[k] + c[k]/interval(k) + sqr!(c2, c, k) + + return nothing + end + + @inline function asin!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, + k::Int) where {T<:Real} + if k == 0 + a0 = constant_term(a) + @inbounds c[0] = asin( a0 ) + @inbounds r[0] = sqrt( one(a0) - a0^2 ) + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = interval(k-1) * r[1] * c[k-1] + else + @inbounds mul!(c[k], interval(k-1) * r[1], c[k-1]) + end + @inbounds for i in 2:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * r[i] * c[k-i] + else + mul!(c[k], interval(k-i) * r[i], c[k-i]) + end + end + sqrt!(r, one(a0)-a^2, k) + @inbounds c[k] = (a[k] - c[k]/interval(k)) / constant_term(r) + return nothing end - return c - end - @eval function atanh(a::$T{Interval{S}}) where {S<:Real} - order = a.order - a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) - aux = atanh(a0) - aa = one(aux) * a - aa[0] = one(aux) * a0 - c = $T( aux, order) - r = $T(one(aux) - a0^2, order) - TS._isthinzero(constant_term(r)) && throw(DomainError(a, - """Series expansion of atanh(x) diverges at x = ±1.""")) - - for k in eachindex(a) - TS.atanh!(c, aa, r, k) + @inline function acos!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, + k::Int) where {T<:Real} + if k == 0 + a0 = constant_term(a) + @inbounds c[0] = acos( a0 ) + @inbounds r[0] = sqrt( one(a0) - a0^2 ) + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = interval(k-1) * r[1] * c[k-1] + else + @inbounds mul!(c[k], interval(k-1) * r[1], c[k-1]) + end + @inbounds for i in 2:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * r[i] * c[k-i] + else + mul!(c[k], interval(k-i) * r[i], c[k-i]) + end + end + sqrt!(r, one(a0)-a^2, k) + @inbounds c[k] = -(a[k] + c[k]/interval(k)) / constant_term(r) + return nothing + end + + @inline function atan!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, + k::Int) where {T<:Real} + if k == 0 + a0 = constant_term(a) + @inbounds c[0] = atan( a0 ) + @inbounds r[0] = one(a0) + a0^2 + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = interval(k-1) * r[1] * c[k-1] + else + @inbounds mul!(c[k], interval(k-1) * r[1], c[k-1]) + end + @inbounds for i in 2:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * r[i] * c[k-i] + else + mul!(c[k], interval(k-i) * r[i], c[k-i]) + end + end + @inbounds sqr!(r, a, k) + @inbounds c[k] = (a[k] - c[k]/interval(k)) / constant_term(r) + return nothing + end + + @inline function sinhcosh!(s::$T{Interval{T}}, c::$T{Interval{T}}, a::$T{Interval{T}}, + k::Int) where {T<:Real} + if k == 0 + @inbounds s[0] = sinh( constant_term(a) ) + @inbounds c[0] = cosh( constant_term(a) ) + return nothing + end + + x = a[1] + if $T == Taylor1 + @inbounds s[k] = x * c[k-1] + @inbounds c[k] = x * s[k-1] + else + @inbounds mul!(s[k], x, c[k-1]) + @inbounds mul!(c[k], x, s[k-1]) + end + @inbounds for i = 2:k + x = interval(i) * a[i] + if $T == Taylor1 + s[k] += x * c[k-i] + c[k] += x * s[k-i] + else + mul!(s[k], x, c[k-i]) + mul!(c[k], x, s[k-i]) + end + end + s[k] = s[k] / interval(k) + c[k] = c[k] / interval(k) + return nothing + end + + @inline function tanh!(c::$T{Interval{T}}, a::$T{Interval{T}}, c2::$T{Interval{T}}, + k::Int) where {T<:Real} + if k == 0 + @inbounds aux = tanh( constant_term(a) ) + @inbounds c[0] = aux + @inbounds c2[0] = aux^2 + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = k * a[k] * c2[0] + else + @inbounds mul!(c[k], k * a[k], c2[0]) + end + @inbounds for i = 1:k-1 + if $T == Taylor1 + c[k] += (k-i) * a[k-i] * c2[i] + else + mul!(c[k], (k-i) * a[k-i], c2[i]) + end + end + @inbounds c[k] = a[k] - c[k]/k + sqr!(c2, c, k) + + return nothing + end + + @inline function asinh!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, k::Int) where {T<:Real} + if k == 0 + a0 = constant_term(a) + @inbounds c[0] = asinh( a0 ) + @inbounds r[0] = sqrt( a0^2 + one(a0) ) + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = interval(k-1) * r[1] * c[k-1] + else + @inbounds mul!(c[k], interval(k-1) * r[1], c[k-1]) + end + @inbounds for i in 2:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * r[i] * c[k-i] + else + mul!(c[k], interval(k-i) * r[i], c[k-i]) + end + end + sqrt!(r, a^2+one(a0), k) + @inbounds c[k] = (a[k] - c[k]/interval(k)) / constant_term(r) + return nothing + end + + @inline function acosh!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, k::Int) where {T<:Real} + if k == 0 + a0 = constant_term(a) + @inbounds c[0] = acosh( a0 ) + @inbounds r[0] = sqrt( a0^2 - one(a0) ) + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = interval(k-1) * r[1] * c[k-1] + else + @inbounds mul!(c[k], interval(k-1) * r[1], c[k-1]) + end + @inbounds for i in 2:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * r[i] * c[k-i] + else + mul!(c[k], interval(k-i) * r[i], c[k-i]) + end + end + sqrt!(r, a^2-one(a0), k) + @inbounds c[k] = (a[k] - c[k]/interval(k)) / constant_term(r) + return nothing + end + + @inline function atanh!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, k::Int) where {T<:Real} + if k == 0 + a0 = constant_term(a) + @inbounds c[0] = atanh( a0 ) + @inbounds r[0] = one(a0) - a0^2 + return nothing + end + + if $T == Taylor1 + @inbounds c[k] = interval(k-1) * r[1] * c[k-1] + else + @inbounds mul!(c[k], interval(k-1) * r[1], c[k-1]) + end + @inbounds for i in 2:k-1 + if $T == Taylor1 + c[k] += interval(k-i) * r[i] * c[k-i] + else + mul!(c[k], interval(k-i) * r[i], c[k-i]) + end + end + @inbounds sqr!(r, a, k) + @inbounds c[k] = (a[k] + c[k]/interval(k)) / constant_term(r) + return nothing end - return c end end @@ -688,93 +1056,7 @@ function _normalize(a::TaylorN, I::AbstractVector{Interval{T}}, ::Val{false}) wh return a(x) end -# Printing-related methods -# function TS.pretty_print(a::Taylor1{Interval{T}}) where {T<:Real} -# z = zero(a[0]) -# var = TS._params_Taylor1_.var_name -# space = string(" ") -# bigO = TS.bigOnotation[end] ? -# string("+ 𝒪(", var, TS.superscriptify(a.order+1), ")") : -# string("") -# TS._isthinzero(a) && return string(space, z, space, bigO) -# strout::String = space -# ifirst = true -# for i in eachindex(a) -# monom::String = i==0 ? string("") : -# i==1 ? string(" ", var) : string(" ", var, TS.superscriptify(i)) -# @inbounds c = a[i] -# TS._isthinzero(c) && continue -# cadena = TS.numbr2str(c, ifirst) -# strout = string(strout, cadena, monom, space) -# ifirst = false -# end -# strout = strout * bigO -# return strout -# end - -# function TS.pretty_print(a::HomogeneousPolynomial{Interval{T}}) where -# {T<:Real} -# z = zero(a[1]) -# space = string(" ") -# TS._isthinzero(a) && return string(space, z) -# strout::String = TS.homogPol2str(a) -# return strout -# end - -# function TS.pretty_print(a::TaylorN{Interval{T}}) where {T<:Real} -# z = zero(a[0]) -# space = string("") -# bigO::String = TS.bigOnotation[end] ? -# string(" + 𝒪(‖x‖", TS.superscriptify(a.order+1), ")") : -# string("") -# TS._isthinzero(a) && return string(space, z, space, bigO) -# strout::String = space -# ifirst = true -# for ord in eachindex(a) -# pol = a[ord] -# TS._isthinzero(pol) && continue -# cadena::String = TS.homogPol2str( pol ) -# strsgn = (ifirst || ord == 0) ? -# string("") : string(" +") -# strout = string( strout, strsgn, cadena) -# ifirst = false -# end -# strout = strout * bigO -# strout -# end - - -# function TS.homogPol2str(a::HomogeneousPolynomial{Interval{T}}) where -# {T<:Real} -# numVars = get_numvars() -# order = a.order -# z = zero(a.coeffs[1]) -# space = string(" ") -# strout::String = space -# ifirst = true -# iIndices = zeros(Int, numVars) -# for pos = 1:TS.size_table[order+1] -# monom::String = string("") -# @inbounds iIndices[:] = TS.coeff_table[order+1][pos] -# for ivar = 1:numVars -# powivar = iIndices[ivar] -# if powivar == 1 -# monom = string(monom, TS.name_taylorNvar(ivar)) -# elseif powivar > 1 -# monom = string(monom, TS.name_taylorNvar(ivar), -# TS.superscriptify(powivar)) -# end -# end -# @inbounds c = a[pos] -# TS._isthinzero(c) && continue -# cadena = TS.numbr2str(c, ifirst) -# strout = string(strout, cadena, monom, space) -# ifirst = false -# end -# return strout[1:prevind(strout, end)] -# end - - +# Printing-related methods numbr2str function TS.numbr2str(zz::Interval, ifirst::Bool=false) TS._isthinzero(zz) && return string( zz ) plusmin = ifelse( ifirst, string(""), string("+ ") ) From a796d8755a25ab45999d801d44d631ef17e64987 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Wed, 13 Dec 2023 08:28:33 -0600 Subject: [PATCH 05/13] Bump IA to v0.22 and fix a test --- Project.toml | 2 +- test/intervals.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index e135ba06..fc84b7ab 100644 --- a/Project.toml +++ b/Project.toml @@ -23,7 +23,7 @@ TaylorSeriesSAExt = "StaticArrays" [compat] Aqua = "0.8" -IntervalArithmetic = "0.15 - 0.20" +IntervalArithmetic = "0.22" JLD2 = "0.4" LinearAlgebra = "<0.0.1, 1" Markdown = "<0.0.1, 1" diff --git a/test/intervals.jl b/test/intervals.jl index e853fbb5..5ab13dc9 100644 --- a/test/intervals.jl +++ b/test/intervals.jl @@ -48,7 +48,7 @@ using Test @test isequal_interval(evaluate(p4(x,y), [a, -b]), p4(a, -b)) @test isequal_interval((p5(x,y))([a, b]), p5(a, b)) @test issubset_interval((a-b)^4, ((x-y)^4)([a, b])) - @test isequal_interval((((x-y)^4)[4])([a, b]), interval(-39, 81)) + @test isequal_interval((((x-y)^4)[4])([a, b]), interval(-64, 81)) p4n = normalize_taylor(p4(x,y), [a, b], true) @test issubset_interval(interval(0, 16), p4n([b, b])) From 1d14066303ff803d120b439ddacf3449b21782b2 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Wed, 17 Jan 2024 18:15:20 -0600 Subject: [PATCH 06/13] Fix a test, and catch empty interval warnings --- test/intervals.jl | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/test/intervals.jl b/test/intervals.jl index 5ab13dc9..16e4e18f 100644 --- a/test/intervals.jl +++ b/test/intervals.jl @@ -10,14 +10,14 @@ using Test a = interval(1, 2) b = interval(-1, 1) c = interval(0, 1) - 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 - # ifour = interval(4) - # ifive = interval(5) - # isix = interval(6) - # iten = interval(10) - # p4(x, a) = x^4 + ifour*a*x^3 + isix*a^2*x^2 + ifour*a^3*x + a^4 - # p5(x, a) = x^5 + ifive*a*x^4 + iten*a^2*x^3 + iten*a^3*x^2 + ifive*a^4*x + a^5 + # 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 + ifour = interval(4) + ifive = interval(5) + isix = interval(6) + iten = interval(10) + p4(x, a) = x^4 + ifour*a*x^3 + isix*a^2*x^2 + ifour*a^3*x + a^4 + p5(x, a) = x^5 + ifive*a*x^4 + iten*a^2*x^3 + iten*a^3*x^2 + ifive*a^4*x + a^5 ti = Taylor1(Interval{Float64}, 10) x, y = set_variables(Interval{Float64}, "x y") @@ -48,7 +48,7 @@ using Test @test isequal_interval(evaluate(p4(x,y), [a, -b]), p4(a, -b)) @test isequal_interval((p5(x,y))([a, b]), p5(a, b)) @test issubset_interval((a-b)^4, ((x-y)^4)([a, b])) - @test isequal_interval((((x-y)^4)[4])([a, b]), interval(-64, 81)) + @test isequal_interval((((x-y)^4)[4])([a, b]), interval(-39, 81)) p4n = normalize_taylor(p4(x,y), [a, b], true) @test issubset_interval(interval(0, 16), p4n([b, b])) @@ -159,8 +159,8 @@ using Test @test_throws DomainError acosh(interval(0.0, 1.0) + x) @test acosh(interval(0.0, 2.0) + x) == acosh(interval(1.0, 2.0) + x) # atanh defined on interval(-1,1) - @test_throws DomainError atanh(interval(1.0, 1.0) + ti) + @test_logs (:warn, "invalid interval, empty interval is returned") @test_throws DomainError atanh(interval(1.0, 1.0) + ti) @test atanh(interval(-2.0, 0.0) + ti) == atanh(interval(-1.0, 0.0) + ti) - @test_throws DomainError atanh(interval(1.0, 1.0) + y) + @test_logs (:warn, "invalid interval, empty interval is returned") @test_throws DomainError atanh(interval(1.0, 1.0) + y) @test atanh(interval(-2.0, 0.0) + y) == atanh(interval(-1.0, 0.0) + y) end From 9b35db238f11f5cc72a4956d46f4df498defb09a Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Wed, 31 Jan 2024 18:27:18 -0600 Subject: [PATCH 07/13] Small fixes and clean-up --- ext/TaylorSeriesIAExt.jl | 302 +++++++-------------------------------- src/power.jl | 6 +- 2 files changed, 55 insertions(+), 253 deletions(-) diff --git a/ext/TaylorSeriesIAExt.jl b/ext/TaylorSeriesIAExt.jl index a901ddff..3986d3de 100644 --- a/ext/TaylorSeriesIAExt.jl +++ b/ext/TaylorSeriesIAExt.jl @@ -41,13 +41,13 @@ for I in (:Interval, :ComplexI) return Base.power_by_squaring(a, n) end - function TS.square(a::$T{$I{T}}) where {T<:Real} - c = $T( constant_term(a)^2, a.order) - for k in 1:a.order - TS.sqr!(c, a, k) - end - return c - end + # function TS.square(a::$T{$I{T}}) where {T<:Real} + # c = $T( constant_term(a)^2, a.order) + # for k in 1:a.order + # TS.sqr!(c, a, k) + # end + # return c + # end ^(a::$T{$I{T}}, r::Rational) where {T<:Real} = a^float(r) end @@ -57,14 +57,12 @@ end function ^(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real} a0 = intersect_interval(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 TS.square(aa) r == 1/2 && return sqrt(aa) - l0 = findfirst(a) lnull = trunc(Int, r*l0 ) if (a.order-lnull < 0) || (lnull > a.order) @@ -75,86 +73,29 @@ function ^(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real} for k = 0:c_order TS.pow!(c, aa, c, r, k) end - return c end -# function ^(a::Taylor1{Complex{Interval{T}}}, r::S) where {T<:Real, S<:Real} -# a0 = intersect_interval(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 TS.square(aa) -# r == 1/2 && return sqrt(aa) - -# 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)) -# c = Taylor1(zero(aux), c_order) -# for k = 0:c_order -# TS.pow!(c, aa, r, k) -# end - -# return c -# end function ^(a::TaylorN{Interval{T}}, r::S) where {T<:Real, S<:Real} - isinteger(r) && return a^round(Int,r) + isinteger(r) && return a^Int(r) a0 = intersect_interval(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 TS.square(aa) r == 1/2 && return sqrt(aa) - # isinteger(r) && return aa^round(Int,r) - - # @assert !iszero(a0) TS._isthinzero(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) end - return c end -# function ^(a::TaylorN{Complex{Interval{T}}}, r::S) where {T<:Real, S<:Real} -# isinteger(r) && return a^round(Int,r) -# a0 = intersect_interval(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 TS.square(aa) -# r == 1/2 && return sqrt(aa) -# # isinteger(r) && return aa^round(Int,r) - -# # @assert !iszero(a0) -# TS._isthinzero(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, r, ord) -# end - -# return c -# end # sqr! @@ -166,7 +107,6 @@ for T = (:Taylor1, :TaylorN) @inbounds c[0] = constant_term(a)^2 return nothing end - kodd = k%2 kend = div(k - 2 + kodd, 2) @inbounds for i = 0:kend @@ -178,42 +118,13 @@ for T = (:Taylor1, :TaylorN) end @inbounds c[k] = interval(2) * c[k] kodd == 1 && return nothing - if $T == Taylor1 @inbounds c[k] += a[div(k,2)]^2 else TS.sqr!(c[k], a[div(k,2)]) end - return nothing end - # @inline function TS.sqr!(c::$T{Complex{Interval{T}}}, - # a::$T{Complex{Interval{T}}}, k::Int) where {T<:Real} - # if k == 0 - # @inbounds c[0] = constant_term(a)^2 - # return nothing - # end - - # kodd = k%2 - # kend = div(k - 2 + kodd, 2) - # @inbounds for i = 0:kend - # if $T == Taylor1 - # c[k] += a[i] * a[k-i] - # else - # TS.mul!(c[k], a[i], a[k-i]) - # end - # end - # @inbounds c[k] = interval(2) * c[k] - # kodd == 1 && return nothing - - # if $T == Taylor1 - # @inbounds c[k] += a[div(k,2)]^2 - # else - # TS.sqr!(c[k], a[div(k,2)]) - # end - - # return nothing - # end end end @inline function TS.sqr!(c::HomogeneousPolynomial{Interval{T}}, @@ -242,58 +153,9 @@ end return nothing end -# @inline function TS.sqr!(c::HomogeneousPolynomial{Complex{Interval{T}}}, -# a::HomogeneousPolynomial{Complex{Interval{T}}}) where {T<:Real} -# iszero(a) && return nothing - -# @inbounds num_coeffs_a = TS.size_table[a.order+1] - -# @inbounds posTb = TS.pos_table[c.order+1] -# @inbounds idxTb = TS.index_table[a.order+1] - -# @inbounds for na = 1:num_coeffs_a -# ca = a[na] -# isthinzero(ca) && continue -# inda = idxTb[na] -# pos = posTb[2*inda] -# c[pos] += ca^2 -# @inbounds for nb = na+1:num_coeffs_a -# cb = a[nb] -# isthinzero(cb) && continue -# indb = idxTb[nb] -# pos = posTb[inda+indb] -# c[pos] += interval(2) * ca * cb -# end -# end -# return nothing -# end - - -function sqrt(a::Taylor1{Interval{T}}) where {T<:Real} - # First non-zero coefficient - l0nz = findfirst(a) - aux = sqrt(zero( constant_term(a) )) - if l0nz < 0 - return Taylor1(aux, a.order) - elseif l0nz%2 == 1 # l0nz must be pair - throw(DomainError(a, - """First non-vanishing Taylor1 coefficient must correspond - to an **even power** in order to expand `sqrt` around 0.""")) - end - # The last l0nz coefficients are set to zero. - lnull = l0nz >> 1 # integer division by 2 - c_order = l0nz == 0 ? a.order : a.order >> 1 - c = Taylor1( aux, c_order ) - @inbounds c[lnull] = sqrt( a[l0nz] ) - aa = one(aux) * a - for k = lnull+1:c_order - TS.sqrt!(c, aa, k, lnull) - end - return c -end -# function sqrt(a::Taylor1{Complex{Interval{T}}}) where {T<:Real} +# function sqrt(a::Taylor1{Interval{T}}) where {T<:Real} # # First non-zero coefficient # l0nz = findfirst(a) # aux = sqrt(zero( constant_term(a) )) @@ -304,7 +166,6 @@ end # """First non-vanishing Taylor1 coefficient must correspond # to an **even power** in order to expand `sqrt` around 0.""")) # end - # # The last l0nz coefficients are set to zero. # lnull = l0nz >> 1 # integer division by 2 # c_order = l0nz == 0 ? a.order : a.order >> 1 @@ -314,29 +175,11 @@ end # for k = lnull+1:c_order # TS.sqrt!(c, aa, k, lnull) # end - # return c # end - -function sqrt(a::TaylorN{Interval{T}}) where {T<:Real} - @inbounds p0 = sqrt( constant_term(a) ) - if TS._isthinzero(p0) - throw(DomainError(a, - """The 0-th order TaylorN coefficient must be non-zero - in order to expand `sqrt` around 0.""")) - end - - c = TaylorN( p0, a.order) - aa = one(p0)*a - for k in 1:a.order - TS.sqrt!(c, aa, k) - end - - return c -end -# function sqrt(a::TaylorN{Complex{Interval{T}}}) where {T<:Real} +# function sqrt(a::TaylorN{Interval{T}}) where {T<:Real} # @inbounds p0 = sqrt( constant_term(a) ) -# if isthinzero(p0) +# if TS._isthinzero(p0) # throw(DomainError(a, # """The 0-th order TaylorN coefficient must be non-zero # in order to expand `sqrt` around 0.""")) @@ -377,39 +220,20 @@ end return nothing end -# @inline function sqrt!(c::TaylorN{Interval{T}}, a::TaylorN{Interval{T}}, -# k::Int) where {T<:Real} -# if k == 0 -# @inbounds c[0] = sqrt( constant_term(a) ) -# return nothing -# end - -# kodd = k%2 -# kend = div(k - 2 + kodd, 2) -# @inbounds for i = 1:kend -# mul!(c[k], c[i], c[k-i]) -# end -# @inbounds aux = a[k] - interval(2) * c[k] -# if kodd == 0 -# @inbounds aux = aux - c[kend+1]^2 -# end -# @inbounds c[k] = aux / (interval(2) * constant_term(c)) - -# return nothing -# end # several math functions for T in (:Taylor1, :TaylorN) @eval begin function log(a::$T{Interval{S}}) where {S<:Real} - TS._isthinzero(constant_term(a)) && throw(DomainError(a, - """The 0-th order coefficient must be non-zero in order to expand `log` around 0.""")) a0 = intersect_interval(constant_term(a), interval(zero(S), S(Inf))) - order = a.order aux = log(a0) - aa = one(aux) * a - aa[0] = one(aux) * a0 + isempty_interval(aux) && throw(DomainError(a, + """The 0-th order coefficient must be positive in order to expand `log` around 0.""")) + order = a.order + uno = one(aux) + aa = uno * a + aa[0] = uno * a0 c = $T( aux, order ) for k in eachindex(a) TS.log!(c, aa, k) @@ -420,16 +244,15 @@ for T in (:Taylor1, :TaylorN) function asin(a::$T{Interval{S}}) where {S<:Real} a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) a0sqr = a0^2 - isequal_interval(a0sqr, one(a0)) && - throw(DomainError(a, + uno = one(a0) + isequal_interval(a0sqr, uno) && throw(DomainError(a, """Series expansion of asin(x) diverges at x = ±1.""")) - order = a.order aux = asin(a0) - aa = one(aux) * a - aa[0] = one(aux) * a0 + aa = uno * a + aa[0] = uno * a0 c = $T( aux, order ) - r = $T( sqrt(1 - a0sqr), order ) + r = $T( sqrt(uno - a0sqr), order ) for k in eachindex(a) TS.asin!(c, aa, r, k) end @@ -439,16 +262,15 @@ for T in (:Taylor1, :TaylorN) function acos(a::$T{Interval{S}}) where {S<:Real} a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) a0sqr = a0^2 - isequal_interval(a0sqr, one(a0)) && - throw(DomainError(a, + uno = one(a0) + isequal_interval(a0sqr, uno) && throw(DomainError(a, """Series expansion of acos(x) diverges at x = ±1.""")) - order = a.order aux = acos(a0) - aa = one(aux) * a - aa[0] = one(aux) * a0 + aa = uno * a + aa[0] = uno * a0 c = $T( aux, order ) - r = $T( sqrt(1 - a0sqr), order ) + r = $T( sqrt(uno - a0sqr), order ) for k in eachindex(a) TS.acos!(c, aa, r, k) end @@ -458,15 +280,15 @@ for T in (:Taylor1, :TaylorN) function acosh(a::$T{Interval{S}}) where {S<:Real} a0 = intersect_interval(constant_term(a), interval(one(S), S(Inf))) a0sqr = a0^2 - isequal_interval(a0sqr, one(a0)) && throw(DomainError(a, + uno = one(a0) + isequal_interval(a0sqr, uno) && throw(DomainError(a, """Series expansion of acosh(x) diverges at x = ±1.""")) - order = a.order aux = acosh(a0) - aa = one(aux) * a - aa[0] = one(aux) * a0 + aa = uno * a + aa[0] = uno * a0 c = $T( aux, order ) - r = $T( sqrt(a0sqr - 1), order ) + r = $T( sqrt(a0sqr - uno), order ) for k in eachindex(a) TS.acosh!(c, aa, r, k) end @@ -476,14 +298,14 @@ for T in (:Taylor1, :TaylorN) function atanh(a::$T{Interval{S}}) where {S<:Real} order = a.order a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) + uno = one(a0) aux = atanh(a0) - aa = one(aux) * a - aa[0] = one(aux) * a0 + aa = uno * a + aa[0] = uno * a0 c = $T( aux, order) - r = $T(one(aux) - a0^2, order) + r = $T(uno - a0^2, order) TS._isthinzero(constant_term(r)) && throw(DomainError(a, """Series expansion of atanh(x) diverges at x = ±1.""")) - for k in eachindex(a) TS.atanh!(c, aa, r, k) end @@ -542,7 +364,6 @@ for T in (:Taylor1, :TaylorN) @inbounds c[1] = a[1] / constant_term(a) return nothing end - if $T == Taylor1 @inbounds c[k] = interval(k-1) * a[1] * c[k-1] else @@ -569,7 +390,6 @@ for T in (:Taylor1, :TaylorN) @inbounds c[1] = a[1] / a0p1 return nothing end - if $T == Taylor1 @inbounds c[k] = interval(k-1) * a[1] * c[k-1] else @@ -593,7 +413,6 @@ for T in (:Taylor1, :TaylorN) @inbounds s[0], c[0] = sincos( a0 ) return nothing end - x = a[1] if $T == Taylor1 @inbounds s[k] = x * c[k-1] @@ -612,7 +431,6 @@ for T in (:Taylor1, :TaylorN) mul!(c[k], -x, s[k-i]) end end - @inbounds s[k] = s[k] / interval(k) @inbounds c[k] = c[k] / interval(k) return nothing @@ -626,7 +444,6 @@ for T in (:Taylor1, :TaylorN) @inbounds c2[0] = aux^2 return nothing end - if $T == Taylor1 @inbounds c[k] = interval(k) * a[k] * c2[0] else @@ -641,7 +458,6 @@ for T in (:Taylor1, :TaylorN) end @inbounds c[k] = a[k] + c[k]/interval(k) sqr!(c2, c, k) - return nothing end @@ -653,7 +469,6 @@ for T in (:Taylor1, :TaylorN) @inbounds r[0] = sqrt( one(a0) - a0^2 ) return nothing end - if $T == Taylor1 @inbounds c[k] = interval(k-1) * r[1] * c[k-1] else @@ -679,7 +494,6 @@ for T in (:Taylor1, :TaylorN) @inbounds r[0] = sqrt( one(a0) - a0^2 ) return nothing end - if $T == Taylor1 @inbounds c[k] = interval(k-1) * r[1] * c[k-1] else @@ -705,7 +519,6 @@ for T in (:Taylor1, :TaylorN) @inbounds r[0] = one(a0) + a0^2 return nothing end - if $T == Taylor1 @inbounds c[k] = interval(k-1) * r[1] * c[k-1] else @@ -730,7 +543,6 @@ for T in (:Taylor1, :TaylorN) @inbounds c[0] = cosh( constant_term(a) ) return nothing end - x = a[1] if $T == Taylor1 @inbounds s[k] = x * c[k-1] @@ -762,7 +574,6 @@ for T in (:Taylor1, :TaylorN) @inbounds c2[0] = aux^2 return nothing end - if $T == Taylor1 @inbounds c[k] = k * a[k] * c2[0] else @@ -777,7 +588,6 @@ for T in (:Taylor1, :TaylorN) end @inbounds c[k] = a[k] - c[k]/k sqr!(c2, c, k) - return nothing end @@ -788,7 +598,6 @@ for T in (:Taylor1, :TaylorN) @inbounds r[0] = sqrt( a0^2 + one(a0) ) return nothing end - if $T == Taylor1 @inbounds c[k] = interval(k-1) * r[1] * c[k-1] else @@ -813,7 +622,6 @@ for T in (:Taylor1, :TaylorN) @inbounds r[0] = sqrt( a0^2 - one(a0) ) return nothing end - if $T == Taylor1 @inbounds c[k] = interval(k-1) * r[1] * c[k-1] else @@ -838,7 +646,6 @@ for T in (:Taylor1, :TaylorN) @inbounds r[0] = one(a0) - a0^2 return nothing end - if $T == Taylor1 @inbounds c[k] = interval(k-1) * r[1] * c[k-1] else @@ -879,7 +686,8 @@ function evaluate(a::Taylor1, dx::Interval{S}) where {S<:Real} return sum_even + sum_odd*dx end -function evaluate(a::TaylorN, dx::AbstractVector{Interval{T}}) where {T<:Real} +function evaluate(a::TaylorN{T}, dx::AbstractVector{Interval{S}}) where + {T<:Real, S<:Real} @assert length(dx) == get_numvars() suma = zero(constant_term(a)) + interval(zero(T)) @inbounds for homPol in reverse(eachindex(a)) @@ -909,7 +717,8 @@ function evaluate(a::Taylor1{TaylorN{T}}, dx::Interval{S}) where {T<:Real, S<:Re end -function evaluate(a::HomogeneousPolynomial, dx::AbstractVector{Interval{T}}) where {T<:Real} +function evaluate(a::HomogeneousPolynomial{T}, + dx::AbstractVector{Interval{S}}) where {T<:Real, S<:Real} @assert length(dx) == get_numvars() all(isequal_interval.(dx, interval(-one(T), one(T)))) && return _evaluate(a, dx, Val(true)) @@ -918,33 +727,28 @@ function evaluate(a::HomogeneousPolynomial, dx::AbstractVector{Interval{T}}) whe return _evaluate(a, dx) end -evaluate(a::TaylorN{Taylor1}, vals::AbstractVector{Interval{T}}) where {T<:Real} = - _evaluate(a, (vals...,), Val(false)) +evaluate(a::TaylorN{Taylor1{T}}, vals::AbstractVector{Interval{S}}) where + {T<:Real, S<:Real} = _evaluate(a, (vals...,), Val(false)) # _evaluate -function _evaluate(a::HomogeneousPolynomial, dx::AbstractVector{Interval{T}}) where - {T<:Real} +function _evaluate(a::HomogeneousPolynomial{T}, + dx::AbstractVector{Interval{S}}) where {T<:Real, S<:Real} a.order == 0 && return a[1] + interval(zero(T)) - ct = TS.coeff_table[a.order+1] @inbounds suma = a[1]*interval(zero(T)) - for (i, a_coeff) in enumerate(a.coeffs) TS._isthinzero(a_coeff) && continue @inbounds tmp = prod(dx .^ ct[i]) suma += a_coeff * tmp end - return suma end -function _evaluate(a::HomogeneousPolynomial, dx::AbstractVector{Interval{T}}, - ::Val{true} ) where {T<:Real} +function _evaluate(a::HomogeneousPolynomial{T}, dx::AbstractVector{Interval{S}}, + ::Val{true} ) where {T<:Real, S<:Real} a.order == 0 && return a[1] + interval(zero(T)) - ct = TS.coeff_table[a.order+1] @inbounds suma = a[1]*interval(zero(T)) - Ieven = interval(zero(T), one(T)) for (i, a_coeff) in enumerate(a.coeffs) TS._isthinzero(a_coeff) && continue @@ -962,10 +766,9 @@ function _evaluate(a::HomogeneousPolynomial, dx::AbstractVector{Interval{T}}, return suma end -function _evaluate(a::HomogeneousPolynomial, dx::AbstractVector{Interval{T}}, - ::Val{false} ) where {T<:Real} +function _evaluate(a::HomogeneousPolynomial{T}, dx::AbstractVector{Interval{S}}, + ::Val{false} ) where {T<:Real, S<:Real} a.order == 0 && return a[1] + interval(zero(T)) - @inbounds suma = zero(a[1])*dx[1] @inbounds for homPol in a.coeffs suma += homPol*dx[1] @@ -973,7 +776,8 @@ function _evaluate(a::HomogeneousPolynomial, dx::AbstractVector{Interval{T}}, return suma end -function _evaluate(a::TaylorN, vals::NTuple{N,TaylorN{Interval{T}}}) where {N,T<:Real} +function _evaluate(a::TaylorN{T}, vals::NTuple{N,TaylorN{Interval{S}}}) where + {N, T<:Real, S<:Real} @assert get_numvars() == N R = promote_type(TS.numtype(a), typeof(vals[1])) a_length = length(a) @@ -984,17 +788,15 @@ function _evaluate(a::TaylorN, vals::NTuple{N,TaylorN{Interval{T}}}) where {N,T< return suma end -function _evaluate(a::HomogeneousPolynomial, - vals::NTuple{N,TaylorN{Interval{T}}}) where {N,T<:Real} +function _evaluate(a::HomogeneousPolynomial{T}, + vals::NTuple{N,TaylorN{Interval{S}}}) where {N, T<:Real, S<:Real} ct = TS.coeff_table[a.order+1] suma = zero(a[1])*vals[1] - for (i, a_coeff) in enumerate(a.coeffs) TS._isthinzero(a_coeff) && continue @inbounds tmp = prod( vals .^ ct[i] ) suma += a_coeff * tmp end - return suma end diff --git a/src/power.jl b/src/power.jl index 67ca6726..721a43ad 100644 --- a/src/power.jl +++ b/src/power.jl @@ -593,13 +593,13 @@ Returns `c += a*a` with no allocation; all parameters are `HomogeneousPolynomial @inbounds for na = 1:num_coeffs_a ca = a[na] - iszero(ca) && continue + _isthinzero(ca) && continue inda = idxTb[na] pos = posTb[2*inda] c[pos] += ca^2 @inbounds for nb = na+1:num_coeffs_a cb = a[nb] - iszero(cb) && continue + _isthinzero(cb) && continue indb = idxTb[nb] pos = posTb[inda+indb] c[pos] += 2 * ca * cb @@ -638,7 +638,7 @@ end function sqrt(a::TaylorN) @inbounds p0 = sqrt( constant_term(a) ) - if iszero(p0) + if TS._isthinzero(p0) throw(DomainError(a, """The 0-th order TaylorN coefficient must be non-zero in order to expand `sqrt` around 0.""")) From 39e33372ae74ceb5996c09b9865ee5012ae4d379 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Fri, 9 Feb 2024 11:42:50 -0600 Subject: [PATCH 08/13] Fixes in functions for Interval's --- ext/TaylorSeriesIAExt.jl | 439 +++++++++++++++++++++++---------------- src/power.jl | 2 +- test/intervals.jl | 30 ++- 3 files changed, 283 insertions(+), 188 deletions(-) diff --git a/ext/TaylorSeriesIAExt.jl b/ext/TaylorSeriesIAExt.jl index 3986d3de..02b13fc8 100644 --- a/ext/TaylorSeriesIAExt.jl +++ b/ext/TaylorSeriesIAExt.jl @@ -8,32 +8,59 @@ import TaylorSeries: evaluate, _evaluate, normalize_taylor isdefined(Base, :get_extension) ? (using IntervalArithmetic) : (using ..IntervalArithmetic) +const NumTypes = IntervalArithmetic.NumTypes + +""" + intersect_interval_nonstd(x, y) + +Returns the intersection of the intervals `x` and `y`, considered as (extended) +sets of real numbers. That is, the set that contains the points common in `x` +and `y`. + +This function is similar to [`intersect_interval`](@ref), but allows for higher +than `trv` decorations depending on `issubset_interval(x, y)`: if +`issubset_interval(x,y) == true` the decoration is that of the intersection +of the bare intervals, otherwise it is `trv` (Section 11.7.1). +""" +_intersect_domain_nonstd(x::BareInterval, y::BareInterval) = intersect_interval(x, y) + +function _intersect_domain_nonstd(x::Interval{T}, y::Interval{T}) where {T<:NumTypes} + d = ifelse(issubset_interval(x, y), decoration(x), trv) + x = intersect_interval(x, y) + return IntervalArithmetic._unsafe_interval(bareinterval(x), d, isguaranteed(x)) +end + +_intersect_domain_nonstd(x::Interval, y::Interval) = + _intersect_domain_nonstd(promote(x, y)...) + # Some functions require special interval functions (isequal_interval, isthinzero) for I in (:Interval, :ComplexI) @eval begin TS._isthinzero(x::$I{T}) where {T<:Real} = isthinzero(x) - function ==(a::Taylor1{$I{T}}, b::Taylor1{$I{S}}) where {T<:Real, S<:Real} + function ==(a::Taylor1{$I{T}}, b::Taylor1{$I{S}}) where {T<:NumTypes, S<:NumTypes} if a.order != b.order a, b = fixorder(a, b) end return all(isequal_interval.(a.coeffs, b.coeffs)) end + function ==(a::HomogeneousPolynomial{$I{T}}, - b::HomogeneousPolynomial{$I{S}}) where {T<:Real, S<:Real} + b::HomogeneousPolynomial{$I{S}}) where {T<:NumTypes, S<:NumTypes} a.order == b.order && return all(isequal_interval.(a.coeffs, b.coeffs)) return all(TS._isthinzero.(a.coeffs)) && all(TS._isthinzero.(b.coeffs)) end - iszero(a::Taylor1{$I{T}}) where {T<:Real} = all(TS._isthinzero.(a.coeffs)) - iszero(a::HomogeneousPolynomial{$I{T}}) where {T<:Real} = + iszero(a::Taylor1{$I{T}}) where {T<:NumTypes} = all(TS._isthinzero.(a.coeffs)) + + iszero(a::HomogeneousPolynomial{$I{T}}) where {T<:NumTypes} = all(TS._isthinzero.(a.coeffs)) end # Methods related to power, sqr, sqrt, ... for T in (:Taylor1, :TaylorN) @eval begin - function ^(a::$T{$I{T}}, n::S) where {T<:Real, S<:Integer} + function ^(a::$T{$I{T}}, n::S) where {T<:NumTypes, S<:Integer} n == 0 && return one(a) n == 1 && return copy(a) n == 2 && return TS.square(a) @@ -41,25 +68,19 @@ for I in (:Interval, :ComplexI) return Base.power_by_squaring(a, n) end - # function TS.square(a::$T{$I{T}}) where {T<:Real} - # c = $T( constant_term(a)^2, a.order) - # for k in 1:a.order - # TS.sqr!(c, a, k) - # end - # return c - # end - - ^(a::$T{$I{T}}, r::Rational) where {T<:Real} = a^float(r) + ^(a::$T{$I{T}}, r::Rational) where {T<:NumTypes} = a^float(r) end end end -function ^(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real} +function ^(a::Taylor1{Interval{T}}, r::S) where {T<:NumTypes, S<:Real} a0 = intersect_interval(constant_term(a), interval(zero(T), T(Inf))) aux = one(a0)^r iszero(r) && return Taylor1(aux, a.order) - aa = one(aux) * a + # aa = one(aux) * a + aa = convert(Taylor1{typeof(aux)}, a) aa[0] = one(aux) * a0 + @show(aa) r == 1 && return aa r == 2 && return TS.square(aa) r == 1/2 && return sqrt(aa) @@ -70,13 +91,13 @@ function ^(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real} end 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) + for k in eachindex(c) + TS.pow!(c, aa, r, k) end return c end -function ^(a::TaylorN{Interval{T}}, r::S) where {T<:Real, S<:Real} +function ^(a::TaylorN{Interval{T}}, r::S) where {T<:NumTypes, S<:Real} isinteger(r) && return a^Int(r) a0 = intersect_interval(constant_term(a), interval(zero(T), T(Inf))) a0r = a0^r @@ -102,7 +123,7 @@ end for T = (:Taylor1, :TaylorN) @eval begin @inline function TS.sqr!(c::$T{Interval{T}}, a::$T{Interval{T}}, - k::Int) where {T<:Real} + k::Int) where {T<:NumTypes} if k == 0 @inbounds c[0] = constant_term(a)^2 return nothing @@ -128,14 +149,11 @@ for T = (:Taylor1, :TaylorN) end end @inline function TS.sqr!(c::HomogeneousPolynomial{Interval{T}}, - a::HomogeneousPolynomial{Interval{T}}) where {T<:Real} + a::HomogeneousPolynomial{Interval{T}}) where {T<:NumTypes} iszero(a) && return nothing - @inbounds num_coeffs_a = TS.size_table[a.order+1] - @inbounds posTb = TS.pos_table[c.order+1] @inbounds idxTb = TS.index_table[a.order+1] - @inbounds for na = 1:num_coeffs_a ca = a[na] TS._isthinzero(ca) && continue @@ -150,52 +168,60 @@ end c[pos] += interval(2) * ca * cb end end - return nothing end -# function sqrt(a::Taylor1{Interval{T}}) where {T<:Real} -# # First non-zero coefficient -# l0nz = findfirst(a) -# aux = sqrt(zero( constant_term(a) )) -# if l0nz < 0 -# return Taylor1(aux, a.order) -# elseif l0nz%2 == 1 # l0nz must be pair -# throw(DomainError(a, -# """First non-vanishing Taylor1 coefficient must correspond -# to an **even power** in order to expand `sqrt` around 0.""")) -# end -# # The last l0nz coefficients are set to zero. -# lnull = l0nz >> 1 # integer division by 2 -# c_order = l0nz == 0 ? a.order : a.order >> 1 -# c = Taylor1( aux, c_order ) -# @inbounds c[lnull] = sqrt( a[l0nz] ) -# aa = one(aux) * a -# for k = lnull+1:c_order -# TS.sqrt!(c, aa, k, lnull) -# end -# return c -# end -# function sqrt(a::TaylorN{Interval{T}}) where {T<:Real} -# @inbounds p0 = sqrt( constant_term(a) ) -# if TS._isthinzero(p0) -# throw(DomainError(a, -# """The 0-th order TaylorN coefficient must be non-zero -# in order to expand `sqrt` around 0.""")) -# end - -# c = TaylorN( p0, a.order) -# aa = one(p0)*a -# for k in 1:a.order -# TS.sqrt!(c, aa, k) -# end - -# return c -# end - -@inline function sqrt!(c::Taylor1{Interval{T}}, a::Taylor1{Interval{T}}, - k::Int, k0::Int=0) where {T<:Real} +function sqrt(a::Taylor1{Interval{T}}) where {T<:NumTypes} + domain = interval(zero(T), typemax(T)) + a0 = _intersect_domain_nonstd(constant_term(a), domain) + aux = sqrt(a0) + isempty_interval(aux) && throw(DomainError(a, + """The 0-th order coefficient must have a positive part + in order to expand `sqrt`.""")) + # First non-zero coefficient + aa = convert(Taylor1{typeof(aux)}, a) + aa[0] = one(aux)*a0 + l0nz = findfirst(aa) + order = a.order + if l0nz < 0 + return Taylor1(zero(aux), order) + elseif l0nz%2 == 1 # l0nz must be pair + throw(DomainError(aa, + """First non-vanishing Taylor1 coefficient must correspond + to a **even power** in order to expand `sqrt`.""")) + end + # The last l0nz coefficients are set to zero. + lnull = l0nz >> 1 # integer division by 2 + c_order = l0nz == 0 ? order : order >> 1 + c = Taylor1( zero(aux), c_order ) + @inbounds c[lnull] = aux + for k = lnull+1:c_order + TS.sqrt!(c, aa, k, lnull) + end + return c +end + +function sqrt(a::TaylorN{Interval{T}}) where {T<:NumTypes} + domain = interval(zero(T), typemax(T)) + a0 = _intersect_domain_nonstd(constant_term(a), domain) + aux = sqrt(a0) + (isempty_interval(aux) || TS._isthinzero(a0)) && throw(DomainError(a, + """The 0-th order coefficient must have a positive part + in order to expand `sqrt`.""")) + # First non-zero coefficient + aa = convert(TaylorN{typeof(aux)}, a) + aa[0] = one(aux)*a0 + order = a.order + c = TaylorN( zero(aux), order) + for k in 1:order + TS.sqrt!(c, aa, k) + end + return c +end + +@inline function TS.sqrt!(c::Taylor1{Interval{T}}, a::Taylor1{Interval{T}}, + k::Int, k0::Int=0) where {T<:AbstractFloat} if k == k0 @inbounds c[k] = sqrt(a[2*k0]) return nothing @@ -217,23 +243,22 @@ end @inbounds aux = aux - c[kend+k0+1]^2 end @inbounds c[k] = aux / (interval(2) * c[k0]) - return nothing end - +# Faltan métodos de sqrt! para TaylorN y Taylor1{TaylorN} # several math functions for T in (:Taylor1, :TaylorN) @eval begin - function log(a::$T{Interval{S}}) where {S<:Real} - a0 = intersect_interval(constant_term(a), interval(zero(S), S(Inf))) + function log(a::$T{Interval{T}}) where {T<:NumTypes} + domain = interval(zero(T), typemax(T)) + a0 = _intersect_domain_nonstd(constant_term(a), domain) aux = log(a0) isempty_interval(aux) && throw(DomainError(a, - """The 0-th order coefficient must be positive in order to expand `log` around 0.""")) + """The 0-th order coefficient must be positive in order to expand `log`.""")) + aa = convert($T{typeof(aux)}, a) + aa[0] = one(aux) * a0 order = a.order - uno = one(aux) - aa = uno * a - aa[0] = uno * a0 c = $T( aux, order ) for k in eachindex(a) TS.log!(c, aa, k) @@ -241,15 +266,18 @@ for T in (:Taylor1, :TaylorN) return c end - function asin(a::$T{Interval{S}}) where {S<:Real} - a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) + function asin(a::$T{Interval{T}}) where {T<:NumTypes} + domain = interval(-one(T), one(T)) + a0 = _intersect_domain_nonstd(constant_term(a), domain) + aux = asin(a0) + isempty_interval(aux) && throw(DomainError(a, + """The 0-th order coefficient must have a non-empty intersection with $domain.""")) a0sqr = a0^2 - uno = one(a0) + uno = one(aux) isequal_interval(a0sqr, uno) && throw(DomainError(a, """Series expansion of asin(x) diverges at x = ±1.""")) order = a.order - aux = asin(a0) - aa = uno * a + aa = convert($T{typeof(aux)}, a) aa[0] = uno * a0 c = $T( aux, order ) r = $T( sqrt(uno - a0sqr), order ) @@ -259,15 +287,18 @@ for T in (:Taylor1, :TaylorN) return c end - function acos(a::$T{Interval{S}}) where {S<:Real} - a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) + function acos(a::$T{Interval{T}}) where {T<:NumTypes} + domain = interval(-one(T), one(T)) + a0 = _intersect_domain_nonstd(constant_term(a), domain) + aux = acos(a0) + isempty_interval(aux) && throw(DomainError(a, + """The 0-th order coefficient must have a non-empty intersection with $domain.""")) a0sqr = a0^2 uno = one(a0) isequal_interval(a0sqr, uno) && throw(DomainError(a, """Series expansion of acos(x) diverges at x = ±1.""")) order = a.order - aux = acos(a0) - aa = uno * a + aa = convert($T{typeof(aux)}, a) aa[0] = uno * a0 c = $T( aux, order ) r = $T( sqrt(uno - a0sqr), order ) @@ -277,15 +308,18 @@ for T in (:Taylor1, :TaylorN) return c end - function acosh(a::$T{Interval{S}}) where {S<:Real} - a0 = intersect_interval(constant_term(a), interval(one(S), S(Inf))) + function acosh(a::$T{Interval{T}}) where {T<:NumTypes} + domain = interval(one(T), typemax(T)) + a0 = _intersect_domain_nonstd(constant_term(a), domain) + aux = acosh(a0) + isempty_interval(aux) && throw(DomainError(a, + """The 0-th order coefficient must have a non-empty intersection with $domain.""")) a0sqr = a0^2 uno = one(a0) isequal_interval(a0sqr, uno) && throw(DomainError(a, """Series expansion of acosh(x) diverges at x = ±1.""")) order = a.order - aux = acosh(a0) - aa = uno * a + aa = convert($T{typeof(aux)}, a) aa[0] = uno * a0 c = $T( aux, order ) r = $T( sqrt(a0sqr - uno), order ) @@ -295,12 +329,15 @@ for T in (:Taylor1, :TaylorN) return c end - function atanh(a::$T{Interval{S}}) where {S<:Real} + function atanh(a::$T{Interval{T}}) where {T<:NumTypes} + domain = interval(-one(T), one(T)) + a0 = _intersect_domain_nonstd(constant_term(a), domain) + aux = atanh(a0) + isempty_interval(aux) && throw(DomainError(a, + """The 0-th order coefficient must have a non-empty intersection with $domain.""")) order = a.order - a0 = intersect_interval(constant_term(a), interval(-one(S), one(S))) uno = one(a0) - aux = atanh(a0) - aa = uno * a + aa = convert($T{typeof(aux)}, a) aa[0] = uno * a0 c = $T( aux, order) r = $T(uno - a0^2, order) @@ -313,7 +350,7 @@ for T in (:Taylor1, :TaylorN) end # Some internal functions - @inline function exp!(c::$T{Interval{T}}, a::$T{Interval{T}}, k::Int) where {T<:Real} + @inline function TS.exp!(c::$T{Interval{T}}, a::$T{Interval{T}}, k::Int) where {T<:AbstractFloat} if k == 0 @inbounds c[0] = exp(constant_term(a)) return nothing @@ -321,20 +358,20 @@ for T in (:Taylor1, :TaylorN) if $T == Taylor1 @inbounds c[k] = interval(k) * a[k] * c[0] else - @inbounds mul!(c[k], interval(k) * a[k], c[0]) + @inbounds TS.mul!(c[k], interval(k) * a[k], c[0]) end @inbounds for i = 1:k-1 if $T == Taylor1 c[k] += interval(k-i) * a[k-i] * c[i] else - mul!(c[k], interval(k-i) * a[k-i], c[i]) + TS.mul!(c[k], interval(k-i) * a[k-i], c[i]) end end @inbounds c[k] = c[k] / interval(k) return nothing end - @inline function expm1!(c::$T{Interval{T}}, a::$T{Interval{T}}, k::Int) where {T<:Real} + @inline function TS.expm1!(c::$T{Interval{T}}, a::$T{Interval{T}}, k::Int) where {T<:AbstractFloat} if k == 0 @inbounds c[0] = expm1(constant_term(a)) return nothing @@ -343,20 +380,20 @@ for T in (:Taylor1, :TaylorN) if $T == Taylor1 @inbounds c[k] = interval(k) * a[k] * c0 else - @inbounds mul!(c[k], interval(k) * a[k], c0) + @inbounds TS.mul!(c[k], interval(k) * a[k], c0) end @inbounds for i = 1:k-1 if $T == Taylor1 c[k] += interval(k-i) * a[k-i] * c[i] else - mul!(c[k], interval(k-i) * a[k-i], c[i]) + TS.mul!(c[k], interval(k-i) * a[k-i], c[i]) end end @inbounds c[k] = c[k] / interval(k) return nothing end - @inline function log!(c::$T{Interval{T}}, a::$T{T}, k::Int) where {T<:Real} + @inline function TS.log!(c::$T{Interval{T}}, a::$T{Interval{T}}, k::Int) where {T<:AbstractFloat} if k == 0 @inbounds c[0] = log(constant_term(a)) return nothing @@ -367,20 +404,20 @@ for T in (:Taylor1, :TaylorN) if $T == Taylor1 @inbounds c[k] = interval(k-1) * a[1] * c[k-1] else - @inbounds mul!(c[k], interval(k-1)*a[1], c[k-1]) + @inbounds TS.mul!(c[k], interval(k-1)*a[1], c[k-1]) end @inbounds for i = 2:k-1 if $T == Taylor1 c[k] += interval(k-i) * a[i] * c[k-i] else - mul!(c[k], interval(k-i)*a[i], c[k-i]) + TS.mul!(c[k], interval(k-i)*a[i], c[k-i]) end end @inbounds c[k] = (a[k] - c[k]/interval(k)) / constant_term(a) return nothing end - @inline function log1p!(c::$T{Interval{T}}, a::$T{Interval{T}}, k::Int) where {T<:Real} + @inline function TS.log1p!(c::$T{Interval{T}}, a::$T{Interval{T}}, k::Int) where {T<:AbstractFloat} a0 = constant_term(a) a0p1 = a0+one(a0) if k == 0 @@ -393,21 +430,21 @@ for T in (:Taylor1, :TaylorN) if $T == Taylor1 @inbounds c[k] = interval(k-1) * a[1] * c[k-1] else - @inbounds mul!(c[k], interval(k-1)*a[1], c[k-1]) + @inbounds TS.mul!(c[k], interval(k-1)*a[1], c[k-1]) end @inbounds for i = 2:k-1 if $T == Taylor1 c[k] += interval(k-i) * a[i] * c[k-i] else - mul!(c[k], interval(k-i)*a[i], c[k-i]) + TS.mul!(c[k], interval(k-i)*a[i], c[k-i]) end end @inbounds c[k] = (a[k] - c[k]/interval(k)) / a0p1 return nothing end - @inline function sincos!(s::$T{Interval{T}}, c::$T{Interval{T}}, a::$T{Interval{T}}, - k::Int) where {T<:Real} + @inline function TS.sincos!(s::$T{Interval{T}}, c::$T{Interval{T}}, a::$T{Interval{T}}, + k::Int) where {T<:AbstractFloat} if k == 0 a0 = constant_term(a) @inbounds s[0], c[0] = sincos( a0 ) @@ -418,8 +455,8 @@ for T in (:Taylor1, :TaylorN) @inbounds s[k] = x * c[k-1] @inbounds c[k] = -x * s[k-1] else - mul!(s[k], x, c[k-1]) - mul!(c[k], -x, s[k-1]) + TS.mul!(s[k], x, c[k-1]) + TS.mul!(c[k], -x, s[k-1]) end @inbounds for i = 2:k x = interval(i) * a[i] @@ -427,8 +464,8 @@ for T in (:Taylor1, :TaylorN) s[k] += x * c[k-i] c[k] -= x * s[k-i] else - mul!(s[k], x, c[k-i]) - mul!(c[k], -x, s[k-i]) + TS.mul!(s[k], x, c[k-i]) + TS.mul!(c[k], -x, s[k-i]) end end @inbounds s[k] = s[k] / interval(k) @@ -436,8 +473,8 @@ for T in (:Taylor1, :TaylorN) return nothing end - @inline function tan!(c::$T{Interval{T}}, a::$T{Interval{T}}, c2::$T{Interval{T}}, - k::Int) where {T<:Real} + @inline function TS.tan!(c::$T{Interval{T}}, a::$T{Interval{T}}, c2::$T{Interval{T}}, + k::Int) where {T<:AbstractFloat} if k == 0 @inbounds aux = tan( constant_term(a) ) @inbounds c[0] = aux @@ -447,22 +484,22 @@ for T in (:Taylor1, :TaylorN) if $T == Taylor1 @inbounds c[k] = interval(k) * a[k] * c2[0] else - @inbounds mul!(c[k], interval(k) * a[k], c2[0]) + @inbounds TS.mul!(c[k], interval(k) * a[k], c2[0]) end @inbounds for i = 1:k-1 if $T == Taylor1 c[k] += interval(k-i) * a[k-i] * c2[i] else - mul!(c[k], interval(k-i) * a[k-i], c2[i]) + TS.mul!(c[k], interval(k-i) * a[k-i], c2[i]) end end @inbounds c[k] = a[k] + c[k]/interval(k) - sqr!(c2, c, k) + TS.sqr!(c2, c, k) return nothing end - @inline function asin!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, - k::Int) where {T<:Real} + @inline function TS.asin!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, + k::Int) where {T<:AbstractFloat} if k == 0 a0 = constant_term(a) @inbounds c[0] = asin( a0 ) @@ -472,22 +509,22 @@ for T in (:Taylor1, :TaylorN) if $T == Taylor1 @inbounds c[k] = interval(k-1) * r[1] * c[k-1] else - @inbounds mul!(c[k], interval(k-1) * r[1], c[k-1]) + @inbounds TS.mul!(c[k], interval(k-1) * r[1], c[k-1]) end @inbounds for i in 2:k-1 if $T == Taylor1 c[k] += interval(k-i) * r[i] * c[k-i] else - mul!(c[k], interval(k-i) * r[i], c[k-i]) + TS.mul!(c[k], interval(k-i) * r[i], c[k-i]) end end - sqrt!(r, one(a0)-a^2, k) + TS.sqrt!(r, one(a[0])-a^2, k) @inbounds c[k] = (a[k] - c[k]/interval(k)) / constant_term(r) return nothing end - @inline function acos!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, - k::Int) where {T<:Real} + @inline function TS.acos!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, + k::Int) where {T<:AbstractFloat} if k == 0 a0 = constant_term(a) @inbounds c[0] = acos( a0 ) @@ -497,22 +534,22 @@ for T in (:Taylor1, :TaylorN) if $T == Taylor1 @inbounds c[k] = interval(k-1) * r[1] * c[k-1] else - @inbounds mul!(c[k], interval(k-1) * r[1], c[k-1]) + @inbounds TS.mul!(c[k], interval(k-1) * r[1], c[k-1]) end @inbounds for i in 2:k-1 if $T == Taylor1 c[k] += interval(k-i) * r[i] * c[k-i] else - mul!(c[k], interval(k-i) * r[i], c[k-i]) + TS.mul!(c[k], interval(k-i) * r[i], c[k-i]) end end - sqrt!(r, one(a0)-a^2, k) + TS.sqrt!(r, one(a[0])-a^2, k) @inbounds c[k] = -(a[k] + c[k]/interval(k)) / constant_term(r) return nothing end - @inline function atan!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, - k::Int) where {T<:Real} + @inline function TS.atan!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, + k::Int) where {T<:AbstractFloat} if k == 0 a0 = constant_term(a) @inbounds c[0] = atan( a0 ) @@ -522,22 +559,22 @@ for T in (:Taylor1, :TaylorN) if $T == Taylor1 @inbounds c[k] = interval(k-1) * r[1] * c[k-1] else - @inbounds mul!(c[k], interval(k-1) * r[1], c[k-1]) + @inbounds TS.mul!(c[k], interval(k-1) * r[1], c[k-1]) end @inbounds for i in 2:k-1 if $T == Taylor1 c[k] += interval(k-i) * r[i] * c[k-i] else - mul!(c[k], interval(k-i) * r[i], c[k-i]) + TS.mul!(c[k], interval(k-i) * r[i], c[k-i]) end end - @inbounds sqr!(r, a, k) + @inbounds TS.sqr!(r, a, k) @inbounds c[k] = (a[k] - c[k]/interval(k)) / constant_term(r) return nothing end - @inline function sinhcosh!(s::$T{Interval{T}}, c::$T{Interval{T}}, a::$T{Interval{T}}, - k::Int) where {T<:Real} + @inline function TS.sinhcosh!(s::$T{Interval{T}}, c::$T{Interval{T}}, a::$T{Interval{T}}, + k::Int) where {T<:AbstractFloat} if k == 0 @inbounds s[0] = sinh( constant_term(a) ) @inbounds c[0] = cosh( constant_term(a) ) @@ -548,8 +585,8 @@ for T in (:Taylor1, :TaylorN) @inbounds s[k] = x * c[k-1] @inbounds c[k] = x * s[k-1] else - @inbounds mul!(s[k], x, c[k-1]) - @inbounds mul!(c[k], x, s[k-1]) + @inbounds TS.mul!(s[k], x, c[k-1]) + @inbounds TS.mul!(c[k], x, s[k-1]) end @inbounds for i = 2:k x = interval(i) * a[i] @@ -557,8 +594,8 @@ for T in (:Taylor1, :TaylorN) s[k] += x * c[k-i] c[k] += x * s[k-i] else - mul!(s[k], x, c[k-i]) - mul!(c[k], x, s[k-i]) + TS.mul!(s[k], x, c[k-i]) + TS.mul!(c[k], x, s[k-i]) end end s[k] = s[k] / interval(k) @@ -566,8 +603,8 @@ for T in (:Taylor1, :TaylorN) return nothing end - @inline function tanh!(c::$T{Interval{T}}, a::$T{Interval{T}}, c2::$T{Interval{T}}, - k::Int) where {T<:Real} + @inline function TS.tanh!(c::$T{Interval{T}}, a::$T{Interval{T}}, c2::$T{Interval{T}}, + k::Int) where {T<:AbstractFloat} if k == 0 @inbounds aux = tanh( constant_term(a) ) @inbounds c[0] = aux @@ -577,21 +614,22 @@ for T in (:Taylor1, :TaylorN) if $T == Taylor1 @inbounds c[k] = k * a[k] * c2[0] else - @inbounds mul!(c[k], k * a[k], c2[0]) + @inbounds TS.mul!(c[k], k * a[k], c2[0]) end @inbounds for i = 1:k-1 if $T == Taylor1 c[k] += (k-i) * a[k-i] * c2[i] else - mul!(c[k], (k-i) * a[k-i], c2[i]) + TS.mul!(c[k], (k-i) * a[k-i], c2[i]) end end @inbounds c[k] = a[k] - c[k]/k - sqr!(c2, c, k) + TS.sqr!(c2, c, k) return nothing end - @inline function asinh!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, k::Int) where {T<:Real} + @inline function TS.asinh!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, + k::Int) where {T<:AbstractFloat} if k == 0 a0 = constant_term(a) @inbounds c[0] = asinh( a0 ) @@ -601,21 +639,22 @@ for T in (:Taylor1, :TaylorN) if $T == Taylor1 @inbounds c[k] = interval(k-1) * r[1] * c[k-1] else - @inbounds mul!(c[k], interval(k-1) * r[1], c[k-1]) + @inbounds TS.mul!(c[k], interval(k-1) * r[1], c[k-1]) end @inbounds for i in 2:k-1 if $T == Taylor1 c[k] += interval(k-i) * r[i] * c[k-i] else - mul!(c[k], interval(k-i) * r[i], c[k-i]) + TS.mul!(c[k], interval(k-i) * r[i], c[k-i]) end end - sqrt!(r, a^2+one(a0), k) + TS.sqrt!(r, a^2+one(a[0]), k) @inbounds c[k] = (a[k] - c[k]/interval(k)) / constant_term(r) return nothing end - @inline function acosh!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, k::Int) where {T<:Real} + @inline function TS.acosh!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, + k::Int) where {T<:AbstractFloat} if k == 0 a0 = constant_term(a) @inbounds c[0] = acosh( a0 ) @@ -625,21 +664,22 @@ for T in (:Taylor1, :TaylorN) if $T == Taylor1 @inbounds c[k] = interval(k-1) * r[1] * c[k-1] else - @inbounds mul!(c[k], interval(k-1) * r[1], c[k-1]) + @inbounds TS.mul!(c[k], interval(k-1) * r[1], c[k-1]) end @inbounds for i in 2:k-1 if $T == Taylor1 c[k] += interval(k-i) * r[i] * c[k-i] else - mul!(c[k], interval(k-i) * r[i], c[k-i]) + TS.mul!(c[k], interval(k-i) * r[i], c[k-i]) end end - sqrt!(r, a^2-one(a0), k) + TS.sqrt!(r, a^2-one(a[0]), k) @inbounds c[k] = (a[k] - c[k]/interval(k)) / constant_term(r) return nothing end - @inline function atanh!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, k::Int) where {T<:Real} + @inline function TS.atanh!(c::$T{Interval{T}}, a::$T{Interval{T}}, r::$T{Interval{T}}, + k::Int) where {T<:AbstractFloat} if k == 0 a0 = constant_term(a) @inbounds c[0] = atanh( a0 ) @@ -649,16 +689,16 @@ for T in (:Taylor1, :TaylorN) if $T == Taylor1 @inbounds c[k] = interval(k-1) * r[1] * c[k-1] else - @inbounds mul!(c[k], interval(k-1) * r[1], c[k-1]) + @inbounds TS.mul!(c[k], interval(k-1) * r[1], c[k-1]) end @inbounds for i in 2:k-1 if $T == Taylor1 c[k] += interval(k-i) * r[i] * c[k-i] else - mul!(c[k], interval(k-i) * r[i], c[k-i]) + TS.mul!(c[k], interval(k-i) * r[i], c[k-i]) end end - @inbounds sqr!(r, a, k) + @inbounds TS.sqr!(r, a, k) @inbounds c[k] = (a[k] + c[k]/interval(k)) / constant_term(r) return nothing end @@ -666,7 +706,7 @@ for T in (:Taylor1, :TaylorN) end -function evaluate(a::Taylor1, dx::Interval{S}) where {S<:Real} +function evaluate(a::Taylor1{T}, dx::Interval{S}) where {T<:Real, S<:NumTypes} order = a.order uno = one(dx) dx2 = dx^2 @@ -687,7 +727,7 @@ function evaluate(a::Taylor1, dx::Interval{S}) where {S<:Real} end function evaluate(a::TaylorN{T}, dx::AbstractVector{Interval{S}}) where - {T<:Real, S<:Real} + {T<:Real, S<:NumTypes} @assert length(dx) == get_numvars() suma = zero(constant_term(a)) + interval(zero(T)) @inbounds for homPol in reverse(eachindex(a)) @@ -718,21 +758,21 @@ end function evaluate(a::HomogeneousPolynomial{T}, - dx::AbstractVector{Interval{S}}) where {T<:Real, S<:Real} + dx::AbstractVector{Interval{S}}) where {T<:Real, S<:NumTypes} @assert length(dx) == get_numvars() all(isequal_interval.(dx, interval(-one(T), one(T)))) && - return _evaluate(a, dx, Val(true)) + return TS._evaluate(a, dx, Val(true)) all(isequal_interval.(dx, interval(zero(T), one(T)))) && - return _evaluate(a, dx, Val(false)) - return _evaluate(a, dx) + return TS._evaluate(a, dx, Val(false)) + return TS._evaluate(a, dx) end evaluate(a::TaylorN{Taylor1{T}}, vals::AbstractVector{Interval{S}}) where - {T<:Real, S<:Real} = _evaluate(a, (vals...,), Val(false)) + {T<:Real, S<:NumTypes} = TS._evaluate(a, (vals...,), Val(false)) # _evaluate -function _evaluate(a::HomogeneousPolynomial{T}, - dx::AbstractVector{Interval{S}}) where {T<:Real, S<:Real} +function TS._evaluate(a::HomogeneousPolynomial{T}, + dx::AbstractVector{Interval{S}}) where {T<:Real, S<:NumTypes} a.order == 0 && return a[1] + interval(zero(T)) ct = TS.coeff_table[a.order+1] @inbounds suma = a[1]*interval(zero(T)) @@ -744,8 +784,8 @@ function _evaluate(a::HomogeneousPolynomial{T}, return suma end -function _evaluate(a::HomogeneousPolynomial{T}, dx::AbstractVector{Interval{S}}, - ::Val{true} ) where {T<:Real, S<:Real} +function TS._evaluate(a::HomogeneousPolynomial{T}, dx::AbstractVector{Interval{S}}, + ::Val{true} ) where {T<:Real, S<:NumTypes} a.order == 0 && return a[1] + interval(zero(T)) ct = TS.coeff_table[a.order+1] @inbounds suma = a[1]*interval(zero(T)) @@ -766,8 +806,8 @@ function _evaluate(a::HomogeneousPolynomial{T}, dx::AbstractVector{Interval{S}}, return suma end -function _evaluate(a::HomogeneousPolynomial{T}, dx::AbstractVector{Interval{S}}, - ::Val{false} ) where {T<:Real, S<:Real} +function TS._evaluate(a::HomogeneousPolynomial{T}, dx::AbstractVector{Interval{S}}, + ::Val{false} ) where {T<:Real, S<:NumTypes} a.order == 0 && return a[1] + interval(zero(T)) @inbounds suma = zero(a[1])*dx[1] @inbounds for homPol in a.coeffs @@ -776,20 +816,20 @@ function _evaluate(a::HomogeneousPolynomial{T}, dx::AbstractVector{Interval{S}}, return suma end -function _evaluate(a::TaylorN{T}, vals::NTuple{N,TaylorN{Interval{S}}}) where - {N, T<:Real, S<:Real} +function TS._evaluate(a::TaylorN{T}, vals::NTuple{N,TaylorN{Interval{S}}}) where + {N, T<:Real, S<:NumTypes} @assert get_numvars() == N R = promote_type(TS.numtype(a), typeof(vals[1])) a_length = length(a) suma = zeros(R, a_length) @inbounds for homPol in 1:a_length - suma[homPol] = _evaluate(a.coeffs[homPol], vals) + suma[homPol] = TS._evaluate(a.coeffs[homPol], vals) end return suma end -function _evaluate(a::HomogeneousPolynomial{T}, - vals::NTuple{N,TaylorN{Interval{S}}}) where {N, T<:Real, S<:Real} +function TS._evaluate(a::HomogeneousPolynomial{T}, + vals::NTuple{N,TaylorN{Interval{S}}}) where {N, T<:Real, S<:NumTypes} ct = TS.coeff_table[a.order+1] suma = zero(a[1])*vals[1] for (i, a_coeff) in enumerate(a.coeffs) @@ -808,7 +848,7 @@ Normalizes `a::Taylor1` such that the interval `I` is mapped by an affine transformation to the interval `-1..1` (`symI=true`) or to `0..1` (`symI=false`). """ -normalize_taylor(a::Taylor1, I::Interval{T}, symI::Bool=true) where {T<:Real} = +normalize_taylor(a::Taylor1, I::Interval{T}, symI::Bool=true) where {T<:NumTypes} = _normalize(a, I, Val(symI)) """ @@ -819,27 +859,43 @@ are mapped by an affine transformation to the intervals `-1..1` (`symI=true`) or to `0..1` (`symI=false`). """ normalize_taylor(a::TaylorN, I::AbstractVector{Interval{T}}, - symI::Bool=true) where {T<:Real} = _normalize(a, I, Val(symI)) + symI::Bool=true) where {T<:NumTypes} = _normalize(a, I, Val(symI)) # I -> -1..1 -function _normalize(a::Taylor1, I::Interval{T}, ::Val{true}) where {T<:Real} +function _normalize(a::Taylor1, I::Interval{T}, ::Val{true}) where {T<:NumTypes} order = get_order(a) - t = Taylor1(T, order) + t = Taylor1(TS.numtype(a), order) tnew = mid(I) + t*radius(I) return a(tnew) end +function _normalize(a::Taylor1{Interval{T}}, I::Interval{T}, ::Val{true}) where + {T<:NumTypes} + order = get_order(a) + t = Taylor1(TS.numtype(a), order) + tnew = interval(mid(I)) + t*interval(radius(I)) + return a(tnew) +end + # I -> 0..1 -function _normalize(a::Taylor1, I::Interval{T}, ::Val{false}) where {T<:Real} +function _normalize(a::Taylor1, I::Interval{T}, ::Val{false}) where {T<:NumTypes} order = get_order(a) - t = Taylor1(T, order) + t = Taylor1(TS.numtype(a), order) tnew = inf(I) + t*diam(I) return a(tnew) end +function _normalize(a::Taylor1{Interval{T}}, I::Interval{T}, ::Val{false}) where + {T<:NumTypes} + order = get_order(a) + t = Taylor1(TS.numtype(a), order) + tnew = interval(inf(I)) + t*interval(diam(I)) + return a(tnew) +end + # I -> IntervalBox(-1..1, Val(N)) function _normalize(a::TaylorN, I::AbstractVector{Interval{T}}, - ::Val{true}) where {T<:Real} + ::Val{true}) where {T<:NumTypes} order = get_order(a) x = Vector{typeof(a)}(undef, length(I)) @inbounds for ind in eachindex(x) @@ -848,8 +904,20 @@ function _normalize(a::TaylorN, I::AbstractVector{Interval{T}}, return a(x) end +function _normalize(a::TaylorN{Interval{T}}, I::AbstractVector{Interval{T}}, + ::Val{true}) where {T<:NumTypes} + order = get_order(a) + x = Vector{typeof(a)}(undef, length(I)) + @inbounds for ind in eachindex(x) + x[ind] = interval(mid(I[ind])) + + TaylorN(ind, order=order)*interval(radius(I[ind])) + end + return a(x) +end + # I -> IntervalBox(0..1, Val(N)) -function _normalize(a::TaylorN, I::AbstractVector{Interval{T}}, ::Val{false}) where {T<:Real} +function _normalize(a::TaylorN, I::AbstractVector{Interval{T}}, + ::Val{false}) where {T<:NumTypes} order = get_order(a) x = Vector{typeof(a)}(undef, length(I)) @inbounds for ind in eachindex(x) @@ -858,12 +926,24 @@ function _normalize(a::TaylorN, I::AbstractVector{Interval{T}}, ::Val{false}) wh return a(x) end +function _normalize(a::TaylorN{Interval{T}}, I::AbstractVector{Interval{T}}, + ::Val{false}) where {T<:NumTypes} + order = get_order(a) + x = Vector{typeof(a)}(undef, length(I)) + @inbounds for ind in eachindex(x) + x[ind] = interval(inf(I[ind])) + + TaylorN(ind, order=order)*interval(diam(I[ind])) + end + return a(x) +end + # Printing-related methods numbr2str function TS.numbr2str(zz::Interval, ifirst::Bool=false) TS._isthinzero(zz) && return string( zz ) plusmin = ifelse( ifirst, string(""), string("+ ") ) return string(plusmin, zz) end + function TS.numbr2str(zz::ComplexI, ifirst::Bool=false) zT = zero(zz.re) TS._isthinzero(zz) && return string(zT) @@ -875,5 +955,4 @@ function TS.numbr2str(zz::ComplexI, ifirst::Bool=false) return cadena end - -end \ No newline at end of file +end diff --git a/src/power.jl b/src/power.jl index 721a43ad..8de73cf5 100644 --- a/src/power.jl +++ b/src/power.jl @@ -179,7 +179,7 @@ end # TODO: get rid of allocations function ^(a::TaylorN, r::S) where {S<:Real} a0 = constant_term(a) - aux = one(a0^r) + aux = one(a0)^r iszero(r) && return TaylorN(aux, a.order) aa = aux*a diff --git a/test/intervals.jl b/test/intervals.jl index 16e4e18f..e8b57967 100644 --- a/test/intervals.jl +++ b/test/intervals.jl @@ -40,7 +40,11 @@ using Test @test p5(x,-a) == (x-a)^5 @test p4(x,-b) == (x-b)^4 for ind in eachindex(p5(x,-b)) - @test all(issubset_interval.((p5(x,-b)[ind]).coeffs, (((x-b)^5)[ind]).coeffs)) + r1 = p5(x,-b)[ind] + r2 = ((x-b)^5)[ind] + @test all(issubset_interval.(r1.coeffs, r2.coeffs)) + @test all(isguaranteed.(getfield(r1, :coeffs))) + @test all(isguaranteed.(getfield(r2, :coeffs))) end @@ -72,7 +76,7 @@ using Test ii = interval(0, 6) t = Taylor1(4) - f(x) = 0.1 * x^3 - 0.5*x^2 + 1 + f(x) = 0.1 * x^3 - 0.5 * x^2 + 1 ft = f(t) f1 = normalize_taylor(ft, ii, true) f2 = normalize_taylor(ft, ii, false) @@ -92,13 +96,25 @@ using Test ii = c t = Taylor1(5) g(x) = 1 - x^4 + x^5 - gt = g(t) - g1 = normalize_taylor(gt, c, true) + gt1 = g(t) + gn1 = normalize_taylor(gt1, c, true) @test issubset_interval(interval(g(4/5), 1), g(ii)) + @test issubset_interval(interval(g(4/5), 1), gt1(ii)) + @test issubset_interval(interval(g(4/5), 1), gn1(b)) + @test isinterior(gn1(b), g(ii)) + @test diam(gn1(b)) < diam(gt1(ii)) + # + ti = Taylor1(typeof(c), 5) + gg(x) = interval(1) - x^4 + x^5 + gt = gg(ti) + @test all(isguaranteed.(getfield(gt, :coeffs))) + gn2 = normalize_taylor(gt, c, true) + @test all(isguaranteed.(getfield(gn2, :coeffs))) + @test issubset_interval(interval(g(4/5), 1), gg(ii)) @test issubset_interval(interval(g(4/5), 1), gt(ii)) - @test issubset_interval(interval(g(4/5), 1), g1(b)) - @test isinterior(g1(b), g(ii)) - @test diam(g1(b)) < diam(gt(ii)) + @test issubset_interval(interval(g(4/5), 1), gn2(b)) + @test isinterior(gn2(b), g(ii)) + @test diam(gn2(b)) < diam(gt(ii)) # Test display for Taylor1{Complex{Interval{T}}} From b8f8625ebd2219baee3b5c0ca602e14791582fb9 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Fri, 9 Feb 2024 11:47:09 -0600 Subject: [PATCH 09/13] Require Julia 1.9 at least (IA compatibility) --- .github/workflows/ci.yml | 2 +- Project.toml | 2 +- src/TaylorSeries.jl | 6 +----- src/arithmetic.jl | 4 ++-- 4 files changed, 5 insertions(+), 9 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index dc3f3b11..c446a847 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -33,7 +33,7 @@ jobs: strategy: fail-fast: false matrix: - julia-version: ['1.6', '1', 'nightly'] + julia-version: ['1.9', '1', 'nightly'] julia-arch: [x64] os: [ubuntu-latest, macOS-latest, windows-latest] # # 32-bit Julia binaries are not available on macOS diff --git a/Project.toml b/Project.toml index fc84b7ab..033dfaba 100644 --- a/Project.toml +++ b/Project.toml @@ -32,7 +32,7 @@ Requires = "0.5.2, 1" SparseArrays = "<0.0.1, 1" StaticArrays = "1" Test = "<0.0.1, 1" -julia = "1" +julia = "1.9" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl index 2e508e66..7b1b02ef 100644 --- a/src/TaylorSeries.jl +++ b/src/TaylorSeries.jl @@ -26,11 +26,7 @@ end using LinearAlgebra: norm, mul!, lu, lu!, LinearAlgebra.lutype, LinearAlgebra.copy_oftype, - LinearAlgebra.issuccess - -if VERSION >= v"1.7.0-DEV.1188" - using LinearAlgebra: NoPivot, RowMaximum -end + LinearAlgebra.issuccess, NoPivot, RowMaximum import LinearAlgebra: norm, mul!, lu diff --git a/src/arithmetic.jl b/src/arithmetic.jl index e44a856a..988e014f 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -1262,8 +1262,8 @@ end # end # see https://github.com/JuliaLang/julia/pull/40623 -const LU_RowMaximum = VERSION >= v"1.7.0-DEV.1188" ? RowMaximum() : Val(true) -const LU_NoPivot = VERSION >= v"1.7.0-DEV.1188" ? NoPivot() : Val(false) +const LU_RowMaximum = RowMaximum() +const LU_NoPivot = NoPivot() # Adapted from (Julia v1.2) stdlib/v1.2/LinearAlgebra/src/lu.jl#240-253 # and (Julia v1.4.0-dev) stdlib/LinearAlgebra/v1.4/src/lu.jl#270-274, From ffb4ebabe3fecb325b775b48dc1e3eea24032c4f Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Sun, 10 Mar 2024 10:45:32 -0600 Subject: [PATCH 10/13] Fix docs --- docs/src/api.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/src/api.md b/docs/src/api.md index c33907b4..1ce7e81e 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -93,6 +93,7 @@ asinh! acosh! atanh! differentiate! +_isthinzero _internalmutfunc_call _dict_unary_ops _dict_binary_calls From c7d80199b08f475e19fa23f67df07be31dabeee5 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Fri, 15 Mar 2024 15:06:31 -0600 Subject: [PATCH 11/13] Fix test --- test/intervals.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/intervals.jl b/test/intervals.jl index e8b57967..2b960a26 100644 --- a/test/intervals.jl +++ b/test/intervals.jl @@ -137,10 +137,10 @@ using Test # Iss 351 (inspired by a test in ReachabilityAnalysis) - p1 = Taylor1([0 .. 0, (0 .. 0.1) + (0 .. 0.01) * y], 4) - p2 = Taylor1([0 .. 0, (0 .. 0.5) + (0 .. 0.02) * x + (0 .. 0.03) * y], 4) - @test evaluate([p1, p2], 0 .. 1) == [p1[1], p2[1]] - @test typeof(p1(0 .. 1)) == TaylorN{Interval{Float64}} + p1 = Taylor1([interval(0, 0), interval(0, 0.1) + interval(0, 0.01) * y], 4) + p2 = Taylor1([interval(0, 0), interval(0, 0.5) + interval(0, 0.02) * x + interval(0, 0.03) * y], 4) + @test evaluate([p1, p2], interval(0, 1)) == [p1[1], p2[1]] + @test typeof(p1(interval(0, 1))) == TaylorN{Interval{Float64}} # Tests related to Iss #311 # `sqrt` and `pow` defined on Interval(0,Inf) From 7d3ca435a1f4a94431fef03e76a7559afd7bd607 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Sat, 25 May 2024 10:29:18 -0600 Subject: [PATCH 12/13] Fix tests again, and clean-up --- ext/TaylorSeriesIAExt.jl | 1 - test/intervals.jl | 6 ++++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/ext/TaylorSeriesIAExt.jl b/ext/TaylorSeriesIAExt.jl index 02b13fc8..d2f59bf8 100644 --- a/ext/TaylorSeriesIAExt.jl +++ b/ext/TaylorSeriesIAExt.jl @@ -80,7 +80,6 @@ function ^(a::Taylor1{Interval{T}}, r::S) where {T<:NumTypes, S<:Real} # aa = one(aux) * a aa = convert(Taylor1{typeof(aux)}, a) aa[0] = one(aux) * a0 - @show(aa) r == 1 && return aa r == 2 && return TS.square(aa) r == 1/2 && return sqrt(aa) diff --git a/test/intervals.jl b/test/intervals.jl index 2b960a26..8625e158 100644 --- a/test/intervals.jl +++ b/test/intervals.jl @@ -6,6 +6,8 @@ using TaylorSeries, IntervalArithmetic using Test # eeuler = Base.MathConstants.e +setdisplay(:full) + @testset "Tests Taylor1 and TaylorN expansions over Intervals" begin a = interval(1, 2) b = interval(-1, 1) @@ -175,8 +177,8 @@ using Test @test_throws DomainError acosh(interval(0.0, 1.0) + x) @test acosh(interval(0.0, 2.0) + x) == acosh(interval(1.0, 2.0) + x) # atanh defined on interval(-1,1) - @test_logs (:warn, "invalid interval, empty interval is returned") @test_throws DomainError atanh(interval(1.0, 1.0) + ti) + @test_logs (:warn, "ill-formed bare interval [a, b] with a = Inf, b = Inf. Empty interval is returned") @test_throws DomainError atanh(interval(1.0, 1.0) + ti) @test atanh(interval(-2.0, 0.0) + ti) == atanh(interval(-1.0, 0.0) + ti) - @test_logs (:warn, "invalid interval, empty interval is returned") @test_throws DomainError atanh(interval(1.0, 1.0) + y) + @test_logs (:warn, "ill-formed bare interval [a, b] with a = Inf, b = Inf. Empty interval is returned") @test_throws DomainError atanh(interval(1.0, 1.0) + y) @test atanh(interval(-2.0, 0.0) + y) == atanh(interval(-1.0, 0.0) + y) end From d2a70fcfe759b0b2bdb2ed155cc420c334cdb3e5 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Wed, 31 Jul 2024 15:18:55 -0600 Subject: [PATCH 13/13] Fixes related to specialized methods pow!, sqr! and accsqr! --- ext/TaylorSeriesIAExt.jl | 52 +++++++++++++++++++++++++++++++++------- src/power.jl | 9 ++++--- test/intervals.jl | 10 ++++---- 3 files changed, 52 insertions(+), 19 deletions(-) diff --git a/ext/TaylorSeriesIAExt.jl b/ext/TaylorSeriesIAExt.jl index d2f59bf8..5d5cf2d5 100644 --- a/ext/TaylorSeriesIAExt.jl +++ b/ext/TaylorSeriesIAExt.jl @@ -91,7 +91,7 @@ function ^(a::Taylor1{Interval{T}}, r::S) where {T<:NumTypes, S<:Real} c_order = l0 == 0 ? a.order : min(a.order, trunc(Int,r*a.order)) c = Taylor1(zero(aux), c_order) for k in eachindex(c) - TS.pow!(c, aa, r, k) + TS.pow!(c, aa, c, r, k) end return c end @@ -124,30 +124,64 @@ for T = (:Taylor1, :TaylorN) @inline function TS.sqr!(c::$T{Interval{T}}, a::$T{Interval{T}}, k::Int) where {T<:NumTypes} if k == 0 - @inbounds c[0] = constant_term(a)^2 + TS.sqr_orderzero!(c, a) return nothing end + # Sanity + TS.zero!(c, k) + # Recursion formula kodd = k%2 kend = div(k - 2 + kodd, 2) - @inbounds for i = 0:kend - if $T == Taylor1 + if $T == Taylor1 + @inbounds for i = 0:kend c[k] += a[i] * a[k-i] - else + end + else + @inbounds for i = 0:kend TS.mul!(c[k], a[i], a[k-i]) end end - @inbounds c[k] = interval(2) * c[k] + @inbounds TS.mul!(c, interval(2), c, k) kodd == 1 && return nothing if $T == Taylor1 - @inbounds c[k] += a[div(k,2)]^2 + @inbounds c[k] += a[k >> 1]^2 else - TS.sqr!(c[k], a[div(k,2)]) + TS.accsqr!(c[k], a[k >> 1]) end return nothing end + + @inline function TS.sqr!(c::$T{Interval{T}}, k::Int) where {T<:NumTypes} + if k == 0 + TS.sqr_orderzero!(c, c) + return nothing + end + # Recursion formula + kodd = k%2 + kend = (k - 2 + kodd) >> 1 + if $T == Taylor1 + (kend ≥ 0) && ( @inbounds c[k] = c[0] * c[k] ) + @inbounds for i = 1:kend + c[k] += c[i] * c[k-i] + end + @inbounds c[k] = interval(2) * c[k] + (kodd == 0) && ( @inbounds c[k] += c[k >> 1]^2 ) + else + (kend ≥ 0) && ( @inbounds TS.mul!(c, c[0][1], c, k) ) + @inbounds for i = 1:kend + TS.mul!(c[k], c[i], c[k-i]) + end + @inbounds TS.mul!(c, interval(2), c, k) + if (kodd == 0) + TS.accsqr!(c[k], c[k >> 1]) + end + end + + return nothing + end end end -@inline function TS.sqr!(c::HomogeneousPolynomial{Interval{T}}, +@inline function TS.accsqr!(c::HomogeneousPolynomial{Interval{T}}, a::HomogeneousPolynomial{Interval{T}}) where {T<:NumTypes} iszero(a) && return nothing @inbounds num_coeffs_a = TS.size_table[a.order+1] diff --git a/src/power.jl b/src/power.jl index 8de73cf5..22cbc937 100644 --- a/src/power.jl +++ b/src/power.jl @@ -231,10 +231,10 @@ end # Homogeneous coefficients for real power @doc doc""" - pow!(c, a, r::Real, k::Int) + pow!(c, a, aux, r::Real, k::Int) Update the `k`-th expansion coefficient `c[k]` of `c = a^r`, for -both `c` and `a` either `Taylor1` or `TaylorN`. +both `c`, `a` and `aux` either `Taylor1` or `TaylorN`. The coefficients are given by @@ -381,9 +381,8 @@ end # The recursion formula for i = 0:ordT-lnull-1 ((i+lnull) > a.order || (l0+kprime-i > a.order)) && continue - aux = r*(kprime-i) - i - # res[ordT] += aux*res[i+lnull]*a[l0+kprime-i] - @inbounds mul_scalar!(res[ordT], aux, res[i+lnull], a[l0+kprime-i]) + aaux = r*(kprime-i) - i + @inbounds mul_scalar!(res[ordT], aaux, res[i+lnull], a[l0+kprime-i]) end # res[ordT] /= a[l0]*kprime @inbounds div_scalar!(res[ordT], 1/kprime, a[l0]) diff --git a/test/intervals.jl b/test/intervals.jl index 8625e158..612aa6c5 100644 --- a/test/intervals.jl +++ b/test/intervals.jl @@ -41,12 +41,12 @@ setdisplay(:full) @test p4(x,-a) == (x-a)^4 @test p5(x,-a) == (x-a)^5 @test p4(x,-b) == (x-b)^4 + r1 = p5(x,-b) + r2 = (x-b)^5 for ind in eachindex(p5(x,-b)) - r1 = p5(x,-b)[ind] - r2 = ((x-b)^5)[ind] - @test all(issubset_interval.(r1.coeffs, r2.coeffs)) - @test all(isguaranteed.(getfield(r1, :coeffs))) - @test all(isguaranteed.(getfield(r2, :coeffs))) + @test all(issubset_interval.(r1[ind].coeffs, r2[ind].coeffs)) + @test all(isguaranteed.(getfield(r1[ind], :coeffs))) + @test all(isguaranteed.(getfield(r2[ind], :coeffs))) end