diff --git a/Project.toml b/Project.toml index 69e45edb..f2c4149d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "TaylorSeries" uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git" -version = "0.14.0" +version = "0.15.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/docs/src/api.md b/docs/src/api.md index 20a8c4db..f8a0ae90 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -50,6 +50,7 @@ inverse abs norm isapprox +isless isfinite displayBigO use_show_default diff --git a/docs/src/userguide.md b/docs/src/userguide.md index 4a88411b..30446a07 100644 --- a/docs/src/userguide.md +++ b/docs/src/userguide.md @@ -105,6 +105,16 @@ of order 5; this is be consistent with the number of known coefficients of the returned series, since the result corresponds to factorize `t` in the numerator to cancel the same factor in the denominator. +`Taylor1` is also equiped with a total order, provided by overloading [`isless`](@ref). +The ordering is consistent with the interpretation that there are infinitessimal +elements in the algebra; for details see M. Berz, "Automatic Differentiation as +Nonarchimedean Analysis", Computer Arithmetic and Enclosure Methods, (1992), Elsevier, +439-450. This is illustrated by: + +```@repl userguide +1 > t > 2*t^2 > 100*t^3 >= 0 +``` + If no valid Taylor expansion can be computed an error is thrown, for instance, when a derivative is not defined at the expansion point, or it simply diverges. @@ -336,12 +346,29 @@ Note that some of the arithmetic operations have been extended for [`HomogeneousPolynomial`](@ref); division, for instance, is not extended. The same convention used for `Taylor1` objects is used when combining `TaylorN` polynomials of different order. +Both `HomogeneousPolynomial` and `TaylorN` are equiped with a total *lexicographical* +order, provided by overloading [`isless`](@ref). +The ordering is consistent with the interpretation that there are infinitessimal +elements in the algebra; for details see M. Berz, "Automatic Differentiation as +Nonarchimedean Analysis", Computer Arithmetic and Enclosure Methods, (1992), Elsevier, +439-450. +Essentially, the lexicographic order works as follows: smaller order monomials +are *larger* than higher order monomials; when the order is the same, *larger* monomials +appear before in the hash-tables; the function [`show_monomials`](@ref). +```@repl userguide +x, y = set_variables("x y", order=10); + +show_monomials(2) +``` +Then, the following then holds: +```@repl userguide +0 < 1e8 * y^2 < x*y < x^2 < y < x/1e8 < 1.0 +``` The elementary functions have also been implemented, again by computing their coefficients recursively: ```@repl userguide -x, y = set_variables("x y", order=10); exy = exp(x+y) ``` diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl index 4703edd5..2ebf3292 100644 --- a/src/TaylorSeries.jl +++ b/src/TaylorSeries.jl @@ -39,7 +39,7 @@ import Base: ==, +, -, *, /, ^ import Base: iterate, size, eachindex, firstindex, lastindex, eltype, length, getindex, setindex!, axes, copyto! -import Base: zero, one, zeros, ones, isinf, isnan, iszero, +import Base: zero, one, zeros, ones, isinf, isnan, iszero, isless, convert, promote_rule, promote, show, real, imag, conj, adjoint, rem, mod, mod2pi, abs, abs2, diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 6b83f480..521eabf1 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -10,7 +10,6 @@ ## Equality ## for T in (:Taylor1, :TaylorN) - @eval begin ==(a::$T{T}, b::$T{S}) where {T<:Number,S<:Number} = ==(promote(a,b)...) @@ -34,11 +33,127 @@ function ==(a::HomogeneousPolynomial, b::HomogeneousPolynomial) return iszero(a.coeffs) && iszero(b.coeffs) end +## Total ordering ## +for T in (:Taylor1, :TaylorN) + @eval begin + @inline function isless(a::$T{<:Number}, b::Real) + a0 = constant_term(a) + a0 != b && return isless(a0, b) + nz = findfirst(a-b) + if nz == -1 + return isless(zero(a0), zero(b)) + else + return isless(a[nz], zero(b)) + end + end + @inline function isless(b::Real, a::$T{<:Number}) + a0 = constant_term(a) + a0 != b && return isless(b, a0) + nz = findfirst(b-a) + if nz == -1 + return isless(zero(b), zero(a0)) + else + return isless(zero(b), a[nz]) + end + end + # + @inline isless(a::$T{T}, b::$T{S}) where {T<:Number, S<:Number} = isless(promote(a,b)...) + @inline isless(a::$T{T}, b::$T{T}) where {T<:Number} = isless(a - b, zero(constant_term(a))) + end +end + +@inline function isless(a::HomogeneousPolynomial{<:Number}, b::Real) + orda = get_order(a) + if orda == 0 + return isless(a[1], b) + else + !iszero(b) && return isless(zero(a[1]), b) + nz = max(findfirst(a), 1) + return isless(a[nz], b) + end +end +@inline function isless(b::Real, a::HomogeneousPolynomial{<:Number}) + orda = get_order(a) + if orda == 0 + return isless(b, a[1]) + else + !iszero(b) && return isless(b, zero(a[1])) + nz = max(findfirst(a),1) + return isless(b, a[nz]) + end +end +# +@inline isless(a::HomogeneousPolynomial{T}, b::HomogeneousPolynomial{S}) where {T<:Number, S<:Number} = isless(promote(a,b)...) +@inline function isless(a::HomogeneousPolynomial{T}, b::HomogeneousPolynomial{T}) where {T<:Number} + orda = get_order(a) + ordb = get_order(b) + if orda == ordb + return isless(a-b, zero(a[1])) + elseif orda < ordb + return isless(a, zero(a[1])) + else + return isless(-b, zero(a[1])) + end +end + +# Mixtures +@inline isless(a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries} = + isless(a - b, zero(T)) +@inline function isless(a::HomogeneousPolynomial{Taylor1{T}}, b::HomogeneousPolynomial{Taylor1{T}}) where {T<:NumberNotSeries} + orda = get_order(a) + ordb = get_order(b) + if orda == ordb + return isless(a-b, zero(T)) + elseif orda < ordb + return isless(a, zero(T)) + else + return isless(-b, zero(T)) + end +end +@inline isless(a::TaylorN{Taylor1{T}}, b::TaylorN{Taylor1{T}}) where {T<:NumberNotSeries} = isless(a - b, zero(T)) + +#= TODO: Nested Taylor1s; needs careful thinking; iss #326. The following works: +@inline isless(a::Taylor1{Taylor1{T}}, b::Taylor1{Taylor1{T}}) where {T<:Number} = isless(a - b, zero(T)) +# Is the following correct? +# ti = Taylor1(3) +# to = Taylor1([zero(ti), one(ti)], 9) +# tito = ti * to +# ti > to > 0 # ok +# to^2 < toti < ti^2 # ok +# ti > ti^2 > to # is this ok? +=# + +@doc doc""" + isless(a::Taylor1{<:Real}, b::Real) + isless(a::TaylorN{<:Real}, b::Real) + +Compute `isless` by comparing the `constant_term(a)` and `b`. If they are equal, +returns `a[nz] < 0`, with `nz` the first +non-zero coefficient after the constant term. This defines a total order. + +For many variables, the ordering includes a lexicographical convention in order to be +total. We have opted for the simplest one, where the *larger* variable appears *before* +when the `TaylorN` variables are defined (e.g., through [`set_variables`](@ref)). + +Refs: +- M. Berz, AIP Conference Proceedings 177, 275 (1988); https://doi.org/10.1063/1.37800 +- M. Berz, "Automatic Differentiation as Nonarchimedean Analysis", Computer Arithmetic and + Enclosure Methods, (1992), Elsevier, 439-450. + +--- + + isless(a::Taylor1{<:Real}, b::Taylor1{<:Real}) + isless(a::TaylorN{<:Real}, b::Taylor1{<:Real}) + +Return `isless(a - b, zero(b))`. +""" isless + + +## zero and one ## for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @eval iszero(a::$T) = iszero(a.coeffs) end -## zero and one ## for T in (:Taylor1, :TaylorN) @eval zero(a::$T) = $T(zero.(a.coeffs)) @eval function one(a::$T) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 3ff92d73..f628ac73 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -268,11 +268,35 @@ 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) - isa(first, Nothing) && return -1 + 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 + +# 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 diff --git a/test/manyvariables.jl b/test/manyvariables.jl index 9082890c..e637d8fa 100644 --- a/test/manyvariables.jl +++ b/test/manyvariables.jl @@ -93,17 +93,23 @@ end @test v == [1,2,0] @test_throws AssertionError TS.resize_coeffsHP!(v,1) hpol_v = HomogeneousPolynomial(v) + @test findfirst(hpol_v) == 1 + @test findlast(hpol_v) == 2 hpol_v[3] = 3 @test v == [1,2,3] hpol_v[1:3] = 3 @test v == [3,3,3] hpol_v[1:2:2] = 0 @test v == [0,3,3] + @test findfirst(hpol_v) == 2 + @test findlast(hpol_v) == 3 hpol_v[1:1:2] = [1,2] @test all(hpol_v[1:1:2] .== [1,2]) @test v == [1,2,3] hpol_v[:] = zeros(Int, 3) @test hpol_v == 0 + @test findfirst(hpol_v) == -1 + @test findlast(hpol_v) == -1 tn_v = TaylorN(HomogeneousPolynomial(zeros(Int, 3))) tn_v[0] = 1 @@ -132,9 +138,12 @@ end @test (@inferred real(xH)) == xH xT = TaylorN(xH, 17) yT = TaylorN(Int, 2, order=17) + @test findfirst(xT) == 1 + @test findlast(yT) == 1 @test (@inferred conj(xT)) == (@inferred adjoint(xT)) @test (@inferred real(xT)) == (xT) zeroT = zero( TaylorN([xH],1) ) + @test findfirst(zeroT) == findlast(zeroT) == -1 @test (@inferred imag(xT)) == (zeroT) @test zeroT.coeffs == zeros(HomogeneousPolynomial{Int}, 1) @test size(xH) == (2,) @@ -156,6 +165,16 @@ end @test TS.fixorder(xH,yH) == (xH,yH) @test_throws AssertionError TS.fixorder(zeros(xH,0)[1],yH) + @testset "Lexicographic order" begin + @test HomogeneousPolynomial([2.1]) > 0.5 * xH > yH > xH^2 + @test -1.0*xH < -yH^2 < 0 < HomogeneousPolynomial([1.0]) + @test -3 < HomogeneousPolynomial([-2]) < 0.0 + @test !(zero(yH^2) > 0) + @test 1 > 0.5 * xT > yT > xT^2 > 0 > TaylorN(-1.0, 3) + @test -1 < -0.25 * yT < -xT^2 < 0.0 < TaylorN(1, 3) + @test !(zero(xT) > 0) + @test !(zero(yT^2) < 0) + end @test constant_term(xT) == 0 @test constant_term(uT) == 1.0 @@ -761,6 +780,9 @@ end @test yT[0] == xT[0]*(pi/180) TS.rad2deg!(yT, xT, 0) @test yT[0] == xT[0]*(180/pi) + # Lexicographic tests with 4 vars + @test 1 > dx[1] > dx[2] > dx[3] > dx[4] + @test dx[4]^2 < dx[3]*dx[4] < dx[3]^2 < dx[2]*dx[4] < dx[2]*dx[3] < dx[2]^2 < dx[1]*dx[4] < dx[1]*dx[3] < dx[1]*dx[2] < dx[1]^2 end @testset "Integrate for several variables" begin diff --git a/test/mixtures.jl b/test/mixtures.jl index 891db4e5..3df76031 100644 --- a/test/mixtures.jl +++ b/test/mixtures.jl @@ -14,6 +14,7 @@ using LinearAlgebra xH = HomogeneousPolynomial(Int, 1) yH = HomogeneousPolynomial(Int, 2) tN = Taylor1(TaylorN{Float64}, 3) + @test findfirst(tN) == 1 @test convert(eltype(tN), tN) == tN @test eltype(xH) == HomogeneousPolynomial{Int} @@ -22,6 +23,13 @@ using LinearAlgebra @test TS.numtype(tN) == TaylorN{Float64} @test normalize_taylor(tN) == tN @test tN.order == 3 + + @testset "Lexicographic order: Taylor1{HomogeneousPolynomial{T}} and Taylor1{TaylorN{T}}" begin + @test HomogeneousPolynomial([1]) > xH > yH > 0.0 + @test -xH^2 < -xH*yH < -yH^2 < -xH^3 < -yH^3 < HomogeneousPolynomial([0.0]) + @test 1 ≥ tN > 2*tN^2 > 100*tN^3 > 0 + @test -2*tN < -tN^2 ≤ 0 + end @test string(zero(tN)) == " 0.0 + 𝒪(‖x‖¹) + 𝒪(t⁴)" @test string(tN) == " ( 1.0 + 𝒪(‖x‖¹)) t + 𝒪(t⁴)" @test string(tN + 3Taylor1(Int, 2)) == " ( 4.0 + 𝒪(‖x‖¹)) t + 𝒪(t³)" @@ -43,6 +51,7 @@ using LinearAlgebra t = Taylor1(3) xHt = HomogeneousPolynomial(typeof(t), 1) + @test findfirst(xHt) == 1 @test convert(eltype(xHt), xHt) === xHt @test eltype(xHt) == HomogeneousPolynomial{Taylor1{Float64}} @test TS.numtype(xHt) == Taylor1{Float64} @@ -50,6 +59,7 @@ using LinearAlgebra @test string(xHt) == " ( 1.0 + 𝒪(t¹)) x₁" xHt = HomogeneousPolynomial([one(t), zero(t)]) yHt = HomogeneousPolynomial([zero(t), t]) + @test findfirst(yHt) == 2 @test string(xHt) == " ( 1.0 + 𝒪(t⁴)) x₁" @test string(yHt) == " ( 1.0 t + 𝒪(t⁴)) x₂" @test string(HomogeneousPolynomial([t])) == " ( 1.0 t + 𝒪(t⁴))" @@ -62,10 +72,13 @@ using LinearAlgebra @test (xHt+yHt)([1, 1]) == (xHt+yHt)((1, 1)) tN1 = TaylorN([HomogeneousPolynomial([t]), xHt, yHt^2]) + @test findfirst(tN1) == 0 @test tN1[0] == HomogeneousPolynomial([t]) @test tN1(t,one(t)) == 2t+t^2 + @test findfirst(tN1(t,one(t))) == 1 @test tN1([t,one(t)]) == tN1((t,one(t))) t1N = convert(Taylor1{TaylorN{Float64}}, tN1) + @test findfirst(zero(tN1)) == -1 @test t1N[0] == HomogeneousPolynomial(1) ctN1 = convert(TaylorN{Taylor1{Float64}}, t1N) @test convert(eltype(tN1), tN1) === tN1 @@ -103,6 +116,13 @@ using LinearAlgebra @test !(@inferred isnan(tN1)) @test !(@inferred isinf(tN1)) + @testset "Lexicographic order: HomogeneousPolynomial{Taylor1{T}} and TaylorN{Taylor1{T}}" begin + @test 1 > xHt > yHt > xHt^2 > 0 + @test -xHt^2 < -xHt*yHt < -yHt^2 < -xHt^3 < -yHt^3 < 0 + @test 1 ≥ tN1 > 2*tN1^2 > 100*tN1^3 > 0 + @test -2*tN1 < -tN1^2 ≤ 0 + end + @test mod(tN1+1,1.0) == 0+tN1 @test mod(tN1-1.125,2) == 0.875+tN1 @test (rem(tN1+1.125,1.0))[0][1] == 0.125 + t @@ -332,6 +352,7 @@ end @testset "Tests with nested Taylor1s" begin ti = Taylor1(3) to = Taylor1([zero(ti), one(ti)], 9) + @test findfirst(to) == 1 @test TS.numtype(to) == Taylor1{Float64} @test normalize_taylor(to) == to @test normalize_taylor(Taylor1([zero(to), one(to)], 5)) == Taylor1([zero(to), one(to)], 5) @@ -340,12 +361,16 @@ end @test string(to^2) == " ( 1.0 + 𝒪(t⁴)) t² + 𝒪(t¹⁰)" @test ti + to == Taylor1([ti, one(ti)], 9) tito = ti * to + # The next tests are related to iss #326 + # @test ti > ti^2 > to > 0 + # @test to^2 < toti < ti^2 @test tito == Taylor1([zero(ti), ti], 9) @test tito / to == ti @test get_order(tito/to) == get_order(to)-1 @test tito / ti == to @test get_order(tito/ti) == get_order(to) @test ti^2-to^2 == (ti+to)*(ti-to) + @test findfirst(ti^2-to^2) == 0 @test sin(to) ≈ Taylor1(one(ti) .* sin(Taylor1(10)).coeffs, 9) @test to(1 + ti) == 1 + ti @test to(1 + ti) isa Taylor1{Float64} diff --git a/test/onevariable.jl b/test/onevariable.jl index 71d908b4..33dfa29f 100644 --- a/test/onevariable.jl +++ b/test/onevariable.jl @@ -45,6 +45,11 @@ Base.iszero(::SymbNumber) = false @test axes(t) == () @test axes([t]) == (Base.OneTo(1),) + @testset "Total order" begin + @test 1 + t ≥ 1.0 > 0.5 + t > t^2 ≥ zero(t) + @test -1.0 < -1/1000 - t < -t < -t^2 ≤ 0 + end + v = [1,2] @test typeof(TaylorSeries.resize_coeffs1!(v,3)) == Nothing @test v == [1,2,0,0] @@ -138,14 +143,10 @@ Base.iszero(::SymbNumber) = false @test length.( TaylorSeries.fixorder(zt, Taylor1([1], 1)) ) == (2, 2) @test eltype(TaylorSeries.fixorder(zt,Taylor1([1]))[1]) == Taylor1{Int} @test TS.numtype(TaylorSeries.fixorder(zt,Taylor1([1]))[1]) == Int - @test TaylorSeries.findfirst(t) == 1 - @test TaylorSeries.findfirst(t^2) == 2 - @test TaylorSeries.findfirst(ot) == 0 - @test TaylorSeries.findfirst(zt) == -1 - @test TaylorSeries.findlast(t) == 1 - @test TaylorSeries.findlast(t^2) == 2 - @test TaylorSeries.findlast(ot) == 0 - @test TaylorSeries.findlast(zt) == -1 + @test findfirst(t) == 1 + @test findfirst(t^2) == 2 + @test findfirst(ot) == 0 + @test findfirst(zt) == -1 @test iszero(zero(t)) @test !iszero(one(t)) @test @inferred isinf(Taylor1([typemax(1.0)]))