Further allocation reduction in power computations (#361)

* working version with oop sqr!

* Update Project.toml

* Update src/TaylorSeries.jl

* Update `zero(x)`; add in-place mul! for TaylorN

* Get rid of most allocations in pow! for Taylor1{TaylorN{T}; add power_by_squaring!; add in-place sqr! method

* Remove extra empty line

* Add power_by_squaring methods to avoid method ambiguities (detected by Aqua)

* Remove duplicated return

* Add inbounds

* Address review by @lbenet

* Add extra arg to pow!

* Update TaylorSeries IA extension

* Update pow! and add fix suggested by @lbenet

* Middle-of-the-road approach to suggestion by @lbenet

* Handle pow! cases by dispatch-by-value

* Remove unneeded deepcopy in setindex! method

* Add `power_by_squaring(x, ::Val{3})` methods and add tests

* Update comments

* Bump patch version

* Switch back from dispatch-by-value to if/else for pow!

* Revert change in setindex! overload for Taylor1

* Add test for nested Taylor1s setindex! method

* De-bump patch version

* Another approach to suggestion by @lbenet

* Fix typo

* Fix another typo

Author: PerezHz

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
-        TS.pow!(c, aa, r, k)
+        TS.pow!(c, aa, c, 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
-        TS.pow!(c, aa, r, ord)
+        TS.pow!(c, aa, c, r, ord)
     end
 
     return c

src/arithmetic.jl

@@ -162,7 +162,9 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
 end
 
 for T in (:Taylor1, :TaylorN)
-    @eval zero(a::$T) = $T(zero.(a.coeffs))
+    @eval function zero(a::$T)
+        return $T(zero.(a.coeffs))
+    end
     @eval function one(a::$T)
         b = zero(a)
         b[0] = one(b[0])
@@ -572,6 +574,42 @@ for T in (:Taylor1, :TaylorN)
     end
 end
 
+# in-place product: `a` <- `a*b`
+# this method computes the product `a*b` and saves it back into `a`
+# assumes `a` and `b` are of same order
+function mul!(a::TaylorN{T}, b::TaylorN{T}) where {T<:Number}
+    @inbounds for k in reverse(eachindex(a))
+        mul!(a, a, b[0][1], k)
+        for l in 1:k
+            mul!(a[k], a[k-l], b[l])
+        end
+    end
+    return nothing
+end
+function mul!(a::Taylor1{T}, b::Taylor1{T}) where {T<:Number}
+    @inbounds for k in reverse(eachindex(a))
+        # a[k] <- a[k]*b[0]
+        mul!(a, a, b[0], k)
+        for l in 1:k
+            # a[k] <- a[k] + a[k-l] * b[l]
+            a[k] += a[k-l] * b[l]
+        end
+    end
+    return nothing
+end
+function mul!(a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
+    @inbounds for k in reverse(eachindex(a))
+        mul!(a, a, b[0], k)
+        for l in 1:k
+            # a[k] += a[k-l] * b[l]
+            for m in eachindex(a[k])
+                mul!(a[k], a[k-l], b[l], m)
+            end
+        end
+    end
+    return nothing
+end
+
 function mul!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}},
         ordT::Int) where {T<:NumberNotSeries}
     # Sanity

src/auxiliary.jl

@@ -17,10 +17,10 @@ If the length of `coeffs` is smaller than `order+1`, it resizes
 function resize_coeffs1!(coeffs::Array{T,1}, order::Int) where {T<:Number}
     lencoef = length(coeffs)
     resize!(coeffs, order+1)
+    c1 = coeffs[1]
     if order > lencoef-1
-        z = zero(coeffs[1])
         @simd for ord in lencoef+1:order+1
-            @inbounds coeffs[ord] = z
+            @inbounds coeffs[ord] = zero(c1)
         end
     end
     return nothing
@@ -39,9 +39,9 @@ function resize_coeffsHP!(coeffs::Array{T,1}, order::Int) where {T<:Number}
     @assert order ≤ get_order() && lencoef ≤ num_coeffs
     num_coeffs == lencoef && return nothing
     resize!(coeffs, num_coeffs)
-    z = zero(coeffs[1])
+    c1 = coeffs[1]
     @simd for ord in lencoef+1:num_coeffs
-        @inbounds coeffs[ord] = z
+        @inbounds coeffs[ord] = zero(c1)
     end
     return nothing
 end

src/constructors.jl

@@ -156,7 +156,7 @@ struct TaylorN{T<:Number} <: AbstractSeries{T}
     order   :: Int
 
     function TaylorN{T}(v::Array{HomogeneousPolynomial{T},1}, order::Int) where T<:Number
-        coeffs = zeros(HomogeneousPolynomial{T}, order)
+        coeffs = isempty(v) ? zeros(HomogeneousPolynomial{T}, order) : zeros(v[1], order)
         @inbounds for i in eachindex(v)
             ord = v[i].order
             if ord ≤ order

src/dictmutfunct.jl

@@ -60,7 +60,7 @@ const _dict_binary_ops = Dict(
     :- => (:subst!, (:_res, :_arg1, :_arg2, :_k), :(_res = _arg1 - _arg2)),
     :* => (:mul!, (:_res, :_arg1, :_arg2, :_k), :(_res = _arg1 * _arg2)),
     :/ => (:div!, (:_res, :_arg1, :_arg2, :_k), :(_res = _arg1 / _arg2)),
-    :^ => (:pow!, (:_res, :_arg1, :_arg2, :_k), :(_res = _arg1 ^ float(_arg2))),
+    :^ => (:pow!, (:_res, :_arg1, :_aux, :_arg2, :_k), :(_res = _arg1 ^ float(_arg2)), :(_aux = zero(_arg1))),
 );
 
 """
@@ -196,4 +196,3 @@ internal mutating functions.
 Evaluating the entries generates symbols that represent
 the actual calls to the internal mutating functions.
 """ _dict_binary_calls
-

src/evaluate.jl

@@ -413,7 +413,7 @@ function _evaluate!(suma::TaylorN{T}, a::HomogeneousPolynomial{T}, ind::Int,
         else
             for ordQ in eachindex(val)
                 zero!(vv, ordQ)
-                pow!(vv, val, vpow, ordQ)
+                pow!(vv, val, vv, vpow, ordQ)
             end
         end
         for ordQ in eachindex(suma)

src/functions.jl

@@ -303,6 +303,13 @@ end
     return nothing
 end
 
+@inline function zero!(a::TaylorN{Taylor1{T}}) where {T<:NumberNotSeries}
+    for k in eachindex(a)
+        zero!(a, k)
+    end
+    return nothing
+end
+
 @inline function one!(c::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries}
     zero!(c, k)
     if k == 0
@@ -316,6 +323,13 @@ end
     return nothing
 end
 
+@inline function identity!(c::HomogeneousPolynomial{Taylor1{T}}, a::HomogeneousPolynomial{Taylor1{T}}, k::Int) where {T<:NumberNotSeries}
+    @inbounds for l in eachindex(c[k])
+        identity!(c[k], a[k], l)
+    end
+    return nothing
+end
+
 for T in (:Taylor1, :TaylorN)
     @eval begin
         @inline function identity!(c::$T{T}, a::$T{T}, k::Int) where {T<:NumberNotSeries}