Skip to content

Commit 1575b39

Browse files
PerezHzPerezHz
committed
WIP: add evaluate! method for TaylorNs
1 parent 34d0fea commit 1575b39

File tree

5 files changed

+130
-60
lines changed

5 files changed

+130
-60
lines changed

‎Project.toml‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@ repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
44
version = "0.17.8"
55

66
[deps]
7+
InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
78
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
89
Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
910
Requires = "ae029012-a4dd-5104-9daa-d747884805df"
@@ -12,14 +13,14 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
1213
[weakdeps]
1314
IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
1415
JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
15-
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
1616
RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
17+
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
1718

1819
[extensions]
1920
TaylorSeriesIAExt = "IntervalArithmetic"
2021
TaylorSeriesJLD2Ext = "JLD2"
21-
TaylorSeriesSAExt = "StaticArrays"
2222
TaylorSeriesRATExt = "RecursiveArrayTools"
23+
TaylorSeriesSAExt = "StaticArrays"
2324

2425
[compat]
2526
Aqua = "0.8"

‎src/arithmetic.jl‎

Lines changed: 40 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -162,9 +162,7 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
162162
end
163163

164164
for T in (:Taylor1, :TaylorN)
165-
@eval function zero(a::$T)
166-
return $T(zero.(a.coeffs))
167-
end
165+
@eval zero(a::$T) = $T(zero.(a.coeffs))
168166
@eval function one(a::$T)
169167
b = zero(a)
170168
b[0] = one(b[0])
@@ -539,25 +537,50 @@ for T in (:Taylor1, :TaylorN)
539537
return nothing
540538
end
541539

542-
@eval @inline function mul!(v::$T, a::$T, b::NumberNotSeries, k::Int)
540+
@eval begin
543541
if $T == Taylor1
544-
@inbounds v[k] = a[k] * b
545-
else
546-
@inbounds for i in eachindex(v[k])
547-
v[k][i] = a[k][i] * b
542+
@inline function mul!(v::$T, a::$T, b::NumberNotSeries, k::Int)
543+
@inbounds v[k] = a[k] * b
544+
return nothing
545+
end
546+
@inline function mul!(v::$T, a::NumberNotSeries, b::$T, k::Int)
547+
@inbounds v[k] = a * b[k]
548+
return nothing
549+
end
550+
@inline function muladd!(v::$T, a::$T, b::NumberNotSeries, k::Int)
551+
@inbounds v[k] += a[k] * b
552+
return nothing
553+
end
554+
@inline function muladd!(v::$T, a::NumberNotSeries, b::$T, k::Int)
555+
@inbounds v[k] += a * b[k]
556+
return nothing
548557
end
549-
end
550-
return nothing
551-
end
552-
@eval @inline function mul!(v::$T, a::NumberNotSeries, b::$T, k::Int)
553-
if $T == Taylor1
554-
@inbounds v[k] = a * b[k]
555558
else
556-
@inbounds for i in eachindex(v[k])
557-
v[k][i] = a * b[k][i]
559+
@inline function mul!(v::$T, a::$T, b::NumberNotSeries, k::Int)
560+
@inbounds for i in eachindex(v[k])
561+
v[k][i] = a[k][i] * b
562+
end
563+
return nothing
564+
end
565+
@inline function mul!(v::$T, a::NumberNotSeries, b::$T, k::Int)
566+
@inbounds for i in eachindex(v[k])
567+
v[k][i] = a * b[k][i]
568+
end
569+
return nothing
570+
end
571+
@inline function muladd!(v::$T, a::$T, b::NumberNotSeries, k::Int)
572+
@inbounds for i in eachindex(v[k])
573+
v[k][i] += a[k][i] * b
574+
end
575+
return nothing
576+
end
577+
@inline function muladd!(v::$T, a::NumberNotSeries, b::$T, k::Int)
578+
@inbounds for i in eachindex(v[k])
579+
v[k][i] += a * b[k][i]
580+
end
581+
return nothing
558582
end
559583
end
560-
return nothing
561584
end
562585

563586
@eval @inline function mul!(v::$T, a::$T, b::NumberNotSeries)

‎src/evaluate.jl‎

Lines changed: 68 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -187,43 +187,42 @@ function _evaluate(a::HomogeneousPolynomial{T}, vals::NTuple) where {T}
187187
return suma
188188
end
189189

190-
function _evaluate(a::HomogeneousPolynomial{T}, vals::NTuple{N,<:TaylorN{T}}) where
191-
{N,T<:NumberNotSeries}
192-
# @assert length(vals) == get_numvars()
193-
a.order == 0 && return a[1]*one(vals[1])
190+
function _evaluate!(res::TaylorN{T}, a::HomogeneousPolynomial{T}, vals::NTuple{N,<:TaylorN{T}},
191+
valscache::Vector{TaylorN{T}}, aux::TaylorN{T}) where {N,T<:NumberNotSeries}
194192
ct = coeff_table[a.order+1]
195-
suma = TaylorN(zero(T), vals[1].order)
196-
#
197-
# vv = power_by_squaring.(vals, ct[1])
198-
vv = vals .^ ct[1]
199-
tmp = zero(suma)
200-
aux = one(suma)
193+
for el in eachindex(valscache)
194+
power_by_squaring!(valscache[el], vals[el], aux, ct[1][el])
195+
end
201196
for (i, a_coeff) in enumerate(a.coeffs)
202197
iszero(a_coeff) && continue
203-
# @inbounds vv .= power_by_squaring.(vals, ct[i])
204-
vv .= vals .^ ct[i]
205-
# tmp = prod( vv )
206-
for ord in eachindex(tmp)
207-
@inbounds one!(aux, vv[1], ord)
198+
# valscache .= vals .^ ct[i]
199+
@inbounds for el in eachindex(valscache)
200+
power_by_squaring!(valscache[el], vals[el], aux, ct[i][el])
208201
end
209-
for j in eachindex(vv)
210-
for ord in eachindex(tmp)
211-
zero!(tmp, ord)
212-
@inbounds mul!(tmp, aux, vv[j], ord)
213-
end
214-
for ord in eachindex(tmp)
215-
identity!(aux, tmp, ord)
216-
end
202+
# aux = one(valscache[1])
203+
for ord in eachindex(aux)
204+
@inbounds one!(aux, valscache[1], ord)
217205
end
218-
# suma += a_coeff * tmp
219-
for ord in eachindex(tmp)
220-
for ordQ in eachindex(tmp[ord])
221-
zero!(aux[ord], ordQ)
222-
aux[ord][ordQ] = a_coeff * tmp[ord][ordQ]
223-
suma[ord][ordQ] += aux[ord][ordQ]
224-
end
206+
for j in eachindex(valscache)
207+
# aux *= valscache[j]
208+
mul!(aux, valscache[j])
209+
end
210+
# res += a_coeff * aux
211+
for ord in eachindex(aux)
212+
muladd!(res, a_coeff, aux, ord)
225213
end
226214
end
215+
return nothing
216+
end
217+
218+
function _evaluate(a::HomogeneousPolynomial{T}, vals::NTuple{N,<:TaylorN{T}}) where
219+
{N,T<:NumberNotSeries}
220+
# @assert length(vals) == get_numvars()
221+
a.order == 0 && return a[1]*one(vals[1])
222+
suma = TaylorN(zero(T), vals[1].order)
223+
valscache = [zero(val) for val in vals]
224+
aux = zero(suma)
225+
_evaluate!(suma, a, vals, valscache, aux)
227226
return suma
228227
end
229228

@@ -319,6 +318,17 @@ _evaluate(a::TaylorN{T}, vals::NTuple, ::Val{true}) where {T<:NumberNotSeries} =
319318
_evaluate(a::TaylorN{T}, vals::NTuple, ::Val{false}) where {T<:Number} =
320319
sum( _evaluate(a, vals) )
321320

321+
function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:TaylorN}, ::Val{false}) where {N,T<:Number}
322+
R = promote_type(T, TS.numtype(vals[1]))
323+
res = TaylorN(zero(R), vals[1].order)
324+
valscache = [zero(val) for val in vals]
325+
aux = zero(res)
326+
@inbounds for homPol in eachindex(a)
327+
_evaluate!(res, a[homPol], vals, valscache, aux)
328+
end
329+
return res
330+
end
331+
322332
function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:Number}) where {N,T<:Number}
323333
R = promote_type(T, typeof(vals[1]))
324334
a_length = length(a)
@@ -329,13 +339,20 @@ function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:Number}) where {N,T<:Number}
329339
return suma
330340
end
331341

332-
function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:TaylorN}) where {N,T<:Number}
333-
R = TaylorN{promote_type(T, TS.numtype(vals[1]))}
334-
a_length = length(a)
335-
suma = zeros(R, a_length)
342+
function _evaluate!(res::Vector{TaylorN{T}}, a::TaylorN{T}, vals::NTuple{N,<:TaylorN},
343+
valscache::Vector{TaylorN{T}}, aux::TaylorN{T}) where {N,T<:Number}
336344
@inbounds for homPol in eachindex(a)
337-
suma[homPol+1] = _evaluate(a[homPol], vals)
345+
_evaluate!(res[homPol+1], a[homPol], vals, valscache, aux)
338346
end
347+
return nothing
348+
end
349+
350+
function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:TaylorN}) where {N,T<:Number}
351+
R = promote_type(T, TS.numtype(vals[1]))
352+
suma = [TaylorN(zero(R), vals[1].order) for _ in eachindex(a)]
353+
valscache = [zero(val) for val in vals]
354+
aux = zero(suma[1])
355+
_evaluate!(suma, a, vals, valscache, aux)
339356
return suma
340357
end
341358

@@ -486,6 +503,22 @@ function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{TaylorN{T},1},
486503
return nothing
487504
end
488505

506+
function evaluate!(a::TaylorN{T}, vals::NTuple{N,TaylorN{T}}, dest::TaylorN{T}, valscache::Vector{TaylorN{T}}, aux::TaylorN{T}) where {N,T<:Number}
507+
@inbounds for homPol in eachindex(a)
508+
_evaluate!(dest, a[homPol], vals, valscache, aux)
509+
end
510+
return nothing
511+
end
512+
513+
function evaluate!(a::AbstractArray{TaylorN{T}}, vals::NTuple{N,TaylorN{T}}, dest::AbstractArray{TaylorN{T}}) where {N,T<:Number}
514+
valscache = [zero(val) for val in vals]
515+
aux = zero(dest[1])
516+
for i in eachindex(a)
517+
(!iszero(a[i])) && zero!(dest[i])
518+
evaluate!(a[i], vals, dest[i], valscache, aux)
519+
end
520+
return nothing
521+
end
489522

490523
# In-place Horner methods, used when the result of an evaluation (substitution)
491524
# is Taylor1{}

‎src/functions.jl‎

Lines changed: 11 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -354,12 +354,18 @@ for T in (:Taylor1, :TaylorN)
354354
return nothing
355355
end
356356

357-
@inline function one!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
358-
zero!(c, k)
359-
if k == 0
360-
@inbounds c[0] = one(a[0])
357+
if $T == Taylor1
358+
@inline function one!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
359+
zero!(c, k)
360+
(k == 0) && (@inbounds c[0] = one(constant_term(a)))
361+
return nothing
362+
end
363+
else
364+
@inline function one!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
365+
zero!(c, k)
366+
(k == 0) && (@inbounds c[0][1] = one(constant_term(a)))
367+
return nothing
361368
end
362-
return nothing
363369
end
364370

365371
@inline function abs!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}

‎src/power.jl‎

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -62,8 +62,15 @@ end
6262
# this method assumes `y`, `x` and `aux` are of same order
6363
# TODO: add power_by_squaring! method for HomogeneousPolynomial
6464
for T in (:Taylor1, :TaylorN)
65-
@eval function power_by_squaring!(y::$T{T}, x::$T{T}, aux::$T{T},
65+
@eval @inline function power_by_squaring_0!(y::$T{T}, x::$T{T}) where {T<:NumberNotSeries}
66+
for k in eachindex(y)
67+
one!(y, x, k)
68+
end
69+
return nothing
70+
end
71+
@eval @inline function power_by_squaring!(y::$T{T}, x::$T{T}, aux::$T{T},
6672
p::Integer) where {T<:NumberNotSeries}
73+
(p == 0) && return power_by_squaring_0!(y, x)
6774
t = trailing_zeros(p) + 1
6875
p >>= t
6976
# aux = x

0 commit comments

Comments
 (0)