From f5587a772e99d93428c5201f5422d5629fa4fcaf Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Wed, 31 Jan 2024 13:10:31 -0600 Subject: [PATCH] Improve `Taylor1{TaylorN{T}}` mixtures (#333) * Add tests for Taylor1{TaylorN{T}} * Add fixorder to Taylor1{TaylorN{T}} * Add methods for arithm opers for Taylor1{TaylorN{T}} * Add methods of ^, sqrt, pow!, sqr!, sqrt! for Taylor1{TaylorN{T} * Add specialized methods for most mutating functions for Taylor1{TaylorN{T}} * Implement suggestions by @PerezHz * Fix error in docs * Use TS * Improvements for evaluate * Clean up commented (old) code * Remove some superfluous methods in evaluate * More fixes * Fixes and clean-up of `evaluate` methods * Returning a zero Homogeneouspolynomial is of order 0 * Scattered fixes and clean up * Fix documentation * Remove evaluate methods for TaylorN with Varargs, and fix tests * Fix _evaluate for mixtures, and new tests * Methods of evaluate for Taylor1 involving TaylorN (fixing certain type instability) * Add tests for new evaluate methods * Bump minor version * Remove a `@show` left there --- Project.toml | 2 +- docs/src/api.md | 3 +- ext/TaylorSeriesIAExt.jl | 18 +- src/arithmetic.jl | 273 ++++++++++++++------ src/auxiliary.jl | 34 ++- src/calculus.jl | 44 ++-- src/conversion.jl | 32 +-- src/dictmutfunct.jl | 18 +- src/evaluate.jl | 310 +++++++++++----------- src/functions.jl | 539 +++++++++++++++++++++++++++++++++++++-- src/hash_tables.jl | 4 +- src/other_functions.jl | 16 +- src/power.jl | 309 ++++++++++++++++++---- test/manyvariables.jl | 1 - test/mixtures.jl | 140 ++++++++-- test/mutatingfuncts.jl | 68 ++--- test/onevariable.jl | 108 ++++---- test/runtests.jl | 2 - 18 files changed, 1421 insertions(+), 500 deletions(-) diff --git a/Project.toml b/Project.toml index 24b3803f..b781e891 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.15.4" +version = "0.16.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/docs/src/api.md b/docs/src/api.md index d3cd70e6..ae190e60 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -77,7 +77,7 @@ div! pow! square sqr! -sqr!(::HomogeneousPolynomial{T}, ::HomogeneousPolynomial{T}) where {T<:NumberNotSeriesN} +accsqr! sqrt! exp! log! @@ -105,6 +105,5 @@ _populate_dicts! ```@index Pages = ["api.md"] -Module = ["TaylorSeries"] Order = [:type, :function] ``` diff --git a/ext/TaylorSeriesIAExt.jl b/ext/TaylorSeriesIAExt.jl index 51c9bc3f..61e4b946 100644 --- a/ext/TaylorSeriesIAExt.jl +++ b/ext/TaylorSeriesIAExt.jl @@ -41,7 +41,7 @@ function ^(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real} 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 - TaylorSeries.pow!(c, aa, r, k) + TS.pow!(c, aa, r, k) end return c @@ -66,7 +66,7 @@ function ^(a::TaylorN{Interval{T}}, r::S) where {T<:Real, S<:Real} c = TaylorN( a0r, a.order) for ord in 1:a.order - TaylorSeries.pow!(c, aa, r, ord) + TS.pow!(c, aa, r, ord) end return c @@ -84,7 +84,7 @@ for T in (:Taylor1, :TaylorN) aa[0] = one(aux) * a0 c = $T( aux, order ) for k in eachindex(a) - TaylorSeries.log!(c, aa, k) + TS.log!(c, aa, k) end return c end @@ -101,7 +101,7 @@ for T in (:Taylor1, :TaylorN) c = $T( aux, order ) r = $T( sqrt(1 - a0^2), order ) for k in eachindex(a) - TaylorSeries.asin!(c, aa, r, k) + TS.asin!(c, aa, r, k) end return c end @@ -118,7 +118,7 @@ for T in (:Taylor1, :TaylorN) c = $T( aux, order ) r = $T( sqrt(1 - a0^2), order ) for k in eachindex(a) - TaylorSeries.acos!(c, aa, r, k) + TS.acos!(c, aa, r, k) end return c end @@ -135,7 +135,7 @@ for T in (:Taylor1, :TaylorN) c = $T( aux, order ) r = $T( sqrt(a0^2 - 1), order ) for k in eachindex(a) - TaylorSeries.acosh!(c, aa, r, k) + TS.acosh!(c, aa, r, k) end return c end @@ -152,14 +152,14 @@ for T in (:Taylor1, :TaylorN) """Series expansion of atanh(x) diverges at x = ±1.""")) for k in eachindex(a) - TaylorSeries.atanh!(c, aa, r, k) + TS.atanh!(c, aa, r, k) end return c end end -function evaluate(a::Taylor1, dx::Interval) +function evaluate(a::Taylor1{T}, dx::Interval{S}) where {T<:Real, S<:Real} order = a.order uno = one(dx) dx2 = dx^2 @@ -201,7 +201,7 @@ 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) - ct = TaylorSeries.coeff_table[a.order+1] + ct = TS.coeff_table[a.order+1] @inbounds suma = a[1]*Interval{T}(0,0) Ieven = Interval{T}(0,1) diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 521eabf1..04b0e2f6 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -22,17 +22,18 @@ for T in (:Taylor1, :TaylorN) end end -function ==(a::Taylor1{TaylorN{T}}, b::TaylorN{Taylor1{S}}) where {T,S} - R = promote_type(S,T) - return a == convert(Taylor1{TaylorN{R}},b) +function ==(a::Taylor1{TaylorN{T}}, b::TaylorN{Taylor1{S}}) where {T, S} + R = promote_type(T, S) + return a == convert(Taylor1{TaylorN{R}}, b) end -==(b::TaylorN{Taylor1{S}}, a::Taylor1{TaylorN{T}}) where {T,S} = a == b +==(b::TaylorN{Taylor1{S}}, a::Taylor1{TaylorN{T}}) where {T, S} = a == b function ==(a::HomogeneousPolynomial, b::HomogeneousPolynomial) a.order == b.order && return a.coeffs == b.coeffs return iszero(a.coeffs) && iszero(b.coeffs) end + ## Total ordering ## for T in (:Taylor1, :TaylorN) @eval begin @@ -57,8 +58,10 @@ for T in (:Taylor1, :TaylorN) 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))) + @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 @@ -83,8 +86,10 @@ end 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} +@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 @@ -97,9 +102,10 @@ 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} +@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 @@ -110,7 +116,8 @@ end 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)) +@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)) @@ -145,7 +152,7 @@ Refs: isless(a::Taylor1{<:Real}, b::Taylor1{<:Real}) isless(a::TaylorN{<:Real}, b::Taylor1{<:Real}) -Return `isless(a - b, zero(b))`. +Returns `isless(a - b, zero(b))`. """ isless @@ -163,10 +170,8 @@ for T in (:Taylor1, :TaylorN) end end -function zero(a::HomogeneousPolynomial{T}) where {T<:Number} - v = zero.(a.coeffs) - return HomogeneousPolynomial(v, a.order) -end +zero(a::HomogeneousPolynomial{T}) where {T<:Number} = + HomogeneousPolynomial(zero.(a.coeffs), a.order) function zeros(a::HomogeneousPolynomial{T}, order::Int) where {T<:Number} order == 0 && return [HomogeneousPolynomial([zero(a[1])], 0)] @@ -205,35 +210,39 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!))) for T in (:Taylor1, :TaylorN) @eval begin - ($f)(a::$T{T}, b::$T{S}) where {T<:Number,S<:Number} = - $f(promote(a,b)...) + ($f)(a::$T{T}, b::$T{S}) where {T<:Number, S<:Number} = + ($f)(promote(a, b)...) - function $f(a::$T{T}, b::$T{T}) where {T<:Number} + function ($f)(a::$T{T}, b::$T{T}) where {T<:Number} if a.order != b.order a, b = fixorder(a, b) end - v = similar(a.coeffs) - @__dot__ v = $f(a.coeffs, b.coeffs) - return $T(v, a.order) + c = $T( zero(constant_term(a)), a.order) + for k in eachindex(a) + ($fc)(c, a, b, k) + end + return c end - function $f(a::$T) - v = similar(a.coeffs) - @__dot__ v = $f(a.coeffs) - return $T(v, a.order) + function ($f)(a::$T) + c = $T( zero(constant_term(a)), a.order) + for k in eachindex(a) + ($fc)(c, a, k) + end + return c end - ($f)(a::$T{T}, b::S) where {T<:Number, S<:Number} = $f(promote(a, b)...) + ($f)(a::$T{T}, b::S) where {T<:Number, S<:Number} = ($f)(promote(a, b)...) - function $f(a::$T{T}, b::T) where {T<:Number} + function ($f)(a::$T{T}, b::T) where {T<:Number} coeffs = copy(a.coeffs) @inbounds coeffs[1] = $f(a[0], b) return $T(coeffs, a.order) - end + end # ($f)(b::S, a::$T{T}) where {T<:Number,S<:Number} = $f(promote(b, a)...) - function $f(b::T, a::$T{T}) where {T<:Number} + function ($f)(b::T, a::$T{T}) where {T<:Number} coeffs = similar(a.coeffs) @__dot__ coeffs = ($f)(a.coeffs) @inbounds coeffs[1] = $f(b, a[0]) @@ -245,49 +254,62 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!))) @inbounds v[k] = ($f)(a[k]) return nothing end + function ($fc)(v::$T{T}, a::T, k::Int) where {T<:Number} @inbounds v[k] = k==0 ? ($f)(a) : zero(a) return nothing end + function ($fc)(v::$T, a::$T, b::$T, k::Int) @inbounds v[k] = ($f)(a[k], b[k]) return nothing end + function ($fc)(v::$T, a::$T, b::Number, k::Int) - @inbounds v[k] = k==0 ? ($f)(a[0], b) : ($f)(a[k], zero(b)) + @inbounds v[k] = k==0 ? + ($f)(constant_term(a), b) : ($f)(a[k], zero(b)) return nothing end + function ($fc)(v::$T, a::Number, b::$T, k::Int) - @inbounds v[k] = k==0 ? ($f)(a, b[0]) : ($f)(zero(a), b[k]) + @inbounds v[k] = k==0 ? + ($f)(a, constant_term(b)) : ($f)(zero(a), b[k]) return nothing end end end - @eval $(f)(a::T, b::S) where {T<:Taylor1, S<:TaylorN} = $(f)(promote(a, b)...) - @eval $(f)(a::T, b::S) where {T<:TaylorN, S<:Taylor1} = $(f)(promote(a, b)...) + @eval ($f)(a::T, b::S) where {T<:Taylor1, S<:TaylorN} = ($f)(promote(a, b)...) + @eval ($f)(a::T, b::S) where {T<:TaylorN, S<:Taylor1} = ($f)(promote(a, b)...) @eval begin ($f)(a::HomogeneousPolynomial{T}, b::HomogeneousPolynomial{S}) where - {T<:NumberNotSeriesN,S<:NumberNotSeriesN} = $f(promote(a,b)...) + {T<:NumberNotSeriesN,S<:NumberNotSeriesN} = ($f)(promote(a,b)...) - function $f(a::HomogeneousPolynomial{T}, b::HomogeneousPolynomial{T}) where + function ($f)(a::HomogeneousPolynomial{T}, b::HomogeneousPolynomial{T}) where {T<:NumberNotSeriesN} @assert a.order == b.order v = similar(a.coeffs) - @__dot__ v = $f(a.coeffs, b.coeffs) + @__dot__ v = ($f)(a.coeffs, b.coeffs) return HomogeneousPolynomial(v, a.order) end - function $f(a::HomogeneousPolynomial) + # NOTE add! and subst! for HomogeneousPolynomial's act as += or -= + function ($fc)(res::HomogeneousPolynomial{T}, a::HomogeneousPolynomial{T}, + b::HomogeneousPolynomial{T}, k::Int) where {T<:NumberNotSeriesN} + res[k] += ($f)(a[k], b[k]) + return nothing + end + + function ($f)(a::HomogeneousPolynomial) v = similar(a.coeffs) - @__dot__ v = $f(a.coeffs) + @__dot__ v = ($f)(a.coeffs) return HomogeneousPolynomial(v, a.order) end function ($f)(a::TaylorN{Taylor1{T}}, b::Taylor1{S}) where - {T<:NumberNotSeries,S<:NumberNotSeries} + {T<:NumberNotSeries, S<:NumberNotSeries} @inbounds aux = $f(a[0][1], b) R = TS.numtype(aux) coeffs = Array{HomogeneousPolynomial{Taylor1{R}}}(undef, a.order+1) @@ -297,7 +319,7 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!))) end function ($f)(b::Taylor1{S}, a::TaylorN{Taylor1{T}}) where - {T<:NumberNotSeries,S<:NumberNotSeries} + {T<:NumberNotSeries, S<:NumberNotSeries} @inbounds aux = $f(b, a[0][1]) R = TS.numtype(aux) coeffs = Array{HomogeneousPolynomial{Taylor1{R}}}(undef, a.order+1) @@ -309,22 +331,23 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!))) function ($f)(a::Taylor1{TaylorN{T}}, b::TaylorN{S}) where {T<:NumberNotSeries,S<:NumberNotSeries} @inbounds aux = $f(a[0], b) - R = TS.numtype(aux) - coeffs = Array{TaylorN{R}}(undef, a.order+1) - coeffs .= a.coeffs - @inbounds coeffs[1] = aux - return Taylor1(coeffs, a.order) - end + c = Taylor1( zero(aux), a.order) + for k in eachindex(a) + ($fc)(c, a, b, k) + end + return c + end function ($f)(b::TaylorN{S}, a::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries,S<:NumberNotSeries} @inbounds aux = $f(b, a[0]) - R = TS.numtype(aux) - coeffs = Array{TaylorN{R}}(undef, a.order+1) - @__dot__ coeffs = $f(a.coeffs) - @inbounds coeffs[1] = aux - return Taylor1(coeffs, a.order) + c = Taylor1( zero(aux), a.order) + for k in eachindex(a) + ($fc)(c, a, b, k) + end + return c end + end end @@ -333,14 +356,14 @@ end for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @eval begin - function *(a::T, b::$T{S}) where {T<:NumberNotSeries,S<:NumberNotSeries} + function *(a::T, b::$T{S}) where {T<:NumberNotSeries, S<:NumberNotSeries} @inbounds aux = a * b.coeffs[1] v = Array{typeof(aux)}(undef, length(b.coeffs)) @__dot__ v = a * b.coeffs return $T(v, b.order) end - *(b::$T{S}, a::T) where {T<:NumberNotSeries,S<:NumberNotSeries} = a * b + *(b::$T{S}, a::T) where {T<:NumberNotSeries, S<:NumberNotSeries} = a * b function *(a::T, b::$T{T}) where {T<:Number} v = Array{T}(undef, length(b.coeffs)) @@ -356,8 +379,7 @@ for T in (:HomogeneousPolynomial, :TaylorN) @eval begin function *(a::Taylor1{T}, b::$T{Taylor1{S}}) where - {T<:NumberNotSeries,S<:NumberNotSeries} - + {T<:NumberNotSeries, S<:NumberNotSeries} @inbounds aux = a * b.coeffs[1] R = typeof(aux) coeffs = Array{R}(undef, length(b.coeffs)) @@ -366,9 +388,9 @@ for T in (:HomogeneousPolynomial, :TaylorN) end *(b::$T{Taylor1{R}}, a::Taylor1{T}) where - {T<:NumberNotSeries,R<:NumberNotSeries} = a * b + {T<:NumberNotSeries, R<:NumberNotSeries} = a * b - function *(a::$T{T}, b::Taylor1{$T{S}}) where {T<:NumberNotSeries,S<:NumberNotSeries} + function *(a::$T{T}, b::Taylor1{$T{S}}) where {T<:NumberNotSeries, S<:NumberNotSeries} @inbounds aux = a * b[0] R = typeof(aux) coeffs = Array{R}(undef, length(b.coeffs)) @@ -376,7 +398,7 @@ for T in (:HomogeneousPolynomial, :TaylorN) return Taylor1(coeffs, b.order) end - *(b::Taylor1{$T{S}}, a::$T{T}) where {T<:NumberNotSeries,S<:NumberNotSeries} = a * b + *(b::Taylor1{$T{S}}, a::$T{T}) where {T<:NumberNotSeries, S<:NumberNotSeries} = a * b end end @@ -385,7 +407,7 @@ for (T, W) in ((:Taylor1, :Number), (:TaylorN, :NumberNotSeriesN)) if a.order != b.order a, b = fixorder(a, b) end - c = $T(zero(a[0]), a.order) + c = $T(zero(constant_term(a)), a.order) for ord in eachindex(c) mul!(c, a, b, ord) # updates c[ord] end @@ -402,18 +424,36 @@ function *(a::HomogeneousPolynomial{T}, b::HomogeneousPolynomial{T}) where order = a.order + b.order - order > get_order() && return HomogeneousPolynomial(zero(a[1]), get_order()) + # NOTE: the following returns order 0, but could be get_order(), or get_order(a) + order > get_order() && return HomogeneousPolynomial(zero(a[1]), get_order(a)) res = HomogeneousPolynomial(zero(a[1]), order) mul!(res, a, b) return res end +function *(a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{S}}) where + {T<:NumberNotSeries, S<:NumberNotSeries} + R = promote_type(T,S) + return *(convert(Taylor1{TaylorN{R}}, a), convert(Taylor1{TaylorN{R}}, b)) +end + +function *(a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries} + if (a.order != b.order) || any(get_order.(a.coeffs) .!= get_order.(b.coeffs)) + a, b = fixorder(a, b) + end + res = Taylor1(zero(a[0]), a.order) + for ordT in eachindex(a) + mul!(res, a, b, ordT) + end + return res +end + # Internal multiplication functions for T in (:Taylor1, :TaylorN) + # NOTE: For $T = TaylorN, `mul!` *accumulates* the result of a * b in c[k] @eval @inline function mul!(c::$T{T}, a::$T{T}, b::$T{T}, k::Int) where {T<:Number} - if $T == Taylor1 @inbounds c[k] = a[0] * b[k] else @@ -426,7 +466,6 @@ for T in (:Taylor1, :TaylorN) mul!(c[k], a[i], b[k-i]) end end - return nothing end @@ -440,11 +479,25 @@ for T in (:Taylor1, :TaylorN) end end +function mul!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}}, + ordT::Int) where {T<:NumberNotSeries} + # Sanity + zero!(res, a, ordT) + for k in 0:ordT + @inbounds for ordQ in eachindex(a[ordT]) + mul!(res[ordT], a[k], b[ordT-k], ordQ) + end + end + return nothing +end + + @doc doc""" mul!(c, a, b, k::Int) --> nothing Update the `k`-th expansion coefficient `c[k]` of `c = a * b`, where all `c`, `a`, and `b` are either `Taylor1` or `TaylorN`. +Note that for `TaylorN` the result of `a * b` is accumulated in `c[k]`. The coefficients are given by @@ -458,7 +511,8 @@ c_k = \sum_{j=0}^k a_j b_{k-j}. """ mul!(c, a, b) --> nothing -Return `c = a*b` with no allocation; all arguments are `HomogeneousPolynomial`. +Accumulates in `c[k]` the result of `a*b` with minimum allocation; all arguments +are `HomogeneousPolynomial`. """ @inline function mul!(c::HomogeneousPolynomial, a::HomogeneousPolynomial, @@ -495,7 +549,7 @@ end ## Division ## -function /(a::Taylor1{Rational{T}}, b::S) where {T<:Integer,S<:NumberNotSeries} +function /(a::Taylor1{Rational{T}}, b::S) where {T<:Integer, S<:NumberNotSeries} R = typeof( a[0] // b) v = Array{R}(undef, a.order+1) @__dot__ v = a.coeffs // b @@ -503,7 +557,7 @@ function /(a::Taylor1{Rational{T}}, b::S) where {T<:Integer,S<:NumberNotSeries} end for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) - @eval function /(a::$T{T}, b::S) where {T<:NumberNotSeries,S<:NumberNotSeries} + @eval function /(a::$T{T}, b::S) where {T<:NumberNotSeries, S<:NumberNotSeries} @inbounds aux = a.coeffs[1] / b v = Array{typeof(aux)}(undef, length(a.coeffs)) @__dot__ v = a.coeffs / b @@ -512,15 +566,20 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @eval function /(a::$T{T}, b::T) where {T<:Number} @inbounds aux = a.coeffs[1] / b - v = Array{typeof(aux)}(undef, length(a.coeffs)) - @__dot__ v = a.coeffs / b - return $T(v, a.order) + # v = Array{typeof(aux)}(undef, length(a.coeffs)) + # @__dot__ v = a.coeffs / b + # return $T(v, a.order) + c = $T( zero(aux), a.order ) + for ord in eachindex(c) + div!(c, a, b, ord) # updates c[ord] + end + return c end end for T in (:HomogeneousPolynomial, :TaylorN) @eval function /(b::$T{Taylor1{S}}, a::Taylor1{T}) where - {T<:NumberNotSeries,S<:NumberNotSeries} + {T<:NumberNotSeries, S<:NumberNotSeries} @inbounds aux = b.coeffs[1] / a R = typeof(aux) coeffs = Array{R}(undef, length(b.coeffs)) @@ -528,7 +587,7 @@ for T in (:HomogeneousPolynomial, :TaylorN) return $T(coeffs, b.order) end - @eval function /(b::$T{Taylor1{T}}, a::S) where {T<:NumberNotSeries,S<:NumberNotSeries} + @eval function /(b::$T{Taylor1{T}}, a::S) where {T<:NumberNotSeries, S<:NumberNotSeries} @inbounds aux = b.coeffs[1] / a R = typeof(aux) coeffs = Array{R}(undef, length(b.coeffs)) @@ -537,16 +596,21 @@ for T in (:HomogeneousPolynomial, :TaylorN) end @eval function /(b::Taylor1{$T{S}}, a::$T{T}) where - {T<:NumberNotSeries,S<:NumberNotSeries} + {T<:NumberNotSeries, S<:NumberNotSeries} @inbounds aux = b[0] / a - R = typeof(aux) - coeffs = Array{R}(undef, length(b.coeffs)) - @__dot__ coeffs = b.coeffs / a - return Taylor1(coeffs, b.order) + # R = typeof(aux) + # coeffs = Array{R}(undef, length(b.coeffs)) + # @__dot__ coeffs = b.coeffs / a + # return Taylor1(coeffs, b.order) + v = Taylor1(zero(aux), b.order) + @inbounds for k in eachindex(b) + v[k] = b[k] / a + end + return v end end -/(a::Taylor1{T}, b::Taylor1{S}) where {T<:Number,S<:Number} = /(promote(a,b)...) +/(a::Taylor1{T}, b::Taylor1{S}) where {T<:Number, S<:Number} = /(promote(a,b)...) function /(a::Taylor1{T}, b::Taylor1{T}) where {T<:Number} iszero(a) && !iszero(b) && return zero(a) @@ -566,7 +630,7 @@ function /(a::Taylor1{T}, b::Taylor1{T}) where {T<:Number} end /(a::TaylorN{T}, b::TaylorN{S}) where - {T<:NumberNotSeriesN,S<:NumberNotSeriesN} = /(promote(a,b)...) + {T<:NumberNotSeriesN, S<:NumberNotSeriesN} = /(promote(a,b)...) function /(a::TaylorN{T}, b::TaylorN{T}) where {T<:NumberNotSeriesN} @assert !iszero(constant_term(b)) @@ -585,6 +649,22 @@ function /(a::TaylorN{T}, b::TaylorN{T}) where {T<:NumberNotSeriesN} return c end +function /(a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries} + iszero(a) && !iszero(b) && return zero(a) + if (a.order != b.order) || any(get_order.(a.coeffs) .!= get_order.(b.coeffs)) + a, b = fixorder(a, b) + end + + # order and coefficient of first factorized term + ordfact, cdivfact = divfactorization(a, b) + + res = Taylor1(cdivfact, a.order-ordfact) + for ordT in eachindex(res) + div!(res, a, b, ordT) + end + return res +end + @inline function divfactorization(a1::Taylor1, b1::Taylor1) # order of first factorized term; a1 and b1 assumed to be of the same order @@ -654,6 +734,7 @@ end div!(v::Taylor1, b::NumberNotSeries, a::Taylor1, k::Int) = div!(v::Taylor1, Taylor1(b, a.order), a, k) +# NOTE: Here `div!` *accumulates* the result of a / b in c[k] (k > 0) @inline function div!(c::TaylorN, a::TaylorN, b::TaylorN, k::Int) if k==0 @inbounds c[0] = a[0] / constant_term(b) @@ -667,6 +748,40 @@ div!(v::Taylor1, b::NumberNotSeries, a::Taylor1, k::Int) = return nothing end +# NOTE: Here `div!` *accumulates* the result of a / b in res[k] (k > 0) +@inline function div!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, + b::Taylor1{TaylorN{T}}, ordT::Int) where {T<:NumberNotSeriesN} + + # order and coefficient of first factorized term + ordfact, cdivfact = divfactorization(a, b) + if ordT == 0 + @inbounds res[0] = cdivfact + return nothing + end + zero!(res, a, ordT) + imin = max(0, ordT+ordfact-b.order) + aux = TaylorN(zero(a[ordT][0][1]), a[ordT].order) + for k in imin:ordT-1 + @inbounds for ordQ in eachindex(a[ordT]) + mul!(aux, res[k], b[ordT+ordfact-k], ordQ) + end + end + + if ordT+ordfact ≤ b.order + for ordQ in eachindex(a[ordT]) + subst!(aux, a[ordT+ordfact], aux, ordQ) + end + else + for ordQ in eachindex(a[ordT]) + subst!(aux, aux, ordQ) + end + end + for ordQ in eachindex(res[ordT]) + div!(res[ordT], aux, b[ordfact], ordQ) + end + + return nothing +end diff --git a/src/auxiliary.jl b/src/auxiliary.jl index f628ac73..fd0406f5 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -245,15 +245,22 @@ Returns the type of the elements of the coefficients of `a`. @inline normalize_taylor(a::AbstractSeries) = a +## _minorder +function _minorder(a, b) + minorder, maxorder = minmax(a.order, b.order) + if minorder ≤ 0 + minorder = maxorder + end + return minorder +end + + ## fixorder ## for T in (:Taylor1, :TaylorN) @eval begin - @inline function fixorder(a::$T, b::$T) + @inline function fixorder(a::$T{T}, b::$T{T}) where {T<:Number} a.order == b.order && return a, b - minorder, maxorder = minmax(a.order, b.order) - if minorder ≤ 0 - minorder = maxorder - end + minorder = _minorder(a, b) return $T(copy(a.coeffs), minorder), $T(copy(b.coeffs), minorder) end end @@ -264,6 +271,23 @@ function fixorder(a::HomogeneousPolynomial, b::HomogeneousPolynomial) return a, b end +for T in (:HomogeneousPolynomial, :TaylorN) + @eval function fixorder(a::Taylor1{$T{T}}, b::Taylor1{$T{S}}) where + {T<:NumberNotSeries, S<:NumberNotSeries} + (a.order == b.order) && (all(get_order.(a.coeffs) .== get_order.(b.coeffs))) && return a, b + minordT = _minorder(a, b) + aa = Taylor1(copy(a.coeffs), minordT) + bb = Taylor1(copy(b.coeffs), minordT) + for ind in eachindex(aa) + aa[ind].order == bb[ind].order && continue + minordQ = _minorder(aa[ind], bb[ind]) + aa[ind] = TaylorN(aa[ind].coeffs, minordQ) + bb[ind] = TaylorN(bb[ind].coeffs, minordQ) + end + return aa, bb + end +end + # Finds the first non zero entry; extended to Taylor1 function Base.findfirst(a::Taylor1{T}) where {T<:Number} diff --git a/src/calculus.jl b/src/calculus.jl index 30db4239..d3460230 100644 --- a/src/calculus.jl +++ b/src/calculus.jl @@ -41,7 +41,7 @@ function differentiate!(res::Taylor1, a::Taylor1) for ord in eachindex(res) differentiate!(res, a, ord) end - nothing + return nothing end """ @@ -58,9 +58,14 @@ p_k = (k+1) a_{k+1}. """ function differentiate!(p::Taylor1, a::Taylor1, k::Int) - if k < a.order - @inbounds p[k] = (k+1)*a[k+1] - end + k >= a.order && return nothing + @inbounds p[k] = (k+1)*a[k+1] + return nothing +end +function differentiate!(p::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, k::Int) where + {T<:NumberNotSeries} + k >= a.order && return nothing + @inbounds p[k] = (k+1)*a[k+1] return nothing end @@ -73,16 +78,17 @@ Compute recursively the `Taylor1` polynomial of the n-th derivative of """ function differentiate(a::Taylor1{T}, n::Int) where {T <: Number} if n > a.order - return Taylor1(T, 0) + return Taylor1(zero(T), 0) + elseif n == a.order + return Taylor1(differentiate(n, a), 0) elseif n==0 return a - else - res = differentiate(a) - for i = 2:n - differentiate!(res, res) - end - return Taylor1(view(res.coeffs, 1:a.order-n+1)) end + res = differentiate(a) + for i = 2:n + differentiate!(res, res) + end + return Taylor1(res.coeffs[1:a.order-n+1]) end """ @@ -92,22 +98,24 @@ Return the value of the `n`-th differentiate of the polynomial `a`. """ function differentiate(n::Int, a::Taylor1{T}) where {T<:Number} @assert a.order ≥ n ≥ 0 - factorial( widen(n) ) * a[n] :: T + return factorial( widen(n) ) * a[n] :: T end + ## Integrating ## """ integrate(a, [x]) Return the integral of `a::Taylor1`. The constant of integration (0-th order coefficient) is set to `x`, which is zero if omitted. +Note that the order of the result is `a.order+1`. """ -function integrate(a::Taylor1{T}, x::S) where {T<:Number,S<:Number} +function integrate(a::Taylor1{T}, x::S) where {T<:Number, S<:Number} order = get_order(a) aa = a[0]/1 + zero(x) R = typeof(aa) coeffs = Array{typeof(aa)}(undef, order+1) - fill!(coeffs, zero(aa)) + # fill!(coeffs, zero(aa)) @inbounds for i = 1:order coeffs[i+1] = a[i-1] / i end @@ -130,7 +138,7 @@ function differentiate(a::HomogeneousPolynomial, r::Int) T = TS.numtype(a) a.order == 0 && return HomogeneousPolynomial(zero(a[1]), 0) @inbounds num_coeffs = size_table[a.order] - coeffs = zeros(T,num_coeffs) + coeffs = zeros(T, num_coeffs) @inbounds posTb = pos_table[a.order] @inbounds num_coeffs = size_table[a.order+1] @@ -184,7 +192,7 @@ function differentiate(a::TaylorN, ntup::NTuple{N,Int}) where {N} aa = copy(a) for nvar in 1:get_numvars() - for numder in 1:ntup[nvar] + for _ in 1:ntup[nvar] aa = differentiate(aa, nvar) end end @@ -230,7 +238,7 @@ function gradient(f::TaylorN) end return grad end -const ∇ = TaylorSeries.gradient +const ∇ = TS.gradient """ ``` @@ -296,7 +304,6 @@ function jacobian!(jac::Array{T,2}, vf::Array{TaylorN{T},1}) where {T<:Number} end function jacobian!(jac::Array{T,2}, vf::Array{TaylorN{T},1}, vals::Array{T,1}) where {T<:Number} - numVars = get_numvars() @assert length(vf) == numVars == length(vals) @assert (numVars, numVars) == size(jac) @@ -355,6 +362,7 @@ function integrate(a::HomogeneousPolynomial, r::Int) @assert 1 ≤ r ≤ get_numvars() order_max = get_order() + # NOTE: the following returns order 0, but could be get_order(), or get_order(a) a.order == order_max && return HomogeneousPolynomial(zero(a[1]/1), 0) @inbounds posTb = pos_table[a.order+2] diff --git a/src/conversion.jl b/src/conversion.jl index 59e9f8dc..154f5f17 100644 --- a/src/conversion.jl +++ b/src/conversion.jl @@ -86,7 +86,7 @@ convert(::Type{TaylorN}, b::Number) = TaylorN( [HomogeneousPolynomial([b], 0)], function convert(::Type{TaylorN{Taylor1{T}}}, s::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries} - orderN = get_order() + orderN = maximum(get_order.(s[:])) r = zeros(HomogeneousPolynomial{Taylor1{T}}, orderN) v = zeros(T, s.order+1) @@ -127,28 +127,6 @@ function convert(::Type{Taylor1{TaylorN{T}}}, s::TaylorN{Taylor1{T}}) where {T<: return Taylor1(vT) end -# iss 339: this triggers some invalidations -# function convert(::Type{Array{TaylorN{Taylor1{T}},N}}, -# s::Array{Taylor1{TaylorN{T}},N}) where {T<:NumberNotSeries,N} -# -# v = Array{TaylorN{Taylor1{T}}}(undef, size(s)) -# @simd for ind in eachindex(s) -# @inbounds v[ind] = convert(TaylorN{Taylor1{T}}, s[ind]) -# end -# return v -# end -# -# function convert(::Type{Array{Taylor1{TaylorN{T}}}}, -# s::Array{TaylorN{Taylor1{T}}}) where {T<:NumberNotSeries} -# -# v = Array{Taylor1{TaylorN{T}}}(undef, size(s)) -# @simd for ind in eachindex(s) -# @inbounds v[ind] = convert(Taylor1{TaylorN{T}}, s[ind]) -# end -# return v -# end - - # Promotion promote_rule(::Type{Taylor1{T}}, ::Type{Taylor1{T}}) where {T<:Number} = Taylor1{T} @@ -198,14 +176,6 @@ promote_rule(::Type{TaylorN{T}}, ::Type{S}) where {T<:Number,S<:Number} = TaylorN{promote_type(T,S)} -# iss 339: this triggers some invalidations -# # Order may matter -# promote_rule(::Type{S}, ::Type{T}) where {S<:NumberNotSeries, T<:AbstractSeries} = -# promote_rule(T,S) -# -# # disambiguation with Base -# promote_rule(::Type{Bool}, ::Type{T}) where {T<:AbstractSeries} = promote_rule(T, Bool) - promote_rule(::Type{S}, ::Type{T}) where {S<:AbstractIrrational,T<:AbstractSeries} = promote_rule(T,S) diff --git a/src/dictmutfunct.jl b/src/dictmutfunct.jl index efaf4f93..d5ec11cf 100644 --- a/src/dictmutfunct.jl +++ b/src/dictmutfunct.jl @@ -21,11 +21,11 @@ Contains parameters and expressions that allow a simple programmatic construction for calling the internal mutating functions. """ -struct _InternalMutFuncs - namef :: Symbol # internal name of the function - argsf :: NTuple # arguments - defexpr :: Expr # defining expr - auxexpr :: Expr # auxiliary expr +struct _InternalMutFuncs{N} + namef :: Symbol # internal name of the function + argsf :: NTuple{N,Symbol} # arguments + defexpr :: Expr # defining expr + auxexpr :: Expr # auxiliary expr end # Constructor @@ -183,13 +183,6 @@ internal mutating functions. Evaluating the entries generates expressions that represent the actual calls to the internal mutating functions. """ _dict_unary_calls -# const _dict_unary_calls = Dict{Symbol, NTuple{3,Expr}}() - -# #Populates the constant vector `_dict_unary_calls`. -# for kk in keys(_dict_unary_ops) -# res = _internalmutfunc_call( _InternalMutFuncs(_dict_unary_ops[kk]) ) -# push!(_dict_unary_calls, kk => res ) -# end @doc """ `_dict_binary_calls::Dict{Symbol, NTuple{2,Expr}}` @@ -203,5 +196,4 @@ internal mutating functions. Evaluating the entries generates symbols that represent the actual calls to the internal mutating functions. """ _dict_binary_calls -# const _dict_binary_calls = Dict{Symbol, NTuple{3,Expr}}() diff --git a/src/evaluate.jl b/src/evaluate.jl index 9d117571..0cca2983 100644 --- a/src/evaluate.jl +++ b/src/evaluate.jl @@ -15,18 +15,18 @@ omitted, its value is considered as zero. Note that the syntax `a(dx)` is equivalent to `evaluate(a,dx)`, and `a()` is equivalent to `evaluate(a)`. """ function evaluate(a::Taylor1{T}, dx::T) where {T<:Number} - @inbounds suma = a[end] - @inbounds for k in a.order-1:-1:0 + @inbounds suma = zero(a[end]) + @inbounds for k in reverse(eachindex(a)) suma = suma*dx + a[k] end - suma + return suma end -function evaluate(a::Taylor1{T}, dx::S) where {T<:Number,S<:Number} - suma = a[end]*one(dx) - @inbounds for k in a.order-1:-1:0 +function evaluate(a::Taylor1{T}, dx::S) where {T<:Number, S<:Number} + suma = a[end]*zero(dx) + @inbounds for k in reverse(eachindex(a)) suma = suma*dx + a[k] end - suma + return suma end evaluate(a::Taylor1{T}) where {T<:Number} = a[0] @@ -40,111 +40,86 @@ syntax `x(δt)` is equivalent to `evaluate(x, δt)`, and `x()` is equivalent to `evaluate(x)`. """ evaluate(x::AbstractArray{Taylor1{T}}, δt::S) where - {T<:Number,S<:Number} = evaluate.(x, δt) + {T<:Number, S<:Number} = evaluate.(x, δt) evaluate(a::AbstractArray{Taylor1{T}}) where {T<:Number} = getcoeff.(a, 0) """ - evaluate!(x, δt, x0) - -Evaluates each element of `x::AbstractArray{Taylor1{T}}`, -representing the Taylor expansion for the dependent variables -of an ODE at *time* `δt`. It updates the vector `x0` with the -computed values. -""" -function evaluate!(x::AbstractArray{Taylor1{T}}, δt::S, - x0::AbstractArray{T}) where {T<:Number,S<:Number} - x0 .= evaluate.( x, δt ) - nothing -end - - -""" - evaluate(a, x) + evaluate(a::Taylor1, x::Taylor1) Substitute `x::Taylor1` as independent variable in a `a::Taylor1` polynomial. Note that the syntax `a(x)` is equivalent to `evaluate(a, x)`. """ -evaluate(a::Taylor1{T}, x::Taylor1{S}) where {T<:Number,S<:Number} = - evaluate(promote(a,x)...) +evaluate(a::Taylor1{T}, x::Taylor1{S}) where {T<:Number, S<:Number} = + evaluate(promote(a, x)...) function evaluate(a::Taylor1{T}, x::Taylor1{T}) where {T<:Number} if a.order != x.order a, x = fixorder(a, x) end - @inbounds suma = a[end]*one(x) - @inbounds for k = a.order-1:-1:0 + @inbounds suma = a[end]*zero(x) + @inbounds for k in reverse(eachindex(a)) suma = suma*x + a[k] end - suma + return suma end -function evaluate(a::Taylor1{Taylor1{T}}, x::Taylor1{T}) where {T<:Number} - @inbounds suma = a[end]*one(x) - @inbounds for k = a.order-1:-1:0 +function evaluate(a::Taylor1{Taylor1{T}}, x::Taylor1{T}) where {T<:NumberNotSeriesN} + @inbounds suma = a[end]*zero(x) + @inbounds for k in reverse(eachindex(a)) suma = suma*x + a[k] end - suma + return suma end -function evaluate(a::Taylor1{T}, x::Taylor1{Taylor1{T}}) where {T<:Number} - @inbounds suma = a[end]*one(x) - @inbounds for k = a.order-1:-1:0 +function evaluate(a::Taylor1{T}, x::Taylor1{Taylor1{T}}) where {T<:NumberNotSeriesN} + @inbounds suma = a[end]*zero(x) + @inbounds for k in reverse(eachindex(a)) suma = suma*x + a[k] end - suma -end - -evaluate(p::Taylor1{T}, x::AbstractArray{S}) where {T<:Number,S<:Number} = evaluate.(Ref(p), x) - - -#function-like behavior for Taylor1 -(p::Taylor1)(x) = evaluate(p, x) -(p::Taylor1)() = evaluate(p) - -#function-like behavior for Vector{Taylor1} -if VERSION >= v"1.3" - (p::AbstractArray{Taylor1{T}})(x) where {T<:Number} = evaluate.(p, x) - (p::AbstractArray{Taylor1{T}})() where {T<:Number} = evaluate.(p) -else - (p::Array{Taylor1{T}})(x) where {T<:Number} = evaluate.(p, x) - (p::SubArray{Taylor1{T}})(x) where {T<:Number} = evaluate.(p, x) - (p::Array{Taylor1{T}})() where {T<:Number} = evaluate.(p) - (p::SubArray{Taylor1{T}})() where {T<:Number} = evaluate.(p) + return suma end - -## In place evaluation of multivariable arrays -function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{T,1}, - x0::AbstractArray{T}) where {T<:Number} - @inbounds for i in eachindex(x, x0) - x0[i] = evaluate( x[i], δx ) +evaluate(p::Taylor1{T}, x::AbstractArray{S}) where {T<:Number, S<:Number} = + evaluate.(Ref(p), x) + +function evaluate(a::Taylor1{TaylorN{T}}, dx::S) where + {T<:NumberNotSeries, S<:NumberNotSeries} + @inbounds suma = TaylorN( zero(T)*constant_term(dx), a[0].order ) + @inbounds for k in reverse(eachindex(a)) + for ordQ in eachindex(a[k]) + mul!(suma, suma, dx, ordQ) + add!(suma, suma, a[k], ordQ) + end end - nothing + return suma end -function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{Taylor1{T},1}, - x0::AbstractArray{Taylor1{T}}) where {T<:NumberNotSeriesN} - @inbounds for i in eachindex(x, x0) - x0[i] = evaluate( x[i], δx ) +function evaluate(a::Taylor1{T}, dx::TaylorN{T}) where {T<:NumberNotSeries} + @inbounds suma = TaylorN( zero(T), get_order(dx) ) + @inbounds for k in reverse(eachindex(a)) + suma = suma*dx + a[k] end - nothing + return suma end -function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{TaylorN{T},1}, - x0::AbstractArray{TaylorN{T}}; sorting::Bool=true) where {T<:NumberNotSeriesN} - @inbounds for i in eachindex(x, x0) - x0[i] = _evaluate( x[i], δx, Val(sorting) ) +function evaluate(a::Taylor1{TaylorN{T}}, dx::TaylorN{T}) where {T<:NumberNotSeries} + order = min(minimum(get_order.(a[:])), get_order(dx)) + @inbounds suma = TaylorN( zero(T), order ) + @inbounds for k in reverse(eachindex(a)) + suma = suma*dx + a[k] end - nothing + return suma end -function evaluate!(x::AbstractArray{TaylorN{T}}, δt::T, - x0::AbstractArray{TaylorN{T}}) where {T<:Number} - @inbounds for i in eachindex(x, x0) - x0[i] = evaluate( x[i], δt ) - end - nothing -end +#function-like behavior for Taylor1 +(p::Taylor1)(x) = evaluate(p, x) +(p::Taylor1)() = evaluate(p) + +#function-like behavior for Vector{Taylor1} (asumes Julia version >= 1.6) +(p::Array{Taylor1{T}})(x) where {T<:Number} = evaluate.(p, x) +(p::SubArray{Taylor1{T}})(x) where {T<:Number} = evaluate.(p, x) +(p::Array{Taylor1{T}})() where {T<:Number} = evaluate.(p) +(p::SubArray{Taylor1{T}})() where {T<:Number} = evaluate.(p) """ @@ -154,43 +129,43 @@ Evaluate a `HomogeneousPolynomial` polynomial at `vals`. If `vals` is omitted, it's evaluated at zero. Note that the syntax `a(vals)` is equivalent to `evaluate(a, vals)`; and `a()` is equivalent to `evaluate(a)`. """ -function evaluate(a::HomogeneousPolynomial, vals::NTuple) +function evaluate(a::HomogeneousPolynomial, vals::NTuple{N,<:Number}) where {N} @assert length(vals) == get_numvars() - return _evaluate(a, vals) end -function _evaluate(a::HomogeneousPolynomial, vals::NTuple) +evaluate(a::HomogeneousPolynomial{T}, vals::AbstractArray{S,1} ) where + {T<:Number,S<:NumberNotSeriesN} = evaluate(a, (vals...,)) + +evaluate(a::HomogeneousPolynomial, v, vals::Vararg{Number,N}) where {N} = + evaluate(a, promote(v, vals...,)) + +evaluate(a::HomogeneousPolynomial, v) = evaluate(a, promote(v...,)) + +function evaluate(a::HomogeneousPolynomial{T}) where {T} + a.order == 0 && return a[1] + return zero(a[1]) +end + +# Internal method that avoids checking that the length of `vals` is the appropriate +function _evaluate(a::HomogeneousPolynomial{T}, vals::NTuple) where {T} + # @assert length(vals) == get_numvars() + a.order == 0 && return a[1]*one(vals[1]) ct = coeff_table[a.order+1] suma = zero(a[1])*vals[1] - for (i, a_coeff) in enumerate(a.coeffs) iszero(a_coeff) && continue @inbounds tmp = prod( vals .^ ct[i] ) suma += a_coeff * tmp end - return suma end -evaluate(a::HomogeneousPolynomial{T}, vals::AbstractArray{S,1} ) where - {T<:Number,S<:NumberNotSeriesN} = _evaluate(a, (vals...,)) - -evaluate(a::HomogeneousPolynomial, v, vals::Vararg) = evaluate(a, (v, vals...,)) - -evaluate(a::HomogeneousPolynomial, v) = evaluate(a, [v...]) - -function evaluate(a::HomogeneousPolynomial) - a.order == 0 && return a[1] - zero(a[1]) -end - #function-like behavior for HomogeneousPolynomial (p::HomogeneousPolynomial)(x) = evaluate(p, x) - -(p::HomogeneousPolynomial)(x, v::Vararg{T, N}) where {T,N} = evaluate(p, (x, v...,)) - +(p::HomogeneousPolynomial)(x, v::Vararg{Number,N}) where {N} = + evaluate(p, promote(x, v...,)) (p::HomogeneousPolynomial)() = evaluate(p) @@ -205,73 +180,82 @@ terms that are added. Note that the syntax `a(vals)` is equivalent to `evaluate(a, vals)`; and `a()` is equivalent to -`evaluate(a)`. No extension exists that incorporates -`sorting`. +`evaluate(a)`; use a(b::Bool, x) corresponds to +evaluate(a, x, sorting=b). """ -evaluate(a::TaylorN, vals; sorting::Bool=true) = _evaluate(a, (vals...,), Val(sorting)) - -evaluate(a::TaylorN, vals::NTuple; sorting::Bool=true) = _evaluate(a, vals, Val(sorting)) - -evaluate(a::TaylorN, v, vals::Vararg; sorting::Bool=true) = _evaluate(a, (v, vals...,), Val(sorting)) - -evaluate(a::TaylorN, vals::NTuple{N,<:AbstractSeries}; sorting::Bool=false) where - {N} = _evaluate(a, vals, Val(sorting)) - -evaluate(a::TaylorN{Taylor1}, vals::NTuple{N,<:AbstractSeries}; sorting::Bool=false) where - {N} = _evaluate(a, vals, Val(sorting)) +function evaluate(a::TaylorN, vals::NTuple{N,<:Number}; + sorting::Bool=true) where {N} + @assert get_numvars() == N + return _evaluate(a, vals, Val(sorting)) +end -evaluate(a::TaylorN{Taylor1}, vals::NTuple; sorting::Bool=false) = - _evaluate(a, vals, Val(sorting)) +function evaluate(a::TaylorN, vals::NTuple{N,<:AbstractSeries}; + sorting::Bool=false) where {N} + @assert get_numvars() == N + return _evaluate(a, vals, Val(sorting)) +end -evaluate(a::TaylorN{Taylor1}, v::Number, vals::Vararg; sorting::Bool=false) = - _evaluate(a, (v, vals...,), Val(sorting)) +evaluate(a::TaylorN{T}, vals::AbstractVector{<:Number}; + sorting::Bool=true) where {T<:NumberNotSeries} = + evaluate(a, (vals...,); sorting=sorting) -evaluate(a::TaylorN, vals::AbstractVector{T}) where {S, T<:AbstractSeries{S}} = - _evaluate(a, (vals...,), Val(false)) +evaluate(a::TaylorN{T}, vals::AbstractVector{<:AbstractSeries}; + sorting::Bool=false) where {T<:NumberNotSeries} = + evaluate(a, (vals...,); sorting=sorting) -evaluate(a::TaylorN{Taylor1}, vals::AbstractVector{T}) where {T} = - _evaluate(a, (vals...,), Val(false)) +evaluate(a::TaylorN{Taylor1{T}}, vals::AbstractVector{S}; + sorting::Bool=false) where {T, S} = + evaluate(a, (vals...,); sorting=sorting) -function evaluate(a::TaylorN{T}, s::Symbol, val::S) where {T<:Number, S<:NumberNotSeriesN} +function evaluate(a::TaylorN{T}, s::Symbol, val::S) where + {T<:Number, S<:NumberNotSeriesN} vars = get_variables(T) ind = lookupvar(s) - vars[ind] = val + zero(vars[ind]) - return evaluate(a, vars) + @inbounds vars[ind] = val + zero(vars[ind]) + return _evaluate(a, (vars...,), Val(false)) end -evaluate(a::TaylorN{T}, x::Pair{Symbol,S}) where {T<:NumberNotSeries,S} = +evaluate(a::TaylorN{T}, x::Pair{Symbol,S}) where {T<:NumberNotSeries, S} = evaluate(a, first(x), last(x)) -evaluate(a::TaylorN{T}, x::Pair{Symbol,S}) where {T<:AbstractSeries,S} = +evaluate(a::TaylorN{Taylor1{T}}, x::Pair{Symbol,S}) where {T<:NumberNotSeries, S} = evaluate(a, first(x), last(x)) -evaluate(a::TaylorN{T}) where {T<:Number} = a[0][1] +evaluate(a::TaylorN{T}) where {T<:Number} = constant_term(a) # _evaluate -function _evaluate(a::TaylorN, vals) - @assert get_numvars() == length(vals) - R = promote_type(numtype(a), typeof.(vals)...) +# Returns a vector with the evaluation of the HomogeneousPolynomials +function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:Number}) where {N,T<:Number} + R = promote_type(T, 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) + @inbounds for homPol in eachindex(a) + suma[homPol+1] = _evaluate(a[homPol], vals) end return suma end -function _evaluate(a::TaylorN{T}, vals::NTuple, ::Val{true}) where {T<:NumberNotSeries} - suma = _evaluate(a, vals) - return sum( sort!(suma, by=abs2) ) -end -function _evaluate(a::TaylorN{T}, vals::NTuple, ::Val{false}) where {T<:Number} - suma = _evaluate(a, vals) - return sum( suma ) +function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:TaylorN}) where {N,T<:Number} + R = TaylorN{promote_type(T, TS.numtype(vals[1]))} + a_length = length(a) + suma = zeros(R, a_length) + @inbounds for homPol in eachindex(a) + suma[homPol+1] = _evaluate(a[homPol], vals) + end + return suma end +_evaluate(a::TaylorN{T}, vals::NTuple, ::Val{true}) where {T<:NumberNotSeries} = + sum( sort!(_evaluate(a, vals), by=abs2) ) + +_evaluate(a::TaylorN{T}, vals::NTuple, ::Val{false}) where {T<:Number} = + sum( _evaluate(a, vals) ) + #High-dim array evaluation -function evaluate(A::AbstractArray{TaylorN{T},N}, δx::Vector{S}) where {T<:Number,S<:Number,N} +function evaluate(A::AbstractArray{TaylorN{T},N}, δx::Vector{S}) where + {T<:Number, S<:Number, N} R = promote_type(T,S) return evaluate(convert(Array{TaylorN{R},N},A), convert(Vector{R},δx)) end @@ -287,17 +271,45 @@ evaluate(A::AbstractArray{TaylorN{T}}) where {T<:Number} = evaluate.(A) (p::TaylorN)() = evaluate(p) (p::TaylorN)(s::Symbol, x) = evaluate(p, s, x) (p::TaylorN)(x::Pair) = evaluate(p, first(x), last(x)) -(p::TaylorN)(x, v::Vararg{T, N}) where {T,N} = evaluate(p, (x, v...,)) +(p::TaylorN)(x, v::Vararg{T}) where {T} = evaluate(p, (x, v...,)) (p::TaylorN)(b::Bool, x) = evaluate(p, x, sorting=b) -(p::TaylorN)(b::Bool, x, v::Vararg{T, N}) where {T,N} = evaluate(p, (x, v...,), sorting=b) +(p::TaylorN)(b::Bool, x, v::Vararg{T}) where {T} = evaluate(p, (x, v...,), sorting=b) #function-like behavior for AbstractArray{TaylorN{T}} -if VERSION > v"1.1" - (p::AbstractArray{TaylorN{T}})(x) where {T<:Number} = evaluate(p, x) - (p::AbstractArray{TaylorN{T}})() where {T<:Number} = evaluate(p) -else - (p::Array{TaylorN{T}})(x) where {T<:Number} = evaluate(p, x) - (p::SubArray{TaylorN{T}})(x) where {T<:Number} = evaluate(p, x) - (p::Array{TaylorN{T}})() where {T<:Number} = evaluate(p) - (p::SubArray{TaylorN{T}})() where {T<:Number} = evaluate(p) +(p::Array{TaylorN{T}})(x) where {T<:Number} = evaluate(p, x) +(p::SubArray{TaylorN{T}})(x) where {T<:Number} = evaluate(p, x) +(p::Array{TaylorN{T}})() where {T<:Number} = evaluate(p) +(p::SubArray{TaylorN{T}})() where {T<:Number} = evaluate(p) + + +""" + evaluate!(x, δt, x0) + +Evaluates each element of `x::AbstractArray{Taylor1{T}}`, +representing the Taylor expansion for the dependent variables +of an ODE at *time* `δt`. It updates the vector `x0` with the +computed values. +""" +function evaluate!(x::AbstractArray{Taylor1{T}}, δt::S, + x0::AbstractArray{T}) where {T<:Number, S<:Number} + x0 .= evaluate.( x, δt ) + return nothing +end + + +## In place evaluation of multivariable arrays +function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{T,1}, + x0::AbstractArray{T}) where {T<:Number} + @inbounds for i in eachindex(x, x0) + x0[i] = evaluate( x[i], δx ) + end + return nothing +end + +function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{TaylorN{T},1}, + x0::AbstractArray{TaylorN{T}}; sorting::Bool=true) where {T<:NumberNotSeriesN} + @inbounds for i in eachindex(x, x0) + x0[i] = _evaluate( x[i], δx, Val(sorting) ) + end + return nothing end diff --git a/src/functions.jl b/src/functions.jl index 4696b5d0..4af728bc 100644 --- a/src/functions.jl +++ b/src/functions.jl @@ -34,7 +34,7 @@ for T in (:Taylor1, :TaylorN) function log(a::$T) iszero(constant_term(a)) && throw(DomainError(a, - """The 0-th order coefficient must be non-zero in order to expand `log` around 0.""")) + """The 0-th order coefficient must be non-zero in order to expand `log` around 0.""")) order = a.order aux = log(constant_term(a)) @@ -118,7 +118,7 @@ for T in (:Taylor1, :TaylorN) function acos(a::$T) a0 = constant_term(a) a0^2 == one(a0) && throw(DomainError(a, - """Series expansion of asin(x) diverges at x = ±1.""")) + """Series expansion of asin(x) diverges at x = ±1.""")) order = a.order aux = acos(a0) @@ -247,10 +247,9 @@ for T in (:Taylor1, :TaylorN) end @inline function one!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number} + zero!(c, a, k) if k == 0 @inbounds c[0] = one(a[0]) - else - @inbounds c[k] = zero(a[k]) end return nothing end @@ -263,8 +262,8 @@ for T in (:Taylor1, :TaylorN) return subst!(c, a, k) else throw(DomainError(a, - """The 0th order coefficient must be non-zero - (abs(x) is not differentiable at x=0).""")) + """The 0th order coefficient must be non-zero + (abs(x) is not differentiable at x=0).""")) end return nothing end @@ -276,6 +275,7 @@ for T in (:Taylor1, :TaylorN) @inbounds c[0] = exp(constant_term(a)) return nothing end + zero!(c, a, k) if $T == Taylor1 @inbounds c[k] = k * a[k] * c[0] else @@ -297,6 +297,7 @@ for T in (:Taylor1, :TaylorN) @inbounds c[0] = expm1(constant_term(a)) return nothing end + zero!(c, a, k) c0 = c[0]+one(c[0]) if $T == Taylor1 @inbounds c[k] = k * a[k] * c0 @@ -322,7 +323,7 @@ for T in (:Taylor1, :TaylorN) @inbounds c[1] = a[1] / constant_term(a) return nothing end - + zero!(c, a, k) if $T == Taylor1 @inbounds c[k] = (k-1) * a[1] * c[k-1] else @@ -340,16 +341,19 @@ for T in (:Taylor1, :TaylorN) end @inline function log1p!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number} - a0 = constant_term(a) - a0p1 = a0+one(a0) if k == 0 + a0 = constant_term(a) @inbounds c[0] = log1p(a0) return nothing elseif k == 1 + a0 = constant_term(a) + a0p1 = a0+one(a0) @inbounds c[1] = a[1] / a0p1 return nothing end - + a0 = constant_term(a) + a0p1 = a0+one(a0) + zero!(c, a, k) if $T == Taylor1 @inbounds c[k] = (k-1) * a[1] * c[k-1] else @@ -372,8 +376,9 @@ for T in (:Taylor1, :TaylorN) @inbounds s[0], c[0] = sincos( a0 ) return nothing end - x = a[1] + zero!(s, a, k) + zero!(c, a, k) if $T == Taylor1 @inbounds s[k] = x * c[k-1] @inbounds c[k] = -x * s[k-1] @@ -403,7 +408,11 @@ for T in (:Taylor1, :TaylorN) @inbounds s[0], c[0] = sincospi( a0 ) return nothing end - sincos!(s, c, a*pi, k) + aa = $T(pi * a[0], a.order) + @inbounds for ordQ in eachindex(a) + aa[ordQ] = pi * a[ordQ] + end + sincos!(s, c, aa, k) return nothing end @@ -414,7 +423,7 @@ for T in (:Taylor1, :TaylorN) @inbounds c2[0] = aux^2 return nothing end - + zero!(c, a, k) if $T == Taylor1 @inbounds c[k] = k * a[k] * c2[0] else @@ -440,7 +449,7 @@ for T in (:Taylor1, :TaylorN) @inbounds r[0] = sqrt( 1 - a0^2 ) return nothing end - + zero!(c, a, k) if $T == Taylor1 @inbounds c[k] = (k-1) * r[1] * c[k-1] else @@ -465,7 +474,7 @@ for T in (:Taylor1, :TaylorN) @inbounds r[0] = sqrt( 1 - a0^2 ) return nothing end - + zero!(c, a, k) if $T == Taylor1 @inbounds c[k] = (k-1) * r[1] * c[k-1] else @@ -490,7 +499,7 @@ for T in (:Taylor1, :TaylorN) @inbounds r[0] = 1 + a0^2 return nothing end - + zero!(c, a, k) if $T == Taylor1 @inbounds c[k] = (k-1) * r[1] * c[k-1] else @@ -514,8 +523,9 @@ for T in (:Taylor1, :TaylorN) @inbounds c[0] = cosh( constant_term(a) ) return nothing end - x = a[1] + zero!(s, a, k) + zero!(c, a, k) if $T == Taylor1 @inbounds s[k] = x * c[k-1] @inbounds c[k] = x * s[k-1] @@ -545,7 +555,7 @@ for T in (:Taylor1, :TaylorN) @inbounds c2[0] = aux^2 return nothing end - + zero!(c, a, k) if $T == Taylor1 @inbounds c[k] = k * a[k] * c2[0] else @@ -571,7 +581,7 @@ for T in (:Taylor1, :TaylorN) @inbounds r[0] = sqrt( a0^2 + 1 ) return nothing end - + zero!(c, a, k) if $T == Taylor1 @inbounds c[k] = (k-1) * r[1] * c[k-1] else @@ -596,7 +606,7 @@ for T in (:Taylor1, :TaylorN) @inbounds r[0] = sqrt( a0^2 - 1 ) return nothing end - + zero!(c, a, k) if $T == Taylor1 @inbounds c[k] = (k-1) * r[1] * c[k-1] else @@ -621,7 +631,7 @@ for T in (:Taylor1, :TaylorN) @inbounds r[0] = 1 - a0^2 return nothing end - + zero!(c, a, k) if $T == Taylor1 @inbounds c[k] = (k-1) * r[1] * c[k-1] else @@ -644,6 +654,493 @@ for T in (:Taylor1, :TaylorN) end +# Mutating functions for Taylor1{TaylorN{T}} +@inline function exp!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, + k::Int) where {T<:NumberNotSeries} + if k == 0 + @inbounds for ordQ in eachindex(a[0]) + # zero!(res[0], a[0], ordQ) + exp!(res[0], a[0], ordQ) + end + return nothing + end + # The recursion formula + tmp = TaylorN( zero(a[k][0][1]), a[0].order) + zero!(res, a, k) + for i = 0:k-1 + @inbounds for ordQ in eachindex(a[0]) + tmp[ordQ] = (k-i) * res[i][ordQ] + mul!(res[k], tmp, a[k-i], ordQ) + end + end + div!(res, res, k, k) + return nothing +end + +@inline function expm1!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, + k::Int) where {T<:NumberNotSeries} + if k == 0 + @inbounds for ordQ in eachindex(a[0]) + # zero!(res[0], a[0], ordQ) + expm1!(res[0], a[0], ordQ) + end + return nothing + end + # The recursion formula + tmp = TaylorN( zero(a[k][0][1]), a[0].order) + zero!(res, a, k) + # i=0 term of sum + @inbounds for ordQ in eachindex(a[0]) + one!(tmp, a[0], ordQ) + add!(tmp, res[0], tmp, ordQ) + tmp[ordQ] = k * tmp[ordQ] + # zero!(res[k], a[0], ordQ) + mul!(res[k], a[k], tmp, ordQ) + end + for i = 1:k-1 + @inbounds for ordQ in eachindex(a[0]) + tmp[ordQ] = (k-i) * res[i][ordQ] + mul!(res[k], tmp, a[k-i], ordQ) + end + end + div!(res, res, k, k) + return nothing +end + +@inline function log!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, + k::Int) where {T<:NumberNotSeries} + if k == 0 + @inbounds for ordQ in eachindex(a[0]) + # zero!(res[k], a[0], ordQ) + log!(res[k], a[0], ordQ) + end + return nothing + elseif k == 1 + @inbounds for ordQ in eachindex(a[0]) + zero!(res[k], a[0], ordQ) + div!(res[k], a[1], a[0], ordQ) + end + return nothing + end + # The recursion formula + tmp = TaylorN( zero(a[k][0][1]), a[0].order) + zero!(res, a, k) + for i = 1:k-1 + @inbounds for ordQ in eachindex(a[0]) + tmp[ordQ] = (k-i) * res[k-i][ordQ] + mul!(res[k], tmp, a[i], ordQ) + end + end + div!(res, res, k, k) + @inbounds for ordQ in eachindex(a[0]) + subst!(tmp, a[k], res[k], ordQ) + zero!(res[k], a[0], ordQ) + div!(res[k], tmp, a[0], ordQ) + end + return nothing +end + +@inline function log1p!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, + k::Int) where {T<:NumberNotSeries} + if k == 0 + @inbounds for ordQ in eachindex(a[0]) + # zero!(res[k], a[0], ordQ) + log1p!(res[k], a[0], ordQ) + end + return nothing + end + tmp1 = TaylorN( zero(a[k][0][1]), a[0].order) + zero!(res, a, k) + @inbounds for ordQ in eachindex(a[0]) + # zero!(res[k], a[0], ordQ) + one!(tmp1, a[0], ordQ) + add!(tmp1, tmp1, a[0], ordQ) + end + if k == 1 + @inbounds for ordQ in eachindex(a[0]) + div!(res[k], a[1], tmp1, ordQ) + end + return nothing + end + # The recursion formula + tmp = TaylorN( zero(a[k][0][1]), a[0].order) + for i = 1:k-1 + @inbounds for ordQ in eachindex(a[0]) + tmp[ordQ] = (k-i) * res[k-i][ordQ] + mul!(res[k], tmp, a[i], ordQ) + end + end + div!(res, res, k, k) + @inbounds for ordQ in eachindex(a[0]) + subst!(tmp, a[k], res[k], ordQ) + zero!(res[k], a[0], ordQ) + div!(res[k], tmp, tmp1, ordQ) + end + return nothing +end + +@inline function sincos!(s::Taylor1{TaylorN{T}}, c::Taylor1{TaylorN{T}}, + a::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries} + if k == 0 + @inbounds for ordQ in eachindex(a[0]) + sincos!(s[0], c[0], a[0], ordQ) + end + return nothing + end + # The recursion formula + x = TaylorN( a[1][0][1], a[0].order ) + zero!(s, a, k) + zero!(c, a, k) + @inbounds for i = 1:k + for ordQ in eachindex(a[0]) + x[ordQ] = i * a[i][ordQ] + mul!(s[k], x, c[k-i], ordQ) + mul!(c[k], x, s[k-i], ordQ) + end + end + div!(s, s, k, k) + subst!(c, c, k) + div!(c, c, k, k) + return nothing +end + +@inline function sincospi!(s::Taylor1{TaylorN{T}}, c::Taylor1{TaylorN{T}}, + a::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries} + if k == 0 + @inbounds for ordQ in eachindex(a[0]) + sincospi!(s[0], c[0], a[0], ordQ) + end + return nothing + end + # aa = pi * a + aa = Taylor1(zero(a[0]), a.order) + @inbounds for ordT in eachindex(a) + mul!(aa, pi, a, ordT) + end + sincos!(s, c, aa, k) + return nothing +end + +@inline function tan!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, + res2::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries} + if k == 0 + @inbounds res[0] = tan( a[0] ) + # zero!(res2, res, 0) + sqr!(res2, res, 0) + return nothing + end + # The recursion formula + tmp = TaylorN( zero(a[0][0][1]), a[0].order) + zero!(res, a, k) + for i = 0:k-1 + @inbounds for ordQ in eachindex(a[0]) + tmp[ordQ] = (k-i) * res2[i][ordQ] + mul!(res[k], tmp, a[k-i], ordQ) + end + end + @inbounds for ordQ in eachindex(a[0]) + # zero!(tmp, res[k], ordQ) + tmp[ordQ] = res[k][ordQ] / k + add!(res[k], a[k], tmp, ordQ) + end + # zero!(res2, res, k) + sqr!(res2, res, k) + return nothing +end + +@inline function asin!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, + r::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries} + if k == 0 + @inbounds res[0] = asin( a[0] ) + # r[0] = sqrt(1-a[0]^2) + tmp = TaylorN( zero(a[0][0][1]), a[0].order) + r[0] = square(a[0]) + for ordQ in eachindex(a[0]) + one!(tmp, a[0], ordQ) + subst!(tmp, tmp, r[0], ordQ) + # zero!(r[0], tmp, ordQ) + sqrt!(r[0], tmp, ordQ) + end + return nothing + end + # The recursion formula + tmp = TaylorN( zero(a[0][0][1]), a[0].order) + zero!(res, a, k) + for i in 1:k-1 + @inbounds for ordQ in eachindex(a[0]) + tmp[ordQ] = (k-i) * res[k-i][ordQ] + mul!(res[k], tmp, r[i], ordQ) + end + end + div!(res, res, k, k) + @inbounds for ordQ in eachindex(a[0]) + subst!(tmp, a[k], res[k], ordQ) + zero!(res[k], a[0], ordQ) + div!(res[k], tmp, r[0], ordQ) + end + # Compute auxiliary term s=1-a^2 + s = Taylor1(zero(a[0]), a.order) + for i = 0:k + sqr!(s, a, i) + subst!(s, s, i) + if i == 0 + s[0] = one(s[0]) - s[0] + add!(s, one(s), s, 0) + end + end + # Update aux term r = sqrt(s) = sqrt(1-a^2) + sqrt!(r, s, k) + return nothing +end + +@inline function acos!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, + r::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries} + if k == 0 + @inbounds res[0] = acos( a[0] ) + # r[0] = sqrt(1-a[0]^2) + tmp = TaylorN( zero(a[0][0][1]), a[0].order) + r[0] = square(a[0]) + for ordQ in eachindex(a[0]) + one!(tmp, a[0], ordQ) + subst!(tmp, tmp, r[0], ordQ) + # zero!(r[0], tmp, ordQ) + sqrt!(r[0], tmp, ordQ) + end + return nothing + end + # The recursion formula + tmp = TaylorN( zero(a[0][0][1]), a[0].order) + zero!(res, a, k) + for i in 1:k-1 + @inbounds for ordQ in eachindex(a[0]) + tmp[ordQ] = (k-i) * res[k-i][ordQ] + mul!(res[k], tmp, r[i], ordQ) + end + end + div!(res, res, k, k) + @inbounds for ordQ in eachindex(a[0]) + add!(tmp, a[k], res[k], ordQ) + subst!(tmp, tmp, ordQ) + zero!(res[k], a[0], ordQ) + div!(res[k], tmp, r[0], ordQ) + end + # Compute auxiliary term s=1-a^2 + s = Taylor1(zero(a[0]), a.order) + for i = 0:k + sqr!(s, a, i) + subst!(s, s, i) + if i == 0 + s[0] = one(s[0]) - s[0] + add!(s, one(s), s, 0) + end + end + # Update aux term r = sqrt(s) = sqrt(1-a^2) + sqrt!(r, s, k) + return nothing +end + +@inline function atan!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, + r::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries} + if k == 0 + res[0] = atan( a[0] ) + # zero!(r, a, 0) + sqr!(r, a, 0) + add!(r, r, one(a[0][0][1]), 0) + return nothing + end + # The recursion formula + tmp = TaylorN( zero(a[0][0][1]), a[0].order ) + zero!(res, a, k) + for i in 1:k-1 + @inbounds for ordQ in eachindex(a[0]) + tmp[ordQ] = (k-i) * res[k-i][ordQ] + mul!(res[k], tmp, r[i], ordQ) + end + end + @inbounds for ordQ in eachindex(a[0]) + # zero!(tmp, res[k], ordQ) + tmp[ordQ] = - res[k][ordQ] / k + add!(tmp, a[k], tmp, ordQ) + zero!(res[k], a[0], ordQ) + div!(res[k], tmp, r[0], ordQ) + end + zero!(r, a, k) + sqr!(r, a, k) + return nothing +end + +@inline function sinhcosh!(s::Taylor1{TaylorN{T}}, c::Taylor1{TaylorN{T}}, + a::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries} + if k == 0 + @inbounds for ordQ in eachindex(a[0]) + sinhcosh!(s[0], c[0], a[0], ordQ) + end + return nothing + end + # The recursion formula + x = TaylorN( a[k][0][1], a[0].order ) + zero!(s, a, k) + zero!(c, a, k) + @inbounds for i = 1:k + for ordQ in eachindex(a[0]) + x[ordQ] = i * a[i][ordQ] + mul!(s[k], x, c[k-i], ordQ) + mul!(c[k], x, s[k-i], ordQ) + end + end + div!(s, s, k, k) + div!(c, c, k, k) + return nothing +end + +@inline function tanh!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, + res2::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries} + if k == 0 + @inbounds res[0] = tanh( a[0] ) + # zero!(res2, res, 0) + sqr!(res2, res, 0) + return nothing + end + # The recursion formula + tmp = TaylorN( zero(a[0][0][1]), a[0].order) + zero!(res, a, k) + for i = 0:k-1 + @inbounds for ordQ in eachindex(a[0]) + tmp[ordQ] = (k-i) * res2[i][ordQ] + mul!(res[k], tmp, a[k-i], ordQ) + end + end + @inbounds for ordQ in eachindex(a[0]) + # zero!(tmp, res[k], ordQ) + tmp[ordQ] = res[k][ordQ] / k + subst!(res[k], a[k], tmp, ordQ) + end + zero!(res2, res, k) + sqr!(res2, res, k) + return nothing +end + +@inline function asinh!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, + r::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries} + if k == 0 + @inbounds res[0] = asinh( a[0] ) + # r[0] = sqrt(1+a[0]^2) + tmp = TaylorN( zero(a[0][0][1]), a[0].order) + r[0] = square(a[0]) + for ordQ in eachindex(a[0]) + one!(tmp, a[0], ordQ) + add!(tmp, tmp, r[0], ordQ) + # zero!(r[0], tmp, ordQ) + sqrt!(r[0], tmp, ordQ) + end + return nothing + end + # The recursion formula + tmp = TaylorN( zero(a[0][0][1]), a[0].order) + zero!(res, a, k) + for i in 1:k-1 + @inbounds for ordQ in eachindex(a[0]) + tmp[ordQ] = (k-i) * res[k-i][ordQ] + mul!(res[k], tmp, r[i], ordQ) + end + end + div!(res, res, k, k) + @inbounds for ordQ in eachindex(a[0]) + subst!(tmp, a[k], res[k], ordQ) + zero!(res[k], a[0], ordQ) + div!(res[k], tmp, r[0], ordQ) + end + # Compute auxiliary term s=1+a^2 + s = Taylor1(zero(a[0]), a.order) + for i = 0:k + sqr!(s, a, i) + if i == 0 + s[0] = one(s[0]) + s[0] + add!(s, one(s), s, 0) + end + end + # Update aux term r = sqrt(s) = sqrt(1+a^2) + sqrt!(r, s, k) + return nothing +end + +@inline function acosh!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, + r::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries} + if k == 0 + @inbounds res[0] = acosh( a[0] ) + # r[0] = sqrt(a[0]^2-1) + tmp = TaylorN( zero(a[0][0][1]), a[0].order) + r[0] = square(a[0]) + for ordQ in eachindex(a[0]) + one!(tmp, a[0], ordQ) + subst!(tmp, r[0], tmp, ordQ) + # zero!(r[0], tmp, ordQ) + sqrt!(r[0], tmp, ordQ) + end + return nothing + end + # The recursion formula + tmp = TaylorN( zero(a[0][0][1]), a[0].order) + zero!(res, a, k) + for i in 1:k-1 + @inbounds for ordQ in eachindex(a[0]) + tmp[ordQ] = (k-i) * res[k-i][ordQ] + mul!(res[k], tmp, r[i], ordQ) + end + end + div!(res, res, k, k) + @inbounds for ordQ in eachindex(a[0]) + subst!(tmp, a[k], res[k], ordQ) + zero!(res[k], a[0], ordQ) + div!(res[k], tmp, r[0], ordQ) + end + # Compute auxiliary term s=a^2-1 + s = Taylor1(zero(a[0]), a.order) + for i = 0:k + sqr!(s, a, i) + if i == 0 + s[0] = one(s[0]) + s[0] + subst!(s, s, one(s), 0) + end + end + # Update aux term r = sqrt(s) = sqrt(a^2-1) + sqrt!(r, s, k) + return nothing +end + +@inline function atanh!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, + r::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries} + if k == 0 + res[0] = atanh( a[0] ) + # zero!(r, a, 0) + sqr!(r, a, 0) + subst!(r, one(a[0][0][1]), r, 0) + return nothing + end + # The recursion formula + tmp = TaylorN( zero(a[0][0][1]), a[0].order ) + zero!(res, a, k) + for i in 1:k-1 + @inbounds for ordQ in eachindex(a[0]) + tmp[ordQ] = (k-i) * res[k-i][ordQ] + mul!(res[k], tmp, r[i], ordQ) + end + end + @inbounds for ordQ in eachindex(a[0]) + # zero!(tmp, res[k], ordQ) + tmp[ordQ] = res[k][ordQ] / k + add!(tmp, a[k], tmp, ordQ) + zero!(res[k], a[0], ordQ) + div!(res[k], tmp, r[0], ordQ) + end + zero!(r, a, k) + sqr!(r, a, k) + return nothing +end + + + + @doc doc""" inverse(f) diff --git a/src/hash_tables.jl b/src/hash_tables.jl index 08dbb9dd..4420e306 100644 --- a/src/hash_tables.jl +++ b/src/hash_tables.jl @@ -121,8 +121,8 @@ List the indices and corresponding of a `HomogeneousPolynomial` of degree `ord`. """ function show_monomials(ord::Int) - z = zeros(Int, TaylorSeries.size_table[ord+1]) - for (index, value) in enumerate(TaylorSeries.coeff_table[ord+1]) + z = zeros(Int, TS.size_table[ord+1]) + for (index, value) in enumerate(TS.coeff_table[ord+1]) z[index] = 1 pol = HomogeneousPolynomial(z) println(" $index --> $(homogPol2str(pol)[4:end])") diff --git a/src/other_functions.jl b/src/other_functions.jl index 531c457e..8d225851 100644 --- a/src/other_functions.jl +++ b/src/other_functions.jl @@ -215,40 +215,36 @@ a `TaylorN` expansion will be computed. If the dimension of x0 (`length(x0)`) is different from the variables set for `TaylorN` (`get_numvars()`), an `AssertionError` will be thrown. """ -function taylor_expand(f::Function; order::Int=15) +function taylor_expand(f::F; order::Int=15) where {F} a = Taylor1(order) return f(a) end -function taylor_expand(f::Function, x0::T; order::Int=15) where {T<:Number} - a = Taylor1([x0, one(T)], order) +function taylor_expand(f::F, x0::T; order::Int=15) where {T<:Number, F} + a = Taylor1(T[x0, one(x0)], order) return f(a) end #taylor_expand function for TaylorN -function taylor_expand(f::Function, x0::Vector{T}; order::Int=get_order()) where {T<:Number} +function taylor_expand(f::F, x0::Vector{T}; order::Int=get_order()) where {T<:Number, F} ll = length(x0) @assert ll == get_numvars() && order <= get_order() X = Array{TaylorN{T}}(undef, ll) - for i in eachindex(X) X[i] = x0[i] + TaylorN(T, i, order=order) end - return f( X ) end -function taylor_expand(f::Function, x0::Vararg; order::Int=get_order()) - x0 = promote(x0...) +function taylor_expand(f::F, x1::Vararg{Number,N}; order::Int=get_order()) where {F, N} + x0 = promote(x1...) ll = length(x0) T = eltype(x0[1]) @assert ll == get_numvars() && order <= get_order() X = Array{TaylorN{T}}(undef, ll) - for i in eachindex(X) X[i] = x0[i] + TaylorN(T, i, order=order) end - return f( X... ) end diff --git a/src/power.jl b/src/power.jl index 2d3ba7c7..8026eb44 100644 --- a/src/power.jl +++ b/src/power.jl @@ -54,6 +54,11 @@ for T in (:Taylor1, :TaylorN) @eval ^(a::$T, x::T) where {T<:Complex} = exp( x*log(a) ) end +^(a::Taylor1{TaylorN{T}}, n::Integer) where {T<:NumberNotSeries} = a^float(n) + +^(a::Taylor1{TaylorN{T}}, r::Rational) where {T<:NumberNotSeries} = a^(r.num/r.den) + + # power_by_squaring; slightly modified from base/intfuncs.jl # Licensed under MIT "Expat" @@ -83,8 +88,9 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) end end + ## Real power ## -function ^(a::Taylor1, r::S) where {S<:Real} +function ^(a::Taylor1{T}, r::S) where {T<:Number, S<:Real} a0 = constant_term(a) aux = one(a0)^r @@ -92,16 +98,15 @@ function ^(a::Taylor1, r::S) where {S<:Real} aa = aux*a r == 1 && return aa r == 2 && return square(aa) - r == 1/2 && return sqrt(aa) + r == 0.5 && 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 + (lnull > a.order) && return Taylor1( zero(aux), a.order) + c_order = l0 == 0 ? a.order : min(a.order, trunc(Int,r*a.order)) c = Taylor1(zero(aux), c_order) - for k = 0:c_order + for k in eachindex(c) pow!(c, aa, r, k) end @@ -117,22 +122,47 @@ function ^(a::TaylorN, r::S) where {S<:Real} aa = aux*a r == 1 && return aa r == 2 && return square(aa) - r == one(r)/2 && return sqrt(aa) - isinteger(r) && return aa^round(Int,r) + r == 0.5 && return sqrt(aa) + isinteger(r) && return aa^round(Int,r) # uses power_by_squaring - # @assert !iszero(a0) iszero(a0) && throw(DomainError(a, """The 0-th order TaylorN coefficient must be non-zero in order to expand `^` around 0.""")) - c = TaylorN( a0^r, a.order) - for ord in 1:a.order + c = TaylorN( zero(aux), a.order) + for ord in eachindex(a) pow!(c, aa, r, ord) end return c end +function ^(a::Taylor1{TaylorN{T}}, r::S) where {T<:NumberNotSeries, S<:Real} + a0 = constant_term(a) + aux = one(a0)^r + + iszero(r) && return Taylor1(aux, a.order) + aa = aux*a + r == 1 && return aa + r == 2 && return square(aa) + r == 0.5 && return sqrt(aa) + # Is the following needed? + # isinteger(r) && r > 0 && iszero(constant_term(a[0])) && + # return power_by_squaring(aa, round(Int,r)) + + l0 = findfirst(a) + lnull = trunc(Int, r*l0 ) + (lnull > a.order) && return Taylor1( zero(aux), a.order) + + c_order = l0 == 0 ? a.order : min(a.order, trunc(Int,r*a.order)) + c = Taylor1(zero(aux), c_order) + for k in eachindex(c) + pow!(c, aa, r, k) + end + + return c +end + # Homogeneous coefficients for real power @doc doc""" @@ -152,7 +182,8 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`. """ pow! -@inline function pow!(c::Taylor1{T}, a::Taylor1{T}, r::S, k::Int) where {T<:Number,S<:Real} +@inline function pow!(c::Taylor1{T}, a::Taylor1{T}, r::S, k::Int) where + {T<:Number, S<:Real} if r == 0 return one!(c, a, k) @@ -164,12 +195,12 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`. return sqrt!(c, a, k) end + # Sanity + zero!(c, a, k) + # First non-zero coefficient l0 = findfirst(a) - if l0 < 0 - c[k] = zero(a[0]) - return nothing - end + l0 < 0 && return nothing # The first non-zero coefficient of the result; must be integer !isinteger(r*l0) && throw(DomainError(a, @@ -177,16 +208,10 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`. to raise the Taylor1 polynomial to a non-integer exponent.""")) lnull = trunc(Int, r*l0 ) kprime = k-lnull - if (kprime < 0) || (lnull > a.order) - @inbounds c[k] = zero(a[0]) - return nothing - end + (kprime < 0 || lnull > a.order) && return nothing # Relevant for positive integer r, to avoid round-off errors - if isinteger(r) && r > 0 && (k > r*findlast(a)) - @inbounds c[k] = zero(a[0]) - return nothing - end + isinteger(r) && r > 0 && (k > r*findlast(a)) && return nothing if k == lnull @inbounds c[k] = (a[l0])^r @@ -196,8 +221,6 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`. # The recursion formula if l0+kprime ≤ a.order @inbounds c[k] = r * kprime * c[lnull] * a[l0+kprime] - else - @inbounds c[k] = zero(a[0]) end for i = 1:k-lnull-1 ((i+lnull) > a.order || (l0+kprime-i > a.order)) && continue @@ -209,7 +232,7 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`. end @inline function pow!(c::TaylorN{T}, a::TaylorN{T}, r::S, k::Int) where - {T<:NumberNotSeriesN,S<:Real} + {T<:NumberNotSeriesN, S<:Real} if r == 0 return one!(c, a, k) @@ -226,6 +249,10 @@ end return nothing end + # Sanity + zero!(c, a, k) + + # The recursion formula for i = 0:k-1 aux = r*(k-i) - i mul!(c[k], aux*a[k-i], c[i]) @@ -235,6 +262,60 @@ end return nothing end +@inline function pow!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, r::S, + ordT::Int) where {T<:NumberNotSeries, S<:Real} + + if r == 0 + return one!(res, a, ordT) + elseif r == 1 + return identity!(res, a, ordT) + elseif r == 2 + return sqr!(res, a, ordT) + elseif r == 0.5 + return sqrt!(res, a, ordT) + end + + # Sanity + zero!(res, a, ordT) + + # First non-zero coefficient + l0 = findfirst(a) + l0 < 0 && return nothing + + # The first non-zero coefficient of the result; must be integer + !isinteger(r*l0) && throw(DomainError(a, + """The 0-th order Taylor1 coefficient must be non-zero + to raise the Taylor1 polynomial to a non-integer exponent.""")) + lnull = trunc(Int, r*l0 ) + kprime = ordT-lnull + (kprime < 0 || lnull > a.order) && return nothing + + # Relevant for positive integer r, to avoid round-off errors + isinteger(r) && r > 0 && (ordT > r*findlast(a)) && return nothing + + if ordT == lnull + res[ordT] = a[l0]^r # if r is integer, uses power_by_squaring internally + return nothing + end + + # The recursion formula + tmp = TaylorN( zero(constant_term(a[ordT])), a[0].order) + tmp1 = TaylorN( zero(constant_term(a[ordT])), a[0].order) + for i = 0:ordT-lnull-1 + ((i+lnull) > a.order || (l0+kprime-i > a.order)) && continue + aux = r*(kprime-i) - i + @inbounds for ordQ in eachindex(tmp) + tmp1[ordQ] = aux * res[i+lnull][ordQ] + mul!(tmp, tmp1, a[l0+kprime-i], ordQ) + end + end + @inbounds for ordQ in eachindex(tmp) + tmp1[ordQ] = tmp[ordQ]/kprime + div!(res[ordT], tmp1, a[l0], ordQ) + end + + return nothing +end ## Square ## @@ -255,11 +336,19 @@ for T in (:Taylor1, :TaylorN) end function square(a::HomogeneousPolynomial) - order = 2*a.order - order > get_order() && return HomogeneousPolynomial(zero(a[1]), get_order()) - + order = 2*get_order(a) + # NOTE: the following returns order 0, but could be get_order(), or get_order(a) + order > get_order() && return HomogeneousPolynomial(zero(a[1]), 0) res = HomogeneousPolynomial(zero(a[1]), order) - sqr!(res, a) + accsqr!(res, a) + return res +end + +function square(a::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries} + res = Taylor1(zero(a[0]), a.order) + for ordT in eachindex(a) + sqr!(res, a, ordT) + end return res end @@ -292,8 +381,12 @@ for T = (:Taylor1, :TaylorN) return nothing end + # Sanity + zero!(c, a, k) + + # Recursion formula kodd = k%2 - kend = div(k - 2 + kodd, 2) + kend = (k - 2 + kodd) >> 1 @inbounds for i = 0:kend if $T == Taylor1 c[k] += a[i] * a[k-i] @@ -305,9 +398,9 @@ for T = (:Taylor1, :TaylorN) kodd == 1 && return nothing if $T == Taylor1 - @inbounds c[k] += a[div(k,2)]^2 + @inbounds c[k] += a[k >> 1]^2 else - sqr!(c[k], a[div(k,2)]) + accsqr!(c[k], a[k >> 1]) end return nothing @@ -315,14 +408,47 @@ for T = (:Taylor1, :TaylorN) end end +@inline function sqr!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, + ordT::Int) where {T<:NumberNotSeries} + # Sanity + zero!(res, a, ordT) + if ordT == 0 + @inbounds for ordQ in eachindex(a[0]) + @inbounds sqr!(res[0], a[0], ordQ) + end + return nothing + end + + # Recursion formula + kodd = ordT%2 + kend = (ordT - 2 + kodd) >> 1 + tmp = TaylorN( zero(constant_term(a[0])), a[0].order ) + for i = 0:kend + @inbounds for ordQ in eachindex(a[ordT]) + mul!(tmp, a[i], a[ordT-i], ordQ) + end + end + @inbounds for ordQ in eachindex(a[ordT]) + res[ordT][ordQ] = 2 * tmp[ordQ] + end + kodd == 1 && return nothing + + @inbounds for ordQ in eachindex(a[0]) + sqr!(tmp, a[ordT >> 1], ordQ) + add!(res[ordT], res[ordT], tmp, ordQ) + end + + return nothing +end + """ - sqr!(c, a) + accsqr!(c, a) -Return `c = a*a` with no allocation; all parameters are `HomogeneousPolynomial`. +Returns `c += a*a` with no allocation; all parameters are `HomogeneousPolynomial`. """ -@inline function sqr!(c::HomogeneousPolynomial{T}, a::HomogeneousPolynomial{T}) where +@inline function accsqr!(c::HomogeneousPolynomial{T}, a::HomogeneousPolynomial{T}) where {T<:NumberNotSeriesN} iszero(a) && return nothing @@ -352,26 +478,24 @@ end ## Square root ## -function sqrt(a::Taylor1) - +function sqrt(a::Taylor1{T}) where {T<:Number} # First non-zero coefficient l0nz = findfirst(a) aux = zero(sqrt( constant_term(a) )) if l0nz < 0 return Taylor1(aux, a.order) - elseif l0nz%2 == 1 # l0nz must be pair + elseif isodd(l0nz) # 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.""")) + """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. + # The last l0nz coefficients are dropped. 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 + aa = convert(Taylor1{eltype(aux)}, a) + for k in eachindex(c) sqrt!(c, aa, k, lnull) end @@ -382,19 +506,44 @@ function sqrt(a::TaylorN) @inbounds p0 = sqrt( constant_term(a) ) if iszero(p0) throw(DomainError(a, - """The 0-th order TaylorN coefficient must be non-zero - in order to expand `sqrt` around 0.""")) + """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 + aa = convert(TaylorN{eltype(p0)}, a) + for k in eachindex(c) sqrt!(c, aa, k) end return c end +function sqrt(a::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries} + # First non-zero coefficient + l0nz = findfirst(a) + aux = TaylorN( zero(sqrt(constant_term(a[0]))), a[0].order ) + if l0nz < 0 + return Taylor1( aux, a.order ) + elseif isodd(l0nz) # 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 dropped. + lnull = l0nz >> 1 # integer division by 2 + c_order = l0nz == 0 ? a.order : a.order >> 1 + c = Taylor1( aux, c_order ) + aa = convert(Taylor1{eltype(aux)}, a) + for k in eachindex(c) + sqrt!(c, aa, k, lnull) + end + + return c +end + + # Homogeneous coefficients for the square-root @doc doc""" sqrt!(c, a, k::Int, k0::Int=0) @@ -419,14 +568,20 @@ coefficient, which must be even. """ sqrt! @inline function sqrt!(c::Taylor1{T}, a::Taylor1{T}, k::Int, k0::Int=0) where {T<:Number} + # Sanity + zero!(c, a, k) + + k < k0 && return nothing if k == k0 @inbounds c[k] = sqrt(a[2*k0]) return nothing end + # Recursion formula kodd = (k - k0)%2 - kend = div(k - k0 - 2 + kodd, 2) + # kend = div(k - k0 - 2 + kodd, 2) + kend = (k - k0 - 2 + kodd) >> 1 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] ) @@ -447,14 +602,17 @@ coefficient, which must be even. end @inline function sqrt!(c::TaylorN{T}, a::TaylorN{T}, k::Int) where {T<:NumberNotSeriesN} + # Sanity + zero!(c, a, k) if k == 0 @inbounds c[0] = sqrt( constant_term(a) ) return nothing end + # Recursion formula kodd = k%2 - kend = div(k - 2 + kodd, 2) + kend = (k - 2 + kodd) >> 1 @inbounds for i = 1:kend mul!(c[k], c[i], c[k-i]) end @@ -466,3 +624,50 @@ end return nothing end + +@inline function sqrt!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, ordT::Int, + ordT0::Int=0) where {T<:NumberNotSeries} + # Sanity + zero!(res, a, ordT) + ordT < ordT0 && return nothing + + if ordT == ordT0 + for ordQ in eachindex(a[0]) + sqrt!(res[ordT], a[2*ordT0], ordQ) + end + return nothing + end + + # Recursion formula + @inbounds tmp = TaylorN( zero( constant_term(a[ordT])), a[0].order) + @inbounds tmp1 = TaylorN( zero( constant_term(a[ordT])), a[0].order) + ordT_odd = (ordT - ordT0)%2 + ordT_end = (ordT - ordT0 - 2 + ordT_odd) >> 1 + imax = min(ordT0+ordT_end, a.order) + imin = max(ordT0+1, ordT+ordT0-a.order) + if imin ≤ imax + @inbounds for ordQ in eachindex(a[0]) + mul!(tmp, res[imin], res[ordT+ordT0-imin], ordQ) + end + end + for i = imin+1:imax + @inbounds for ordQ in eachindex(a[0]) + mul!(tmp, res[i], res[ordT+ordT0-i], ordQ) + end + end + + @inbounds for ordQ in eachindex(a[0]) + tmp[ordQ] = -2 * tmp[ordQ] + if ordT+ordT0 ≤ a.order + add!(tmp, a[ordT+ordT0], tmp, ordQ) + end + if ordT_odd == 0 + sqr!(tmp1, res[ordT_end+ordT0+1], ordQ) + subst!(tmp, tmp, tmp1, ordQ) + end + tmp1[ordQ] = 2*res[ordT0][ordQ] + div!(res[ordT], tmp, tmp1, ordQ) + end + + return nothing +end diff --git a/test/manyvariables.jl b/test/manyvariables.jl index 42ed7b7f..6c91d637 100644 --- a/test/manyvariables.jl +++ b/test/manyvariables.jl @@ -426,7 +426,6 @@ end @test imag((exp(yT))^(-1im)') == sin(yT) exy = exp( xT+yT ) @test evaluate(exy) == 1 - @test evaluate(exy, 0.1im, 0.01im) == exp(0.11im) @test exy(0.1im, 0.01im) == exp(0.11im) @test evaluate(exy,(0.1im, 0.01im)) == exp(0.11im) @test exy((0.1im, 0.01im)) == exp(0.11im) diff --git a/test/mixtures.jl b/test/mixtures.jl index ed100739..02792a25 100644 --- a/test/mixtures.jl +++ b/test/mixtures.jl @@ -7,8 +7,8 @@ using Test # using LinearAlgebra @testset "Tests with mixtures of Taylor1 and TaylorN" begin - @test TaylorSeries.NumberNotSeries == Union{Real,Complex} - @test TaylorSeries.NumberNotSeriesN == Union{Real,Complex,Taylor1} + @test TS.NumberNotSeries == Union{Real,Complex} + @test TS.NumberNotSeriesN == Union{Real,Complex,Taylor1} set_variables("x", numvars=2, order=6) xH = HomogeneousPolynomial(Int, 1) @@ -30,6 +30,7 @@ using Test @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³)" @@ -82,12 +83,12 @@ using Test @test t1N[0] == HomogeneousPolynomial(1) ctN1 = convert(TaylorN{Taylor1{Float64}}, t1N) @test convert(eltype(tN1), tN1) === tN1 - @test eltype(xHt) == HomogeneousPolynomial{Taylor1{Float64}} @test eltype(tN1) == TaylorN{Taylor1{Float64}} @test eltype(Taylor1([xH])) == Taylor1{HomogeneousPolynomial{Int}} @test TS.numtype(xHt) == Taylor1{Float64} @test TS.numtype(tN1) == Taylor1{Float64} @test TS.numtype(Taylor1([xH])) == HomogeneousPolynomial{Int} + @test TS.numtype(t1N) == TaylorN{Float64} @test normalize_taylor(tN1) == tN1 @test get_order(HomogeneousPolynomial([Taylor1(1), 1.0+Taylor1(2)])) == 1 @test 3*tN1 == TaylorN([HomogeneousPolynomial([3t]),3xHt,3yHt^2]) @@ -155,10 +156,12 @@ using Test @test t1N() == TaylorN(xH, 2) @test string(evaluate(t1N, 0.0)) == " 1.0 x₁ + 𝒪(‖x‖³)" @test string(evaluate(t1N^2, 1.0)) == " 1.0 + 2.0 x₁ + 1.0 x₁² + 2.0 x₂² + 𝒪(‖x‖³)" - @test string((t1N^2)(1.0)) == " 1.0 + 2.0 x₁ + 1.0 x₁² + 2.0 x₂² + 𝒪(‖x‖³)" + @test string((t1N^(2//1))(1.0)) == " 1.0 + 2.0 x₁ + 1.0 x₁² + 2.0 x₂² + 𝒪(‖x‖³)" v = zeros(TaylorN{Float64},2) - @test evaluate!([t1N, t1N^2], 0.0, v) == nothing + @test isnothing(evaluate!([t1N, t1N^2], 0.0, v)) @test v == [TaylorN(1), TaylorN(1)^2] + @test tN1() == t + @test evaluate(tN1, :x₁ => 1.0) == TaylorN([HomogeneousPolynomial([1.0+t]), zero(xHt), yHt^2]) # Tests for functions of mixtures t1N = Taylor1([zero(TaylorN(Float64,1)), one(TaylorN(Float64,1))], 6) @@ -180,9 +183,18 @@ using Test @test x*fn(cc+t1N) == fn(cc+t)*xN1 @test t*fn(cc+xN1) == fn(cc+x)*t1N end + ee = Taylor1(t1N[0:5], 6) + for ord in eachindex(t1N) + TS.differentiate!(ee, exp(t1N), ord) + end + @test iszero(ee[6]) + @test getcoeff.(ee, 0:5) == getcoeff.(exp(t1N), 0:5) + ee = differentiate(t1N, get_order(t1N)) + @test iszero(ee) + @test iszero(get_order(ee)) vt = zeros(Taylor1{Float64},2) - @test evaluate!([tN1, tN1^2], [t, t], vt) == nothing + @test isnothing(evaluate!([tN1, tN1^2], [t, t], vt)) @test vt == [2t, 4t^2] tint = Taylor1(Int, 10) @@ -199,12 +211,25 @@ using Test @test (1+y)/one(t) == 1 + y @test typeof(y+t) == TaylorN{Taylor1{Float64}} + t = Taylor1(4) + xN, yN = get_variables() + @test evaluate(1.0 + t + t^2, xN) == 1.0 + xN + xN^2 + v1 = [1.0 + t + t^2 + t^4, 1.0 - t^2 + t^3] + @test v1(yN^2) == [1.0 + yN^2 + yN^4, 1.0 - yN^4 + yN^6] + tN = Taylor1([zero(xN), one(xN)], 4) + q1N = 1 + yN*tN + xN*tN^4 + @test q1N(-1.0) == 1.0 - yN + xN + @test q1N(-xN^2) == 1.0 - xN^2*yN + # See #92 and #94 δx, δy = set_variables("δx δy") - xx = 1+Taylor1(δx,5) + xx = 1+Taylor1(δx, 5) + yy = 1+Taylor1(δy, 5) + tt = Taylor1([zero(δx), one(δx)], xx.order) + @test all((xx, yy) .== (TS.fixorder(xx, yy))) @test typeof(xx) == Taylor1{TaylorN{Float64}} @test eltype(xx) == Taylor1{TaylorN{Float64}} - @test TS.numtype(xx) == TaylorN{Float64} + @test TS.numtype(tt) == TaylorN{Float64} @test !(@inferred isnan(xx)) @test !(@inferred isnan(δx)) @test !(@inferred isinf(xx)) @@ -218,7 +243,88 @@ using Test @test xx/xx == one(xx) @test xx*δx + Taylor1(typeof(δx),5) == δx + δx^2 + Taylor1(typeof(δx),5) @test xx/(1+δx) == one(xx) + @test 1/(1-tt) == 1 + tt + tt^2 + tt^3 + tt^4 + tt^5 + @test xx/(1-tt) == xx * (1 + tt + tt^2 + tt^3 + tt^4 + tt^5) + res = xx * tt + @test 1/(1-xx*tt) == 1 + res + res^2 + res^3 + res^4 + res^5 @test typeof(xx+δx) == Taylor1{TaylorN{Float64}} + res = 1/(1+δx) + @test one(xx)/(xx*(1+tt)) == Taylor1([res, -res, res, -res, res, -res]) + res = 1/(1+δx)^2 + @test (xx^2 + yy^2)/(xx*yy) == xx/yy + yy/xx + @test ((xx+yy)*tt)^2/((xx+yy)*tt) == (xx+yy)*tt + @test sqrt(xx) == Taylor1(sqrt(xx[0]), xx.order) + @test xx^0.25 == sqrt(sqrt(xx)) + @test (xx*yy*tt^2)^0.5 == sqrt(xx*yy)*tt + FF(x,y,t) = (1 + x + y + t)^4 + QQ(x,y,t) = x^4.0 + (4*y + (4*t + 4))*x^3 + (6*y^2 + (12*t + 12)*y + (6*t^2 + 12*t + 6))*x^2 + + (4*y^3 + (12*t + 12)*y^2 + (12*t^2 + 24*t + 12)*y + (4*t^3 + 12*t^2 + 12*t + 4))*x + + (y^4 + (4*t + 4)*y^3 + (6*t^2 + 12*t + 6)*y^2 + (4*t^3 + 12*t^2 + 12*t + 4)*y + + (t^4 + 4*t^3 + 6*t^2 + 4*t + 1)) + @test FF(xx, yy, tt) == QQ(xx, yy, tt) + @test FF(tt, yy-1, xx-1) == QQ(xx-1, yy-1, tt) + @test (xx+tt)^4 == xx^4 + 4*xx^3*tt + 6*xx^2*tt^2 + 4*xx*tt^3 + tt^4 + pp = xx*yy*(1+tt)^4 + @test pp^0.25 == sqrt(sqrt(pp)) + @test (xx*yy)^(3/2)*(1+tt+tt^2) == (sqrt(xx*yy))^3*(1+tt+tt^2) + @test sqrt((xx+yy+tt)^3) ≈ (xx+yy+tt)^1.5 + @test (sqrt(xx*(1+tt)))^4 ≈ xx^2 * (1+tt)^2 + @test (sqrt(xx*(1+tt)))^5 ≈ xx^2.5 * (1+tt)^2.5 + @test exp(xx) == exp(xx[0]) + zero(tt) + @test exp(xx+tt) ≈ exp(xx)*exp(tt) + @test norm(exp(xx+tt) - exp(xx)*exp(tt), Inf) < 5e-17 + @test expm1(xx) == expm1(xx[0]) + zero(tt) + @test expm1(xx+tt) ≈ exp(xx+tt)-1 + @test norm(expm1(xx+tt) - (exp(xx+tt)-1), Inf) < 5e-16 + @test log(xx) == log(xx[0]) + zero(tt) + @test log((xx+tt)*(yy+tt)) ≈ log(xx+tt)+log(yy+tt) + @test log(pp) ≈ log(xx)+log(yy)+4*log(1+tt) + @test norm(log(pp) - (log(xx)+log(yy)+4*log(1+tt)), Inf) < 1e-15 + @test log1p(xx) == log1p(xx[0]) + zero(tt) + @test log1p(xx+tt) == log(1+xx+tt) + exp(log(xx+tt)) == log(exp(xx+tt)) + @test exp(log(pp)) ≈ log(exp(pp)) + @test sincos(δx + xx) == sincos(δx+xx[0]) .+ zero(tt) + qq = sincos(pp) + @test exp(im*pp) ≈ qq[2] + im*qq[1] + @test qq[2]^2 ≈ 1 - qq[1]^2 + @test sincospi(δx + xx - 1) == sincospi(δx + xx[0] - 1) .+ zero(tt) + @test all(sincospi(pp) .≈ sincos(pi*pp)) + @test tan(xx) == tan(xx[0]) + zero(tt) + @test tan(xx+tt) ≈ sin(xx+tt)/cos(xx+tt) + @test tan(xx*tt) ≈ sin(xx*tt)/cos(xx*tt) + @test asin(xx-1) == asin(xx[0]-1) + zero(tt) + @test asin(sin(xx+tt)) ≈ xx + tt + @test sin(asin(xx*tt)) == xx * tt + @test_throws DomainError asin(2+xx+tt) + @test acos(xx-1) == acos(xx[0]-1) + zero(tt) + @test acos(cos(xx+tt)) ≈ xx + tt + @test cos(acos(xx*tt)) ≈ xx * tt + @test_throws DomainError acos(2+xx+tt) + @test atan(xx) == atan(xx[0]) + zero(tt) + @test atan(tan(xx+tt)) ≈ xx + tt + @test tan(atan(xx*tt)) == xx * tt + @test TS.sinhcosh(xx) == TS.sinhcosh(xx[0]) .+ zero(tt) + qq = TS.sinhcosh(pp) + @test qq[2]^2 - 1 ≈ qq[1]^2 + @test 2*qq[1] ≈ exp(pp)-exp(-pp) + @test 2*qq[2] ≈ exp(pp)+exp(-pp) + qq = TS.sinhcosh(xx+tt) + @test tanh(xx) == tanh(xx[0]) + zero(tt) + @test tanh(xx+tt) ≈ qq[1]/qq[2] + qq = TS.sinhcosh(xx*tt) + @test tanh(xx*tt) ≈ qq[1]/qq[2] + @test asinh(xx-1) == asinh(xx[0]-1) + zero(tt) + @test asinh(sinh(xx+tt)) ≈ xx + tt + @test sinh(asinh(xx*tt)) == xx * tt + @test acosh(xx+1) == acosh(xx[0]+1) + zero(tt) + @test acosh(cosh(xx+1+tt)) ≈ xx + 1 + tt + @test cosh(acosh(2+xx*tt)) ≈ 2 + xx * tt + @test_throws DomainError acosh(xx+tt) + @test atanh(xx-1) == atanh(xx[0]-1) + zero(tt) + @test atanh(tanh(-1+xx+tt)) ≈ -1 + xx + tt + @test tanh(atanh(xx*tt)) ≈ xx * tt + @test_throws DomainError atanh(xx+tt) #testing evaluate and function-like behavior of Taylor1, TaylorN for mixtures: t = Taylor1(25) @@ -295,15 +401,15 @@ using Test @test norm(-10X+4Y,Inf) == 10. - @test TaylorSeries.rtoldefault(TaylorN{Taylor1{Int}}) == 0 - @test TaylorSeries.rtoldefault(Taylor1{TaylorN{Int}}) == 0 + @test TS.rtoldefault(TaylorN{Taylor1{Int}}) == 0 + @test TS.rtoldefault(Taylor1{TaylorN{Int}}) == 0 for T in (Float64, BigFloat) - @test TaylorSeries.rtoldefault(TaylorN{Taylor1{T}}) == sqrt(eps(T)) - @test TaylorSeries.rtoldefault(Taylor1{TaylorN{T}}) == sqrt(eps(T)) - @test TaylorSeries.real(TaylorN{Taylor1{T}}) == TaylorN{Taylor1{T}} - @test TaylorSeries.real(Taylor1{TaylorN{T}}) == Taylor1{TaylorN{T}} - @test TaylorSeries.real(TaylorN{Taylor1{Complex{T}}}) == TaylorN{Taylor1{T}} - @test TaylorSeries.real(Taylor1{TaylorN{Complex{T}}}) == Taylor1{TaylorN{T}} + @test TS.rtoldefault(TaylorN{Taylor1{T}}) == sqrt(eps(T)) + @test TS.rtoldefault(Taylor1{TaylorN{T}}) == sqrt(eps(T)) + @test TS.real(TaylorN{Taylor1{T}}) == TaylorN{Taylor1{T}} + @test TS.real(Taylor1{TaylorN{T}}) == Taylor1{TaylorN{T}} + @test TS.real(TaylorN{Taylor1{Complex{T}}}) == TaylorN{Taylor1{T}} + @test TS.real(Taylor1{TaylorN{Complex{T}}}) == Taylor1{TaylorN{T}} end rndT1(ord1) = Taylor1(-1 .+ 2rand(ord1+1)) # generates a random Taylor1 with order `ord` @@ -341,7 +447,7 @@ using Test Pv = [rndTN(get_order(), 3), rndTN(get_order(), 3)] Qv = convert.(Taylor1{TaylorN{Float64}}, Pv) - @test TaylorSeries.jacobian(Pv) == TaylorSeries.jacobian(Qv) + @test TS.jacobian(Pv) == TS.jacobian(Qv) @test_throws ArgumentError Taylor1(2) + TaylorN(1) @test_throws ArgumentError Taylor1(2) - TaylorN(1) diff --git a/test/mutatingfuncts.jl b/test/mutatingfuncts.jl index ca569027..98de7745 100644 --- a/test/mutatingfuncts.jl +++ b/test/mutatingfuncts.jl @@ -8,84 +8,84 @@ using Test t1 = Taylor1(6) # Dictionaries with calls - @test length(TaylorSeries._dict_binary_ops) == 5 - # @test length(TaylorSeries._dict_unary_ops) == 22 # why are these tested? + @test length(TS._dict_binary_ops) == 5 + # @test length(TS._dict_unary_ops) == 22 # why are these tested? - @test all([haskey(TaylorSeries._dict_binary_ops, op) + @test all([haskey(TS._dict_binary_ops, op) for op in [:+, :-, :*, :/, :^]]) - @test all([haskey(TaylorSeries._dict_binary_calls, op) + @test all([haskey(TS._dict_binary_calls, op) for op in [:+, :-, :*, :/, :^]]) - for kk in keys(TaylorSeries._dict_binary_ops) - res = TaylorSeries._internalmutfunc_call( - TaylorSeries._InternalMutFuncs(TaylorSeries._dict_binary_ops[kk]) ) - @test TaylorSeries._dict_binary_calls[kk] == res + for kk in keys(TS._dict_binary_ops) + res = TS._internalmutfunc_call( + TS._InternalMutFuncs(TS._dict_binary_ops[kk]) ) + @test TS._dict_binary_calls[kk] == res end - for kk in keys(TaylorSeries._dict_unary_ops) - res = TaylorSeries._internalmutfunc_call( - TaylorSeries._InternalMutFuncs(TaylorSeries._dict_unary_ops[kk]) ) - @test TaylorSeries._dict_unary_calls[kk] == res + for kk in keys(TS._dict_unary_ops) + res = TS._internalmutfunc_call( + TS._InternalMutFuncs(TS._dict_unary_ops[kk]) ) + @test TS._dict_unary_calls[kk] == res end # Some examples t1 = Taylor1(5) t2 = zero(t1) - TaylorSeries.pow!(t2, t1, 2, 2) + TS.pow!(t2, t1, 2, 2) @test t2[2] == 1.0 # res = zero(t1) - TaylorSeries.add!(res, t1, t2, 3) + TS.add!(res, t1, t2, 3) @test res[3] == 0.0 - TaylorSeries.add!(res, 1, t2, 3) + TS.add!(res, 1, t2, 3) @test res[3] == 0.0 - TaylorSeries.add!(res, t2, 3, 0) + TS.add!(res, t2, 3, 0) @test res[0] == 3.0 - TaylorSeries.subst!(res, t1, t2, 2) + TS.subst!(res, t1, t2, 2) @test res[2] == -1.0 - TaylorSeries.subst!(res, t1, 1, 0) + TS.subst!(res, t1, 1, 0) @test res[0] == -1.0 @test res[2] == -1.0 - TaylorSeries.subst!(res, 1, t2, 2) + TS.subst!(res, 1, t2, 2) @test res[2] == -1.0 res[3] = rand() - TaylorSeries.mul!(res, t1, t2, 3) + TS.mul!(res, t1, t2, 3) @test res[3] == 1.0 - TaylorSeries.mul!(res, res, 2, 3) + TS.mul!(res, res, 2, 3) @test res[3] == 2.0 - TaylorSeries.mul!(res, 0.5, res, 3) + TS.mul!(res, 0.5, res, 3) @test res[3] == 1.0 res[0] = rand() - TaylorSeries.div!(res, t2-1, 1+t1, 0) + TS.div!(res, t2-1, 1+t1, 0) res[1] = rand() - TaylorSeries.div!(res, t2-1, 1+t1, 1) + TS.div!(res, t2-1, 1+t1, 1) @test res[0] == (t1-1)[0] @test res[1] == (t1-1)[1] - TaylorSeries.div!(res, res, 2, 0) + TS.div!(res, res, 2, 0) @test res[0] == -0.5 res = zero(t1) - TaylorSeries.identity!(res, t1, 0) + TS.identity!(res, t1, 0) @test res[0] == t1[0] - TaylorSeries.zero!(res, t1, 0) - TaylorSeries.zero!(res, t1, 1) + TS.zero!(res, t1, 0) + TS.zero!(res, t1, 1) @test res[0] == zero(t1[0]) @test res[1] == zero(t1[1]) - TaylorSeries.one!(res, t1, 0) - TaylorSeries.one!(res, t1, 0) + TS.one!(res, t1, 0) + TS.one!(res, t1, 0) @test res[0] == one(t1[0]) @test res[1] == zero(t1[1]) res = zero(t1) - TaylorSeries.abs!(res, -1-t2, 2) + TS.abs!(res, -1-t2, 2) @test res[2] == 1.0 - @test_throws DomainError TaylorSeries.abs!(res, t2, 2) + @test_throws DomainError TS.abs!(res, t2, 2) res = zero(t1) - TaylorSeries.abs2!(res, 1-t1, 1) + TS.abs2!(res, 1-t1, 1) @test res[1] == -2.0 - TaylorSeries.abs2!(res, t1, 2) + TS.abs2!(res, t1, 2) @test res[2] == 1.0 end diff --git a/test/onevariable.jl b/test/onevariable.jl index 9aaa826a..abd231d8 100644 --- a/test/onevariable.jl +++ b/test/onevariable.jl @@ -51,11 +51,11 @@ Base.iszero(::SymbNumber) = false end v = [1,2] - @test typeof(TaylorSeries.resize_coeffs1!(v,3)) == Nothing + @test typeof(TS.resize_coeffs1!(v,3)) == Nothing @test v == [1,2,0,0] - TaylorSeries.resize_coeffs1!(v,0) + TS.resize_coeffs1!(v,0) @test v == [1] - TaylorSeries.resize_coeffs1!(v,3) + TS.resize_coeffs1!(v,3) setindex!(Taylor1(v),3,2) @test v == [1,0,3,0] pol_int = Taylor1(v) @@ -139,10 +139,10 @@ Base.iszero(::SymbNumber) = false @test TS.numtype(promote(1.0+im, zt)[1]) == Complex{Float64} @test length(Taylor1(10)) == 11 - @test length.( TaylorSeries.fixorder(zt, Taylor1([1])) ) == (16, 16) - @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 length.( TS.fixorder(zt, Taylor1([1])) ) == (16, 16) + @test length.( TS.fixorder(zt, Taylor1([1], 1)) ) == (2, 2) + @test eltype(TS.fixorder(zt,Taylor1([1]))[1]) == Taylor1{Int} + @test TS.numtype(TS.fixorder(zt,Taylor1([1]))[1]) == Int @test findfirst(t) == 1 @test findfirst(t^2) == 2 @test findfirst(ot) == 0 @@ -380,86 +380,86 @@ Base.iszero(::SymbNumber) = false ut = 1.0*t tt = zero(ut) - TaylorSeries.one!(tt, ut, 0) + TS.one!(tt, ut, 0) @test tt[0] == 1.0 - TaylorSeries.one!(tt, ut, 1) + TS.one!(tt, ut, 1) @test tt[1] == 0.0 - TaylorSeries.abs!(tt, 1.0+ut, 0) + TS.abs!(tt, 1.0+ut, 0) @test tt[0] == 1.0 - TaylorSeries.add!(tt, ut, ut, 1) + TS.add!(tt, ut, ut, 1) @test tt[1] == 2.0 - TaylorSeries.add!(tt, -3.0, 0) + TS.add!(tt, -3.0, 0) @test tt[0] == -3.0 - TaylorSeries.add!(tt, -3.0, 1) + TS.add!(tt, -3.0, 1) @test tt[1] == 0.0 - TaylorSeries.subst!(tt, ut, ut, 1) + TS.subst!(tt, ut, ut, 1) @test tt[1] == 0.0 - TaylorSeries.subst!(tt, -3.0, 0) + TS.subst!(tt, -3.0, 0) @test tt[0] == 3.0 - TaylorSeries.subst!(tt, -2.5, 1) + TS.subst!(tt, -2.5, 1) @test tt[1] == 0.0 - iind, cind = TaylorSeries.divfactorization(ut, ut) + iind, cind = TS.divfactorization(ut, ut) @test iind == 1 @test cind == 1.0 - TaylorSeries.div!(tt, ut, ut, 0) + TS.div!(tt, ut, ut, 0) @test tt[0] == cind - TaylorSeries.div!(tt, 1+ut, 1+ut, 0) + TS.div!(tt, 1+ut, 1+ut, 0) @test tt[0] == 1.0 - TaylorSeries.div!(tt, 1, 1+ut, 0) + TS.div!(tt, 1, 1+ut, 0) @test tt[0] == 1.0 - TaylorSeries.pow!(tt, 1.0+t, 1.5, 0) + TS.pow!(tt, 1.0+t, 1.5, 0) @test tt[0] == 1.0 - TaylorSeries.pow!(tt, 0.0*t, 1.5, 0) + TS.pow!(tt, 0.0*t, 1.5, 0) @test tt[0] == 0.0 - TaylorSeries.pow!(tt, 0.0+t, 18, 0) + TS.pow!(tt, 0.0+t, 18, 0) @test tt[0] == 0.0 - TaylorSeries.pow!(tt, 1.0+t, 1.5, 0) + TS.pow!(tt, 1.0+t, 1.5, 0) @test tt[0] == 1.0 - TaylorSeries.pow!(tt, 1.0+t, 0.5, 1) + TS.pow!(tt, 1.0+t, 0.5, 1) @test tt[1] == 0.5 - TaylorSeries.pow!(tt, 1.0+t, 0, 0) + TS.pow!(tt, 1.0+t, 0, 0) @test tt[0] == 1.0 - TaylorSeries.pow!(tt, 1.0+t, 1, 1) + TS.pow!(tt, 1.0+t, 1, 1) @test tt[1] == 1.0 tt = zero(ut) - TaylorSeries.pow!(tt, 1.0+t, 2, 0) + TS.pow!(tt, 1.0+t, 2, 0) @test tt[0] == 1.0 - TaylorSeries.pow!(tt, 1.0+t, 2, 1) + TS.pow!(tt, 1.0+t, 2, 1) @test tt[1] == 2.0 - TaylorSeries.pow!(tt, 1.0+t, 2, 2) + TS.pow!(tt, 1.0+t, 2, 2) @test tt[2] == 1.0 - TaylorSeries.sqrt!(tt, 1.0+t, 0, 0) + TS.sqrt!(tt, 1.0+t, 0, 0) @test tt[0] == 1.0 - TaylorSeries.sqrt!(tt, 1.0+t, 0) + TS.sqrt!(tt, 1.0+t, 0) @test tt[0] == 1.0 - TaylorSeries.exp!(tt, 1.0*t, 0) + TS.exp!(tt, 1.0*t, 0) @test tt[0] == exp(t[0]) - TaylorSeries.log!(tt, 1.0+t, 0) + TS.log!(tt, 1.0+t, 0) @test tt[0] == 0.0 - TaylorSeries.log1p!(tt, 0.25+t, 0) + TS.log1p!(tt, 0.25+t, 0) @test tt[0] == log1p(0.25) - TaylorSeries.log1p!(tt, 0.25+t, 1) + TS.log1p!(tt, 0.25+t, 1) @test tt[1] == 1/1.25 ct = zero(ut) - TaylorSeries.sincos!(tt, ct, 1.0*t, 0) + TS.sincos!(tt, ct, 1.0*t, 0) @test tt[0] == sin(t[0]) @test ct[0] == cos(t[0]) - TaylorSeries.tan!(tt, 1.0*t, ct, 0) + TS.tan!(tt, 1.0*t, ct, 0) @test tt[0] == tan(t[0]) @test ct[0] == tan(t[0])^2 - TaylorSeries.asin!(tt, 1.0*t, ct, 0) + TS.asin!(tt, 1.0*t, ct, 0) @test tt[0] == asin(t[0]) @test ct[0] == sqrt(1.0-t[0]^2) - TaylorSeries.acos!(tt, 1.0*t, ct, 0) + TS.acos!(tt, 1.0*t, ct, 0) @test tt[0] == acos(t[0]) @test ct[0] == sqrt(1.0-t[0]^2) - TaylorSeries.atan!(tt, ut, ct, 0) + TS.atan!(tt, ut, ct, 0) @test tt[0] == atan(t[0]) @test ct[0] == 1.0+t[0]^2 - TaylorSeries.sinhcosh!(tt, ct, ut, 0) + TS.sinhcosh!(tt, ct, ut, 0) @test tt[0] == sinh(t[0]) @test ct[0] == cosh(t[0]) - TaylorSeries.tanh!(tt, ut, ct, 0) + TS.tanh!(tt, ut, ct, 0) @test tt[0] == tanh(t[0]) @test ct[0] == tanh(t[0])^2 @@ -543,12 +543,12 @@ Base.iszero(::SymbNumber) = false displayBigO(false) @test string(ta(-3)) == " - 3 + 1 t " @test string(ta(0)^3-3) == " - 3 + 1 t³ " - @test TaylorSeries.pretty_print(ta(3im)) == " ( 0 + 3im ) + ( 1 + 0im ) t " + @test TS.pretty_print(ta(3im)) == " ( 0 + 3im ) + ( 1 + 0im ) t " @test string(Taylor1([1,2,3,4,5], 2)) == string(Taylor1([1,2,3])) displayBigO(true) @test string(ta(-3)) == " - 3 + 1 t + 𝒪(t¹⁶)" @test string(ta(0)^3-3) == " - 3 + 1 t³ + 𝒪(t¹⁶)" - @test TaylorSeries.pretty_print(ta(3im)) == " ( 0 + 3im ) + ( 1 + 0im ) t + 𝒪(t¹⁶)" + @test TS.pretty_print(ta(3im)) == " ( 0 + 3im ) + ( 1 + 0im ) t + 𝒪(t¹⁶)" @test string(Taylor1([1,2,3,4,5], 2)) == string(Taylor1([1,2,3])) a = collect(1:12) @@ -562,11 +562,11 @@ Base.iszero(::SymbNumber) = false @test norm(t_a, Inf) == 12 @test norm(t_C) == norm(complex(3.0,4.0)*a) - @test TaylorSeries.rtoldefault(Taylor1{Int}) == 0 - @test TaylorSeries.rtoldefault(Taylor1{Float64}) == sqrt(eps(Float64)) - @test TaylorSeries.rtoldefault(Taylor1{BigFloat}) == sqrt(eps(BigFloat)) - @test TaylorSeries.real(Taylor1{Float64}) == Taylor1{Float64} - @test TaylorSeries.real(Taylor1{Complex{Float64}}) == Taylor1{Float64} + @test TS.rtoldefault(Taylor1{Int}) == 0 + @test TS.rtoldefault(Taylor1{Float64}) == sqrt(eps(Float64)) + @test TS.rtoldefault(Taylor1{BigFloat}) == sqrt(eps(BigFloat)) + @test TS.real(Taylor1{Float64}) == Taylor1{Float64} + @test TS.real(Taylor1{Complex{Float64}}) == Taylor1{Float64} @test isfinite(t_C) @test isfinite(t_a) @test !isfinite( Taylor1([0, Inf]) ) @@ -596,10 +596,10 @@ Base.iszero(::SymbNumber) = false a = Taylor1(rand(10)) b = Taylor1(rand(10)) c = deepcopy(a) - TaylorSeries.deg2rad!(b, a, 0) + TS.deg2rad!(b, a, 0) @test a == c @test a[0]*(pi/180) == b[0] - # TaylorSeries.deg2rad!.(b, a, [0,1,2]) + # TS.deg2rad!.(b, a, [0,1,2]) # @test a == c # for i in 0:2 # @test a[i]*(pi/180) == b[i] @@ -607,10 +607,10 @@ Base.iszero(::SymbNumber) = false a = Taylor1(rand(10)) b = Taylor1(rand(10)) c = deepcopy(a) - TaylorSeries.rad2deg!(b, a, 0) + TS.rad2deg!(b, a, 0) @test a == c @test a[0]*(180/pi) == b[0] - # TaylorSeries.rad2deg!.(b, a, [0,1,2]) + # TS.rad2deg!.(b, a, [0,1,2]) # @test a == c # for i in 0:2 # @test a[i]*(180/pi) == b[i] diff --git a/test/runtests.jl b/test/runtests.jl index 217ccb15..707c02bd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,8 +15,6 @@ testfiles = ( ) for file in testfiles - # Skipping tests with intervals, since IntervalArithmetic.jl requires Julia v1.1+ - VERSION < v"1.1" && file == "intervals.jl" && continue include(file) end # After `using intervalArithmetic` new ambiguities may arise