Further allocation reduction in power computations

src/arithmetic.jl

@@ -162,7 +162,10 @@ 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))
+        return za
+    end
     @eval function one(a::$T)
         b = zero(a)
         b[0] = one(b[0])
@@ -572,6 +575,19 @@ 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}
+    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!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}},
         ordT::Int) where {T<:NumberNotSeries}
     # Sanity

src/auxiliary.jl

@@ -18,9 +18,8 @@ function resize_coeffs1!(coeffs::Array{T,1}, order::Int) where {T<:Number}
     lencoef = length(coeffs)
     resize!(coeffs, order+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(coeffs[1])
         end
     end
     return nothing
@@ -39,9 +38,8 @@ 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])
     @simd for ord in lencoef+1:num_coeffs
-        @inbounds coeffs[ord] = z
+        @inbounds coeffs[ord] = zero(coeffs[1])
     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/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}

src/power.jl

@@ -58,37 +58,95 @@ end
 
 ^(a::Taylor1{TaylorN{T}}, r::Rational) where {T<:NumberNotSeries} = a^(r.num/r.den)
 
+# in-place form of power_by_squaring
+# this method assumes `y`, `x` and `aux1` are of same order
+function power_by_squaring!(y::TaylorN{T}, x::TaylorN{T}, aux1::TaylorN{T}, p::Integer) where {T<:NumberNotSeries}
+    t = trailing_zeros(p) + 1
+    p >>= t
+    # aux1 = x
+    for k in eachindex(aux1)
+        identity!(aux1, x, k)
+    end
+    while (t -= 1) > 0
+        # aux1 = square(aux1)
+        for k in reverse(eachindex(aux1))
+            sqr!(aux1, k)
+        end
+    end
+    # y = aux1
+    for k in eachindex(y)
+        identity!(y, aux1, k)
+    end
+    while p > 0
+        t = trailing_zeros(p) + 1
+        p >>= t
+        while (t -= 1) ≥ 0
+            # aux1 = square(aux1)
+            for k in reverse(eachindex(aux1))
+                sqr!(aux1, k)
+            end
+        end
+        # y = y * aux1
+        mul!(y, aux1)
+    end
+    return nothing
+end
 
 
 # power_by_squaring; slightly modified from base/intfuncs.jl
 # Licensed under MIT "Expat"
 for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
-    @eval function power_by_squaring(x::$T, p::Integer)
-        if p == 1
-            return copy(x)
-        elseif p == 0
-            return one(x)
-        elseif p == 2
-            return square(x)
-        end
-        t = trailing_zeros(p) + 1
-        p >>= t
-        while (t -= 1) > 0
-            x = square(x)
-        end
-        y = x
-        while p > 0
+    @eval power_by_squaring(x::$T, p::Integer) = power_by_squaring(x, Val(p))
+    @eval power_by_squaring(x::$T, ::Val{0}) = one(x)
+    @eval power_by_squaring(x::$T, ::Val{1}) = copy(x)
+    @eval power_by_squaring(x::$T, ::Val{2}) = square(x)
+    @eval function power_by_squaring(x::$T, ::Val{P}) where P
+        p = P # copy static parameter `P` into local variable `p`
+        if $T == TaylorN
+            y = zero(x)
+            aux1 = zero(x)
+            power_by_squaring!(y, x, aux1, p)
+        else
             t = trailing_zeros(p) + 1
             p >>= t
-            while (t -= 1) ≥ 0
+            while (t -= 1) > 0
                 x = square(x)
             end
-            y *= x
+            y = x
+            while p > 0
+                t = trailing_zeros(p) + 1
+                p >>= t
+                while (t -= 1) ≥ 0
+                    x = square(x)
+                end
+                y *= x
+            end
         end
         return y
     end
 end
 
+power_by_squaring(x::TaylorN{Taylor1{T}}, ::Val{0}) where {T<:NumberNotSeries} = one(x)
+power_by_squaring(x::TaylorN{Taylor1{T}}, ::Val{1}) where {T<:NumberNotSeries} = copy(x)
+power_by_squaring(x::TaylorN{Taylor1{T}}, ::Val{2}) where {T<:NumberNotSeries} = square(x)
+function power_by_squaring(x::TaylorN{Taylor1{T}}, ::Val{P}) where {P, T<:NumberNotSeries}
+    p = P # copy static parameter `P` into local variable `p`
+    t = trailing_zeros(p) + 1
+    p >>= t
+    while (t -= 1) > 0
+        x = square(x)
+    end
+    y = x
+    while p > 0
+        t = trailing_zeros(p) + 1
+        p >>= t
+        while (t -= 1) ≥ 0
+            x = square(x)
+        end
+        y *= x
+    end
+    return y
+end
 
 ## Real power ##
 function ^(a::Taylor1{T}, r::S) where {T<:Number, S<:Real}
@@ -300,7 +358,8 @@ end
     if ordT == lnull
         if isinteger(r)
             # TODO: get rid of allocations here
-            res[ordT] = a[l0]^round(Int,r) # uses power_by_squaring
+            aux1 = deepcopy(res[ordT])
+            power_by_squaring!(res[ordT], a[l0], aux1, round(Int,r))
             return nothing
         end
 
@@ -339,7 +398,7 @@ Return `a^2`; see [`TaylorSeries.sqr!`](@ref).
 
 for T in (:Taylor1, :TaylorN)
     @eval function square(a::$T)
-        c = $T( zero(constant_term(a)), a.order)
+        c = zero(a)
         for k in eachindex(a)
             sqr!(c, a, k)
         end
@@ -447,6 +506,37 @@ for T = (:Taylor1, :TaylorN)
 
             return nothing
         end
+
+        # in-place squaring: given `c`, compute expansion of `c^2` and save back into `c`
+        @inline function sqr!(c::$T{T}, k::Int) where {T<:NumberNotSeries}
+            if k == 0
+                sqr_orderzero!(c, c)
+                return nothing
+            end
+
+            # Recursion formula
+            kodd = k%2
+            kend = (k - 2 + kodd) >> 1
+            if $T == Taylor1
+                (kodd == 0) && ( @inbounds c[k] = (c[k >> 1]^2)/2 )
+                c[k] = c[0] * c[k]
+                @inbounds for i = 1:kend
+                    c[k] += c[i] * c[k-i]
+                end
+                @inbounds c[k] = 2 * c[k]
+            else
+                (kend ≥ 0) && mul!(c, c[0][1], c, k)
+                @inbounds for i = 1:kend
+                    mul!(c[k], c[i], c[k-i])
+                end
+                @inbounds mul!(c, 2, c, k)
+                if (kodd == 0)
+                    accsqr!(c[k], c[k >> 1])
+                end
+            end
+
+            return nothing
+        end
     end
 end