Reduce memory allocation when squaring a nested Taylor1 (#388)

* Pass aux as argument to the nested Taylor1s-specialized mutating method of sqr

* Add tests for the nested Taylor1s-specialized mutating version of pow

* Bump nested Taylor1s mutating methods of subst, div and pow to NumberNotSeriesN as suggested by @lbenet

* Include compatibility with IA v0.23

* Bump patch version

Author: LuEdRaMo

Project.toml

@@ -1,7 +1,7 @@
 name = "TaylorSeries"
 uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
-version = "0.20.6"
+version = "0.20.7"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -22,7 +22,7 @@ TaylorSeriesSAExt = "StaticArrays"
 
 [compat]
 Aqua = "0.8"
-IntervalArithmetic = "1"
+IntervalArithmetic = "0.23, 1"
 JLD2 = "0.5"
 LinearAlgebra = "<0.0.1, 1"
 Markdown = "<0.0.1, 1"

src/arithmetic.jl

@@ -473,7 +473,7 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
     end
 end
 
-function subst!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, k::Int) where {T <: TS.NumberNotSeries}
+function subst!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, k::Int) where {T <: NumberNotSeriesN}
     @inbounds for i in eachindex(v[k])
         v[k][i] = -a[k][i]
     end
@@ -1236,7 +1236,7 @@ function div!(c::Taylor1, a::NumberNotSeries, b::Taylor1)
 end
 
 @inline function div!(c::Taylor1{Taylor1{T}}, a::NumberNotSeries,
-        b::Taylor1{Taylor1{T}}, k::Int) where {T<:NumberNotSeries}
+        b::Taylor1{Taylor1{T}}, k::Int) where {T<:NumberNotSeriesN}
     zero!(c, k)
     iszero(a) && !iszero(b) && return nothing
     # order and coefficient of first factorized term

src/power.jl

@@ -321,11 +321,11 @@ end
     return nothing
 end
 
-@inline function pow!(c::Taylor1{T}, a::Taylor1{T}, aux::Taylor1{T},
+@inline function pow!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, aux::Taylor1{Taylor1{T}},
         r::S, k::Int) where {T<:NumberNotSeriesN, S<:Real}
     (r == 0) && return one!(c, a, k)
     (r == 1) && return identity!(c, a, k)
-    (r == 2) && return sqr!(c, a, k)
+    (r == 2) && return sqr!(c, a, aux[k], k)
     (r == 0.5) && return sqrt!(c, a, k)
     # Sanity
     zero!(c, k)
@@ -553,7 +553,7 @@ end
     return nothing
 end
 
-@inline function sqr!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
+@inline function sqr!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, aux::Taylor1{T},
         k::Int) where {T<:NumberNotSeriesN}
     if k == 0
         sqr_orderzero!(c, a)
@@ -564,7 +564,7 @@ end
     # Recursion formula
     kodd = k%2
     kend = (k - 2 + kodd) >> 1
-    aux = zero(c[k])
+    # aux = zero(c[k])
     @inbounds for i = 0:kend
         for j in eachindex(a[k])
             # c[k] += 2 * a[i] * a[k-i]

test/mixtures.jl

@@ -671,6 +671,21 @@ end
 
             end
         end
+
+        for r in (0, 1, 2, 0.5)
+            local inorder = 6
+            local outorder = 25
+
+            x = Taylor1([Taylor1(rand(inorder+1), inorder) for _ in 1:outorder+1], outorder)
+            y = zero(x)
+            z = zero(x)
+
+            w = x^r
+            for k in eachindex(z)
+                TS.pow!(z, x, y, r, k)
+            end
+            @test norm(z - w, Inf) == 0.0
+        end
     end
 
     # Back to default