diff --git a/Project.toml b/Project.toml index e25d96af..e135ba06 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.17.8" +version = "0.18.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -12,14 +12,14 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [weakdeps] IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" -StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [extensions] TaylorSeriesIAExt = "IntervalArithmetic" TaylorSeriesJLD2Ext = "JLD2" -TaylorSeriesSAExt = "StaticArrays" TaylorSeriesRATExt = "RecursiveArrayTools" +TaylorSeriesSAExt = "StaticArrays" [compat] Aqua = "0.8" diff --git a/README.md b/README.md index e5813543..2670e837 100644 --- a/README.md +++ b/README.md @@ -30,7 +30,7 @@ julia> t = Taylor1(Float64, 5) julia> exp(t) 1.0 + 1.0 t + 0.5 t² + 0.16666666666666666 t³ + 0.041666666666666664 t⁴ + 0.008333333333333333 t⁵ + 𝒪(t⁶) - + julia> log(1 + t) 1.0 t - 0.5 t² + 0.3333333333333333 t³ - 0.25 t⁴ + 0.2 t⁵ + 𝒪(t⁶) ``` @@ -40,7 +40,7 @@ julia> x, y = set_variables("x y", order=2); julia> exp(x + y) 1.0 + 1.0 x + 1.0 y + 0.5 x² + 1.0 x y + 0.5 y² + 𝒪(‖x‖³) - + ``` Differential and integral calculus on Taylor series: ```julia @@ -95,6 +95,6 @@ We thank the participants of the course for putting up with the half-baked material and contributing energy and ideas. We acknowledge financial support from DGAPA-UNAM PAPIME grants -PE-105911 and PE-107114, and DGAPA-PAPIIT grants IG-101113, -IG-100616, and IG-100819. +PE-105911 and PE-107114, and DGAPA-PAPIIT grants IG-101113, +IG-100616, IG-100819 and IG-101122. LB acknowledges support through a *Cátedra Moshinsky* (2013). diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 4c8c0a9a..e44a856a 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -162,9 +162,7 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) end for T in (:Taylor1, :TaylorN) - @eval function zero(a::$T) - return $T(zero.(a.coeffs)) - end + @eval zero(a::$T) = $T(zero.(a.coeffs)) @eval function one(a::$T) b = zero(a) b[0] = one(b[0]) @@ -539,25 +537,50 @@ for T in (:Taylor1, :TaylorN) return nothing end - @eval @inline function mul!(v::$T, a::$T, b::NumberNotSeries, k::Int) + @eval begin if $T == Taylor1 - @inbounds v[k] = a[k] * b - else - @inbounds for i in eachindex(v[k]) - v[k][i] = a[k][i] * b + @inline function mul!(v::$T, a::$T, b::NumberNotSeries, k::Int) + @inbounds v[k] = a[k] * b + return nothing + end + @inline function mul!(v::$T, a::NumberNotSeries, b::$T, k::Int) + @inbounds v[k] = a * b[k] + return nothing + end + @inline function muladd!(v::$T, a::$T, b::NumberNotSeries, k::Int) + @inbounds v[k] += a[k] * b + return nothing + end + @inline function muladd!(v::$T, a::NumberNotSeries, b::$T, k::Int) + @inbounds v[k] += a * b[k] + return nothing end - end - return nothing - end - @eval @inline function mul!(v::$T, a::NumberNotSeries, b::$T, k::Int) - if $T == Taylor1 - @inbounds v[k] = a * b[k] else - @inbounds for i in eachindex(v[k]) - v[k][i] = a * b[k][i] + @inline function mul!(v::$T, a::$T, b::NumberNotSeries, k::Int) + @inbounds for i in eachindex(v[k]) + v[k][i] = a[k][i] * b + end + return nothing + end + @inline function mul!(v::$T, a::NumberNotSeries, b::$T, k::Int) + @inbounds for i in eachindex(v[k]) + v[k][i] = a * b[k][i] + end + return nothing + end + @inline function muladd!(v::$T, a::$T, b::NumberNotSeries, k::Int) + @inbounds for i in eachindex(v[k]) + v[k][i] += a[k][i] * b + end + return nothing + end + @inline function muladd!(v::$T, a::NumberNotSeries, b::$T, k::Int) + @inbounds for i in eachindex(v[k]) + v[k][i] += a * b[k][i] + end + return nothing end end - return nothing end @eval @inline function mul!(v::$T, a::$T, b::NumberNotSeries) @@ -951,7 +974,8 @@ end @inline function div!(c::Taylor1{T}, a::NumberNotSeries, b::Taylor1{T}, k::Int) where {T<:Number} - iszero(a) && !iszero(b) && zero!(c, k) + zero!(c, k) + iszero(a) && !iszero(b) && return nothing # order and coefficient of first factorized term # In this case, since a[k]=0 for k>0, we can simplify to: # ordfact, cdivfact = 0, a/b[0] @@ -970,7 +994,8 @@ end @inline function div!(c::Taylor1{TaylorN{T}}, a::NumberNotSeries, b::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries} - iszero(a) && !iszero(b) && zero!(c, k) + zero!(c, k) + iszero(a) && !iszero(b) && return nothing # order and coefficient of first factorized term # In this case, since a[k]=0 for k>0, we can simplify to: # ordfact, cdivfact = 0, a/b[0] @@ -1141,14 +1166,14 @@ end """Division does not define a Taylor1 polynomial; order k=$(ordfact) => coeff[$(ordfact)]=$(cdivfact).""") ) + zero!(c, k) + if k == 0 # @inbounds c[0] = a[ordfact]/b[ordfact] @inbounds div!(c[0], a[ordfact], b[ordfact]) return nothing end - @inbounds zero!(c, k) - imin = max(0, k+ordfact-b.order) @inbounds mul!(c[k], c[imin], b[k+ordfact-imin]) @inbounds for i = imin+1:k-1 diff --git a/src/evaluate.jl b/src/evaluate.jl index 7cc23a67..6594c3fa 100644 --- a/src/evaluate.jl +++ b/src/evaluate.jl @@ -187,43 +187,42 @@ function _evaluate(a::HomogeneousPolynomial{T}, vals::NTuple) where {T} return suma end -function _evaluate(a::HomogeneousPolynomial{T}, vals::NTuple{N,<:TaylorN{T}}) where - {N,T<:NumberNotSeries} - # @assert length(vals) == get_numvars() - a.order == 0 && return a[1]*one(vals[1]) +function _evaluate!(res::TaylorN{T}, a::HomogeneousPolynomial{T}, vals::NTuple{N,<:TaylorN{T}}, + valscache::Vector{TaylorN{T}}, aux::TaylorN{T}) where {N,T<:NumberNotSeries} ct = coeff_table[a.order+1] - suma = TaylorN(zero(T), vals[1].order) - # - # vv = power_by_squaring.(vals, ct[1]) - vv = vals .^ ct[1] - tmp = zero(suma) - aux = one(suma) + for el in eachindex(valscache) + power_by_squaring!(valscache[el], vals[el], aux, ct[1][el]) + end for (i, a_coeff) in enumerate(a.coeffs) iszero(a_coeff) && continue - # @inbounds vv .= power_by_squaring.(vals, ct[i]) - vv .= vals .^ ct[i] - # tmp = prod( vv ) - for ord in eachindex(tmp) - @inbounds one!(aux, vv[1], ord) + # valscache .= vals .^ ct[i] + @inbounds for el in eachindex(valscache) + power_by_squaring!(valscache[el], vals[el], aux, ct[i][el]) end - for j in eachindex(vv) - for ord in eachindex(tmp) - zero!(tmp, ord) - @inbounds mul!(tmp, aux, vv[j], ord) - end - for ord in eachindex(tmp) - identity!(aux, tmp, ord) - end + # aux = one(valscache[1]) + for ord in eachindex(aux) + @inbounds one!(aux, valscache[1], ord) end - # suma += a_coeff * tmp - for ord in eachindex(tmp) - for ordQ in eachindex(tmp[ord]) - zero!(aux[ord], ordQ) - aux[ord][ordQ] = a_coeff * tmp[ord][ordQ] - suma[ord][ordQ] += aux[ord][ordQ] - end + for j in eachindex(valscache) + # aux *= valscache[j] + mul!(aux, valscache[j]) + end + # res += a_coeff * aux + for ord in eachindex(aux) + muladd!(res, a_coeff, aux, ord) end end + return nothing +end + +function _evaluate(a::HomogeneousPolynomial{T}, vals::NTuple{N,<:TaylorN{T}}) where + {N,T<:NumberNotSeries} + # @assert length(vals) == get_numvars() + a.order == 0 && return a[1]*one(vals[1]) + suma = TaylorN(zero(T), vals[1].order) + valscache = [zero(val) for val in vals] + aux = zero(suma) + _evaluate!(suma, a, vals, valscache, aux) return suma end @@ -319,6 +318,17 @@ _evaluate(a::TaylorN{T}, vals::NTuple, ::Val{true}) where {T<:NumberNotSeries} = _evaluate(a::TaylorN{T}, vals::NTuple, ::Val{false}) where {T<:Number} = sum( _evaluate(a, vals) ) +function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:TaylorN}, ::Val{false}) where {N,T<:Number} + R = promote_type(T, TS.numtype(vals[1])) + res = TaylorN(zero(R), vals[1].order) + valscache = [zero(val) for val in vals] + aux = zero(res) + @inbounds for homPol in eachindex(a) + _evaluate!(res, a[homPol], vals, valscache, aux) + end + return res +end + function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:Number}) where {N,T<:Number} R = promote_type(T, typeof(vals[1])) a_length = length(a) @@ -329,13 +339,20 @@ function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:Number}) where {N,T<:Number} return suma end -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) +function _evaluate!(res::Vector{TaylorN{T}}, a::TaylorN{T}, vals::NTuple{N,<:TaylorN}, + valscache::Vector{TaylorN{T}}, aux::TaylorN{T}) where {N,T<:Number} @inbounds for homPol in eachindex(a) - suma[homPol+1] = _evaluate(a[homPol], vals) + _evaluate!(res[homPol+1], a[homPol], vals, valscache, aux) end + return nothing +end + +function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:TaylorN}) where {N,T<:Number} + R = promote_type(T, TS.numtype(vals[1])) + suma = [TaylorN(zero(R), vals[1].order) for _ in eachindex(a)] + valscache = [zero(val) for val in vals] + aux = zero(suma[1]) + _evaluate!(suma, a, vals, valscache, aux) return suma end @@ -486,6 +503,24 @@ function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{TaylorN{T},1}, return nothing end +function evaluate!(a::TaylorN{T}, vals::NTuple{N,TaylorN{T}}, dest::TaylorN{T}, valscache::Vector{TaylorN{T}}, aux::TaylorN{T}) where {N,T<:Number} + @inbounds for homPol in eachindex(a) + _evaluate!(dest, a[homPol], vals, valscache, aux) + end + return nothing +end + +function evaluate!(a::AbstractArray{TaylorN{T}}, vals::NTuple{N,TaylorN{T}}, dest::AbstractArray{TaylorN{T}}) where {N,T<:Number} + # initialize evaluation cache + valscache = [zero(val) for val in vals] + aux = zero(dest[1]) + # loop over elements of `a` + for i in eachindex(a) + (!iszero(dest[i])) && zero!(dest[i]) + evaluate!(a[i], vals, dest[i], valscache, aux) + end + return nothing +end # In-place Horner methods, used when the result of an evaluation (substitution) # is Taylor1{} diff --git a/src/functions.jl b/src/functions.jl index 3df81f29..bb27b9af 100644 --- a/src/functions.jl +++ b/src/functions.jl @@ -354,12 +354,18 @@ for T in (:Taylor1, :TaylorN) return nothing end - @inline function one!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number} - zero!(c, k) - if k == 0 - @inbounds c[0] = one(a[0]) + if $T == Taylor1 + @inline function one!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number} + zero!(c, k) + (k == 0) && (@inbounds c[0] = one(constant_term(a))) + return nothing + end + else + @inline function one!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number} + zero!(c, k) + (k == 0) && (@inbounds c[0][1] = one(constant_term(a))) + return nothing end - return nothing end @inline function abs!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number} diff --git a/src/power.jl b/src/power.jl index 30bbc0f3..67ca6726 100644 --- a/src/power.jl +++ b/src/power.jl @@ -62,8 +62,15 @@ end # this method assumes `y`, `x` and `aux` are of same order # TODO: add power_by_squaring! method for HomogeneousPolynomial and mixtures for T in (:Taylor1, :TaylorN) - @eval function power_by_squaring!(y::$T{T}, x::$T{T}, aux::$T{T}, + @eval @inline function power_by_squaring_0!(y::$T{T}, x::$T{T}) where {T<:NumberNotSeries} + for k in eachindex(y) + one!(y, x, k) + end + return nothing + end + @eval @inline function power_by_squaring!(y::$T{T}, x::$T{T}, aux::$T{T}, p::Integer) where {T<:NumberNotSeries} + (p == 0) && return power_by_squaring_0!(y, x) t = trailing_zeros(p) + 1 p >>= t # aux = x @@ -272,7 +279,7 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`. isinteger(r) && r > 0 && (k > r*findlast(a)) && return nothing if k == lnull - @inbounds c[k] = (a[l0])^r + @inbounds c[k] = (a[l0])^float(r) return nothing end diff --git a/test/manyvariables.jl b/test/manyvariables.jl index 32b4e88b..be842a61 100644 --- a/test/manyvariables.jl +++ b/test/manyvariables.jl @@ -821,6 +821,33 @@ end @test float(TaylorN{Complex{Rational}}) == float(TaylorN{Complex{Float64}}) end + @testset "Test evaluate! method for arrays of TaylorN" begin + x = set_variables("x", order=2, numvars=4) + function radntn!(y) + for k in eachindex(y) + for l in eachindex(y[k]) + y[k][l] = randn() + end + end + nothing + end + y = zero(x[1]) + radntn!(y) + n = 10 + v = [zero(x[1]) for _ in 1:n] + r = [zero(x[1]) for _ in 1:n] # output vector + radntn!.(v) + x1 = randn(4) .+ x + # warmup + TaylorSeries.evaluate!(v, (x1...,), r) + # call twice to make sure `r` is reset on second call + TaylorSeries.evaluate!(v, (x1...,), r) + r2 = evaluate.(v, Ref(x1)) + @test r == r2 + @test iszero(norm(r-r2, Inf)) + + end + end @testset "Integrate for several variables" begin