Avoid some allocations in mutating methods

Project.toml

@@ -9,6 +9,12 @@ Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
+[weakdeps]
+IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+
+[extensions]
+TaylorSeriesIAExt = "IntervalArithmetic"
+
 [compat]
 Aqua = "0.8"
 IntervalArithmetic = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20"
@@ -19,9 +25,6 @@ SparseArrays = "<0.0.1, 1"
 Test = "<0.0.1, 1"
 julia = "1"
 
-[extensions]
-TaylorSeriesIAExt = "IntervalArithmetic"
-
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
@@ -31,6 +34,3 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["IntervalArithmetic", "LinearAlgebra", "SparseArrays", "Test", "Aqua"]
-
-[weakdeps]
-IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"

src/TaylorSeries.jl

@@ -16,7 +16,6 @@ see also [`HomogeneousPolynomial`](@ref).
 """
 module TaylorSeries
 
-
 using SparseArrays: SparseMatrixCSC
 using Markdown
 

src/arithmetic.jl

@@ -482,7 +482,7 @@ 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)
+    zero!(res, ordT)
     for k in 0:ordT
         @inbounds for ordQ in eachindex(a[ordT])
             mul!(res[ordT], a[k], b[ordT-k], ordQ)
@@ -665,6 +665,19 @@ function /(a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}}) where {T<:NumberNotSe
     return res
 end
 
+function /(a::NumberNotSeries, b::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
+    iszero(a) && !iszero(b) && return zero(a)
+
+    # 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]
+
+    res = Taylor1(a/b[0], b.order)
+    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
@@ -731,8 +744,25 @@ end
     return nothing
 end
 
-div!(v::Taylor1, b::NumberNotSeries, a::Taylor1, k::Int) =
-    div!(v::Taylor1, Taylor1(b, a.order), a, k)
+@inline function div!(c::Taylor1, a::NumberNotSeries,
+        b::Taylor1, k::Int)
+    iszero(a) && !iszero(b) && zero!(c, k)
+    # 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]
+    if k == 0
+        @inbounds c[0] = a/b[0]
+        return nothing
+    end
+
+    imin = max(0, k-b.order)
+    @inbounds c[k] = c[imin] * b[k-imin]
+    @inbounds for i = imin+1:k-1
+        c[k] += c[i] * b[k-i]
+    end
+    @inbounds c[k] = -c[k]/b[0]
+    return nothing
+end
 
 # 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)
@@ -758,7 +788,7 @@ end
         @inbounds res[0] = cdivfact
         return nothing
     end
-    zero!(res, a, ordT)
+    zero!(res, ordT)
     imin = max(0, ordT+ordfact-b.order)
     aux = TaylorN(zero(a[ordT][0][1]), a[ordT].order)
     for k in imin:ordT-1
@@ -783,6 +813,48 @@ end
     return nothing
 end
 
+@inline function mul!(res::Taylor1{TaylorN{T}}, a::NumberNotSeries,
+        b::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries}
+    for l in eachindex(b[k])
+        for m in eachindex(b[k][l])
+            res[k][l][m] = a*b[k][l][m]
+        end
+    end
+    return nothing
+end
+
+mul!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
+    b::NumberNotSeries, k::Int) where {T<:NumberNotSeries} = mul!(res, b, a, k)
+
+@inline function div!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
+        b::NumberNotSeries, k::Int) where {T<:NumberNotSeries}
+    for l in eachindex(a[k])
+        for m in eachindex(a[k][l])
+            res[k][l][m] = a[k][l][m]/b
+        end
+    end
+    return nothing
+end
+
+# ### TODO: use non-allocating div!(::TaylorN,::TaylorN,::TaylorN)
+# @inline function div!(c::Taylor1{TaylorN{T}}, a::NumberNotSeries,
+#         b::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries}
+#     # 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]
+#     if k == 0
+#         @inbounds c[0] = a/b[0]
+#         return nothing
+#     end
+
+#     imin = max(0, k-b.order)
+#     @inbounds c[k] = c[imin] * b[k-imin]
+#     @inbounds for i = imin+1:k-1
+#         c[k] += c[i] * b[k-i]
+#     end
+#     @inbounds c[k] = -c[k]/b[0]
+#     return nothing
+# end