From dbdbcc111d75bbf38e74a4a5f5aea23f9ae8cb12 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Mon, 8 Jul 2024 20:17:02 +0200
Subject: [PATCH 01/26] working version with oop sqr!

---
 Project.toml        |   2 +
 src/TaylorSeries.jl |   2 +
 src/arithmetic.jl   |  40 ++++++++++++-
 src/auxiliary.jl    |   6 +-
 src/constructors.jl |   2 +-
 src/functions.jl    |  14 +++++
 src/power.jl        | 140 ++++++++++++++++++++++++++++++++++++++------
 test/runtests.jl    |  18 +++---
 8 files changed, 190 insertions(+), 34 deletions(-)

diff --git a/Project.toml b/Project.toml
index e25d96af..eed44623 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,6 +8,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
 [weakdeps]
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
@@ -32,6 +33,7 @@ Requires = "0.5.2, 1"
 SparseArrays = "<0.0.1, 1"
 StaticArrays = "1"
 Test = "<0.0.1, 1"
+InteractiveUtils = "<0.0.1, 1"
 julia = "1"
 
 [extras]
diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl
index 2e508e66..82bb84cf 100644
--- a/src/TaylorSeries.jl
+++ b/src/TaylorSeries.jl
@@ -28,6 +28,8 @@ using LinearAlgebra: norm, mul!,
     lu, lu!, LinearAlgebra.lutype, LinearAlgebra.copy_oftype,
     LinearAlgebra.issuccess
 
+using InteractiveUtils
+
 if VERSION >= v"1.7.0-DEV.1188"
     using LinearAlgebra: NoPivot, RowMaximum
 end
diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 2af3e8eb..abfc439e 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -162,7 +162,12 @@ 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))
+        # za = deepcopy(a)
+        # zero!(za)
+        return za
+    end
     @eval function one(a::$T)
         b = zero(a)
         b[0] = one(b[0])
@@ -666,6 +671,39 @@ c, a and b are `HomogeneousPolynomial`.
     return nothing
 end
 
+@inline function mul_non_accumulating!(c::HomogeneousPolynomial, a::HomogeneousPolynomial,
+        b::HomogeneousPolynomial)
+
+    (iszero(b) || iszero(a)) && return nothing
+
+    @inbounds num_coeffs_a = size_table[a.order+1]
+    @inbounds num_coeffs_b = size_table[b.order+1]
+
+    @inbounds posTb = pos_table[c.order+1]
+
+    @inbounds indTa = index_table[a.order+1]
+    @inbounds indTb = index_table[b.order+1]
+
+    count = 0
+    @inbounds for na in 1:num_coeffs_a
+        ca = a[na]
+        # iszero(ca) && continue
+        inda = indTa[na]
+
+        @inbounds for nb in 1:num_coeffs_b
+            cb = b[nb]
+            # iszero(cb) && continue
+            indb = indTb[nb]
+
+            pos = posTb[inda + indb]
+            (count == 0) ? (c[pos] = ca * cb) : (c[pos] += ca * cb)
+            count += 1
+        end
+    end
+
+    return nothing
+end
+
 
 """
     mul_scalar!(c, scalar, a, b) --> nothing
diff --git a/src/auxiliary.jl b/src/auxiliary.jl
index 22c1255c..938120c7 100644
--- a/src/auxiliary.jl
+++ b/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
diff --git a/src/constructors.jl b/src/constructors.jl
index 8c8a5b3e..c3022673 100644
--- a/src/constructors.jl
+++ b/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
diff --git a/src/functions.jl b/src/functions.jl
index d4fc6001..3df81f29 100644
--- a/src/functions.jl
+++ b/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}
diff --git a/src/power.jl b/src/power.jl
index 7e9788d0..8c6a0d19 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -59,31 +59,91 @@ end
 ^(a::Taylor1{TaylorN{T}}, r::Rational) where {T<:NumberNotSeries} = a^(r.num/r.den)
 
 
-
-# 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)
+# in-place form of power_by_squaring
+function power_by_squaring!(y::TaylorN, x::TaylorN, aux1::TaylorN, aux2::TaylorN, p::Integer)
+    # this method assumes `y`, `aux1` and `aux2` are initialized to zero(x)
+    t = trailing_zeros(p) + 1
+    p >>= t
+    # aux1 = x
+    for k in eachindex(aux1)
+        identity!(aux1, x, k)
+    end
+    while (t -= 1) > 0
+        # TODO: square x in-place (in appropriate method)
+        # y = square(aux1)
+        for k in eachindex(y)
+            sqr!(y, aux1, k)
+        end
+        # aux1 = y
+        for k in eachindex(aux1)
+            identity!(aux1, y, 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
-            x = square(x)
+        while (t -= 1) ≥ 0
+            # TODO: square x in-place (in appropriate method)
+            # aux2 = square(aux1)
+            for k in eachindex(aux2)
+                sqr!(aux2, aux1, k)
+            end
+            # aux1 = aux2
+            for k in eachindex(aux1)
+                identity!(aux1, aux2, k)
+            end
+        end
+        # aux2 = y
+        for k in eachindex(aux2)
+            identity!(aux2, y, k)
         end
-        y = x
-        while p > 0
+        # y = aux2 * aux1
+        # if $T == Taylor1
+        #     for k in eachindex(y)
+        #         mul!(y, aux2, aux1, k)
+        #     end
+        # end
+        zero!(y)
+        for k in eachindex(y)
+            mul!(y, aux2, aux1, k)
+        end
+    end
+    return nothing
+end
+
+# power_by_squaring; slightly modified from base/intfuncs.jl
+# Licensed under MIT "Expat"
+for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
+    @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)
+            aux2 = zero(x)
+            power_by_squaring!(y, x, aux1, aux2, 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
@@ -300,7 +360,7 @@ end
     if ordT == lnull
         if isinteger(r)
             # TODO: get rid of allocations here
-            res[ordT] = a[l0]^round(Int,r) # uses power_by_squaring
+            @time res[ordT] = a[l0]^round(Int,r) # uses power_by_squaring
             return nothing
         end
 
@@ -339,7 +399,10 @@ 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 = $T( zero(constant_term(a)), a.order)
+        # c = deepcopy(a)
+        # zero!(c)
+        c = zero(a)
         for k in eachindex(a)
             sqr!(c, a, k)
         end
@@ -447,6 +510,45 @@ 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
+                ck = deepcopy(c[k])
+                (kodd == 0) && ( @inbounds c[k] = (c[k >> 1]^2)/2 )
+                c[k] = c[0] * ck
+                @inbounds for i = 1:kend
+                    c[k] += c[i] * c[k-i]
+                end
+                @inbounds c[k] = 2 * c[k]
+            else
+                # zz = deepcopy(c[k])
+                # zero!(c[k])
+                # @show c[0][1], typeof(c[0][1])
+                # _2c01 = 2c[0][1]
+                # @inbounds mul!(c, _2c01, c, k)
+                # mul!(c[k], zz, c[0])
+                (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])
+                    # mul!(c, 0.5, c, k)
+                end
+            end
+
+            return nothing
+        end
     end
 end
 
diff --git a/test/runtests.jl b/test/runtests.jl
index ef447b14..d041e692 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,16 +6,16 @@ testfiles = (
     "onevariable.jl",
     "manyvariables.jl",
     "mixtures.jl",
-    "mutatingfuncts.jl",
-    "intervals.jl",
-    "broadcasting.jl",
-    "identities_Euler.jl",
-    "fateman40.jl",
-    "staticarrays.jl",
-    "jld2.jl",
-    "rat.jl",
+    # "mutatingfuncts.jl",
+    # "intervals.jl",
+    # "broadcasting.jl",
+    # "identities_Euler.jl",
+    # "fateman40.jl",
+    # "staticarrays.jl",
+    # "jld2.jl",
+    # "rat.jl",
     # Run Aqua tests at the very end
-    "aqua.jl",
+    # "aqua.jl",
     )
 
 for file in testfiles

From a34ce7bde5a109935695dc5095ca2489cb52f6b5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Tue, 9 Jul 2024 15:53:29 +0200
Subject: [PATCH 02/26] Update Project.toml

---
 Project.toml | 2 --
 1 file changed, 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index eed44623..e25d96af 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,7 +8,6 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
 [weakdeps]
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
@@ -33,7 +32,6 @@ Requires = "0.5.2, 1"
 SparseArrays = "<0.0.1, 1"
 StaticArrays = "1"
 Test = "<0.0.1, 1"
-InteractiveUtils = "<0.0.1, 1"
 julia = "1"
 
 [extras]

From c04d9f80018a3cae5980d1a7b0265899d5094d81 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Tue, 9 Jul 2024 15:53:50 +0200
Subject: [PATCH 03/26] Update src/TaylorSeries.jl

---
 src/TaylorSeries.jl | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl
index 82bb84cf..2e508e66 100644
--- a/src/TaylorSeries.jl
+++ b/src/TaylorSeries.jl
@@ -28,8 +28,6 @@ using LinearAlgebra: norm, mul!,
     lu, lu!, LinearAlgebra.lutype, LinearAlgebra.copy_oftype,
     LinearAlgebra.issuccess
 
-using InteractiveUtils
-
 if VERSION >= v"1.7.0-DEV.1188"
     using LinearAlgebra: NoPivot, RowMaximum
 end

From 5d3d57db36300babeefe0069fe946d8aae63ffdb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Tue, 9 Jul 2024 15:54:52 +0200
Subject: [PATCH 04/26] Update `zero(x)`; add in-place mul! for TaylorN

---
 src/arithmetic.jl | 48 +++++++++++++----------------------------------
 1 file changed, 13 insertions(+), 35 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index abfc439e..a93a989c 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -164,8 +164,6 @@ end
 for T in (:Taylor1, :TaylorN)
     @eval function zero(a::$T)
         return $T(zero.(a.coeffs))
-        # za = deepcopy(a)
-        # zero!(za)
         return za
     end
     @eval function one(a::$T)
@@ -577,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
@@ -671,39 +682,6 @@ c, a and b are `HomogeneousPolynomial`.
     return nothing
 end
 
-@inline function mul_non_accumulating!(c::HomogeneousPolynomial, a::HomogeneousPolynomial,
-        b::HomogeneousPolynomial)
-
-    (iszero(b) || iszero(a)) && return nothing
-
-    @inbounds num_coeffs_a = size_table[a.order+1]
-    @inbounds num_coeffs_b = size_table[b.order+1]
-
-    @inbounds posTb = pos_table[c.order+1]
-
-    @inbounds indTa = index_table[a.order+1]
-    @inbounds indTb = index_table[b.order+1]
-
-    count = 0
-    @inbounds for na in 1:num_coeffs_a
-        ca = a[na]
-        # iszero(ca) && continue
-        inda = indTa[na]
-
-        @inbounds for nb in 1:num_coeffs_b
-            cb = b[nb]
-            # iszero(cb) && continue
-            indb = indTb[nb]
-
-            pos = posTb[inda + indb]
-            (count == 0) ? (c[pos] = ca * cb) : (c[pos] += ca * cb)
-            count += 1
-        end
-    end
-
-    return nothing
-end
-
 
 """
     mul_scalar!(c, scalar, a, b) --> nothing

From a9c770b27054e7fcfe6c1b555382ab95581b8e4a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Tue, 9 Jul 2024 15:57:50 +0200
Subject: [PATCH 05/26] Get rid of most allocations in pow! for
 Taylor1{TaylorN{T}; add power_by_squaring!; add in-place sqr! method

---
 src/power.jl     | 80 ++++++++++++++++++++----------------------------
 test/runtests.jl | 18 +++++------
 2 files changed, 42 insertions(+), 56 deletions(-)

diff --git a/src/power.jl b/src/power.jl
index 8c6a0d19..afc7218c 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -60,8 +60,8 @@ end
 
 
 # in-place form of power_by_squaring
-function power_by_squaring!(y::TaylorN, x::TaylorN, aux1::TaylorN, aux2::TaylorN, p::Integer)
-    # this method assumes `y`, `aux1` and `aux2` are initialized to zero(x)
+# 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
@@ -69,14 +69,9 @@ function power_by_squaring!(y::TaylorN, x::TaylorN, aux1::TaylorN, aux2::TaylorN
         identity!(aux1, x, k)
     end
     while (t -= 1) > 0
-        # TODO: square x in-place (in appropriate method)
-        # y = square(aux1)
-        for k in eachindex(y)
-            sqr!(y, aux1, k)
-        end
-        # aux1 = y
-        for k in eachindex(aux1)
-            identity!(aux1, y, k)
+        # aux1 = square(aux1)
+        for k in reverse(eachindex(aux1))
+            sqr!(aux1, k)
         end
     end
     # y = aux1
@@ -87,34 +82,18 @@ function power_by_squaring!(y::TaylorN, x::TaylorN, aux1::TaylorN, aux2::TaylorN
         t = trailing_zeros(p) + 1
         p >>= t
         while (t -= 1) ≥ 0
-            # TODO: square x in-place (in appropriate method)
-            # aux2 = square(aux1)
-            for k in eachindex(aux2)
-                sqr!(aux2, aux1, k)
-            end
-            # aux1 = aux2
-            for k in eachindex(aux1)
-                identity!(aux1, aux2, k)
+            # aux1 = square(aux1)
+            for k in reverse(eachindex(aux1))
+                sqr!(aux1, k)
             end
         end
-        # aux2 = y
-        for k in eachindex(aux2)
-            identity!(aux2, y, k)
-        end
-        # y = aux2 * aux1
-        # if $T == Taylor1
-        #     for k in eachindex(y)
-        #         mul!(y, aux2, aux1, k)
-        #     end
-        # end
-        zero!(y)
-        for k in eachindex(y)
-            mul!(y, aux2, aux1, k)
-        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)
@@ -127,8 +106,7 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
         if $T == TaylorN
             y = zero(x)
             aux1 = zero(x)
-            aux2 = zero(x)
-            power_by_squaring!(y, x, aux1, aux2, p)
+            power_by_squaring!(y, x, aux1, p)
         else
             t = trailing_zeros(p) + 1
             p >>= t
@@ -149,6 +127,24 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
     end
 end
 
+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}
@@ -360,7 +356,8 @@ end
     if ordT == lnull
         if isinteger(r)
             # TODO: get rid of allocations here
-            @time 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
 
@@ -399,9 +396,6 @@ 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 = deepcopy(a)
-        # zero!(c)
         c = zero(a)
         for k in eachindex(a)
             sqr!(c, a, k)
@@ -522,20 +516,13 @@ for T = (:Taylor1, :TaylorN)
             kodd = k%2
             kend = (k - 2 + kodd) >> 1
             if $T == Taylor1
-                ck = deepcopy(c[k])
                 (kodd == 0) && ( @inbounds c[k] = (c[k >> 1]^2)/2 )
-                c[k] = c[0] * ck
+                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
-                # zz = deepcopy(c[k])
-                # zero!(c[k])
-                # @show c[0][1], typeof(c[0][1])
-                # _2c01 = 2c[0][1]
-                # @inbounds mul!(c, _2c01, c, k)
-                # mul!(c[k], zz, c[0])
                 (kend ≥ 0) && mul!(c, c[0][1], c, k)
                 @inbounds for i = 1:kend
                     mul!(c[k], c[i], c[k-i])
@@ -543,7 +530,6 @@ for T = (:Taylor1, :TaylorN)
                 @inbounds mul!(c, 2, c, k)
                 if (kodd == 0)
                     accsqr!(c[k], c[k >> 1])
-                    # mul!(c, 0.5, c, k)
                 end
             end
 
diff --git a/test/runtests.jl b/test/runtests.jl
index d041e692..ef447b14 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,16 +6,16 @@ testfiles = (
     "onevariable.jl",
     "manyvariables.jl",
     "mixtures.jl",
-    # "mutatingfuncts.jl",
-    # "intervals.jl",
-    # "broadcasting.jl",
-    # "identities_Euler.jl",
-    # "fateman40.jl",
-    # "staticarrays.jl",
-    # "jld2.jl",
-    # "rat.jl",
+    "mutatingfuncts.jl",
+    "intervals.jl",
+    "broadcasting.jl",
+    "identities_Euler.jl",
+    "fateman40.jl",
+    "staticarrays.jl",
+    "jld2.jl",
+    "rat.jl",
     # Run Aqua tests at the very end
-    # "aqua.jl",
+    "aqua.jl",
     )
 
 for file in testfiles

From 5bcccf42d6acba54f0ab69084e9767621201b544 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Tue, 9 Jul 2024 16:00:12 +0200
Subject: [PATCH 06/26] Remove extra empty line

---
 src/power.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/power.jl b/src/power.jl
index afc7218c..6c9b2daa 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -58,7 +58,6 @@ 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}

From 6cce90ade660a4a583236a26cc64d93453160e4f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Tue, 9 Jul 2024 16:16:57 +0200
Subject: [PATCH 07/26] Add power_by_squaring methods to avoid method
 ambiguities (detected by Aqua)

---
 src/power.jl | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/power.jl b/src/power.jl
index 6c9b2daa..c8a2597d 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -126,6 +126,9 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
     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

From a205dcd5b763c7578d5e73d8c6a2935d5fc78caa Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Wed, 10 Jul 2024 22:21:10 +0200
Subject: [PATCH 08/26] Remove duplicated return

---
 src/arithmetic.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index a93a989c..e65e1a9b 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -164,7 +164,6 @@ end
 for T in (:Taylor1, :TaylorN)
     @eval function zero(a::$T)
         return $T(zero.(a.coeffs))
-        return za
     end
     @eval function one(a::$T)
         b = zero(a)

From e44811cfe8220fd5e6935922ee5131d32a5b1fe6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Wed, 10 Jul 2024 22:42:45 +0200
Subject: [PATCH 09/26] Add inbounds

---
 src/arithmetic.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index e65e1a9b..28c8fa45 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -578,7 +578,7 @@ end
 # 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))
+    @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])

From 2e4e60af40611b37491440f65e97f43d2cd3cbe5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Thu, 11 Jul 2024 11:03:27 +0200
Subject: [PATCH 10/26] Address review by @lbenet

---
 src/arithmetic.jl | 23 ++++++++++++++++++++
 src/power.jl      | 55 +++++++++++++++++++++++++----------------------
 2 files changed, 52 insertions(+), 26 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 28c8fa45..4c8c0a9a 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -586,6 +586,29 @@ function mul!(a::TaylorN{T}, b::TaylorN{T}) where {T<:Number}
     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}
diff --git a/src/power.jl b/src/power.jl
index c8a2597d..5ebaba27 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -60,36 +60,39 @@ end
 
 # 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
+for T in (:Taylor1, :TaylorN)
+    @eval function power_by_squaring!(y::$T{T}, x::$T{T}, aux1::$T{T},
+            p::Integer) where {T<:NumberNotSeries}
         t = trailing_zeros(p) + 1
         p >>= t
-        while (t -= 1) ≥ 0
+        # 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 = y * aux1
-        mul!(y, aux1)
+        # 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
-    return nothing
 end
 
 
@@ -102,7 +105,7 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
     @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
+        if $T != HomogeneousPolynomial
             y = zero(x)
             aux1 = zero(x)
             power_by_squaring!(y, x, aux1, p)
@@ -518,14 +521,14 @@ for T = (:Taylor1, :TaylorN)
             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]
+                (kend ≥ 0) && ( @inbounds 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]
+                (kodd == 0) && ( @inbounds c[k] = c[k >> 1]^2 )
             else
-                (kend ≥ 0) && mul!(c, c[0][1], c, k)
+                (kend ≥ 0) && ( @inbounds mul!(c, c[0][1], c, k) )
                 @inbounds for i = 1:kend
                     mul!(c[k], c[i], c[k-i])
                 end

From f581b215f3db1c1ac97ab7f4e3d356c8303d99aa Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Thu, 11 Jul 2024 12:28:19 +0200
Subject: [PATCH 11/26] Add extra arg to pow!

---
 src/dictmutfunct.jl    |  3 +--
 src/evaluate.jl        |  2 +-
 src/power.jl           | 54 ++++++++++++++++++++++--------------------
 test/manyvariables.jl  | 14 +++++------
 test/mutatingfuncts.jl |  3 ++-
 test/onevariable.jl    | 22 +++++++++--------
 6 files changed, 51 insertions(+), 47 deletions(-)

diff --git a/src/dictmutfunct.jl b/src/dictmutfunct.jl
index d5ec11cf..87f3b63f 100644
--- a/src/dictmutfunct.jl
+++ b/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, :_arg2, :_arg3, :_k), :(_res = _arg1 ^ float(_arg3))),
 );
 
 """
@@ -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
-
diff --git a/src/evaluate.jl b/src/evaluate.jl
index a3ecc544..7cc23a67 100644
--- a/src/evaluate.jl
+++ b/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)
diff --git a/src/power.jl b/src/power.jl
index 5ebaba27..3f1403e0 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -59,37 +59,37 @@ 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
+# this method assumes `y`, `x` and `aux` are of same order
 for T in (:Taylor1, :TaylorN)
-    @eval function power_by_squaring!(y::$T{T}, x::$T{T}, aux1::$T{T},
+    @eval function power_by_squaring!(y::$T{T}, x::$T{T}, aux::$T{T},
             p::Integer) where {T<:NumberNotSeries}
         t = trailing_zeros(p) + 1
         p >>= t
-        # aux1 = x
-        for k in eachindex(aux1)
-            identity!(aux1, x, k)
+        # aux = x
+        for k in eachindex(aux)
+            identity!(aux, x, k)
         end
         while (t -= 1) > 0
-            # aux1 = square(aux1)
-            for k in reverse(eachindex(aux1))
-                sqr!(aux1, k)
+            # aux = square(aux)
+            for k in reverse(eachindex(aux))
+                sqr!(aux, k)
             end
         end
-        # y = aux1
+        # y = aux
         for k in eachindex(y)
-            identity!(y, aux1, k)
+            identity!(y, aux, 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)
+                # aux = square(aux)
+                for k in reverse(eachindex(aux))
+                    sqr!(aux, k)
                 end
             end
-            # y = y * aux1
-            mul!(y, aux1)
+            # y = y * aux
+            mul!(y, aux)
         end
         return nothing
     end
@@ -107,8 +107,8 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
         p = P # copy static parameter `P` into local variable `p`
         if $T != HomogeneousPolynomial
             y = zero(x)
-            aux1 = zero(x)
-            power_by_squaring!(y, x, aux1, p)
+            aux = zero(x)
+            power_by_squaring!(y, x, aux, p)
         else
             t = trailing_zeros(p) + 1
             p >>= t
@@ -168,8 +168,9 @@ function ^(a::Taylor1{T}, r::S) where {T<:Number, S<:Real}
 
     c_order = l0 == 0 ? a.order : min(a.order, trunc(Int,r*a.order))
     c = Taylor1(zero(aux), c_order)
+    aux0 = deepcopy(c[0])
     for k in eachindex(c)
-        pow!(c, aa, r, k)
+        pow!(c, aa, aux0, r, k)
     end
 
     return c
@@ -193,8 +194,9 @@ function ^(a::TaylorN, r::S) where {S<:Real}
         in order to expand `^` around 0."""))
 
     c = TaylorN( zero(aux), a.order)
+    aux = deepcopy(c)
     for ord in eachindex(a)
-        pow!(c, aa, r, ord)
+        pow!(c, aa, aux, r, ord)
     end
 
     return c
@@ -219,8 +221,9 @@ function ^(a::Taylor1{TaylorN{T}}, r::S) where {T<:NumberNotSeries, S<:Real}
 
     c_order = l0 == 0 ? a.order : min(a.order, trunc(Int,r*a.order))
     c = Taylor1(zero(aux), c_order)
+    aux0 = deepcopy(c[0])
     for k in eachindex(c)
-        pow!(c, aa, r, k)
+        pow!(c, aa, aux0, r, k)
     end
 
     return c
@@ -245,7 +248,7 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`.
 
 """ pow!
 
-@inline function pow!(c::Taylor1{T}, a::Taylor1{T}, r::S, k::Int) where
+@inline function pow!(c::Taylor1{T}, a::Taylor1{T}, ::T, r::S, k::Int) where
         {T<:Number, S<:Real}
 
     if r == 0
@@ -294,7 +297,7 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`.
     return nothing
 end
 
-@inline function pow!(c::TaylorN{T}, a::TaylorN{T}, r::S, k::Int) where
+@inline function pow!(c::TaylorN{T}, a::TaylorN{T}, ::TaylorN{T}, r::S, k::Int) where
         {T<:NumberNotSeriesN, S<:Real}
 
     if r == 0
@@ -327,7 +330,7 @@ end
     return nothing
 end
 
-@inline function pow!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, r::S,
+@inline function pow!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, aux::TaylorN{T}, r::S,
         ordT::Int) where {T<:NumberNotSeries, S<:Real}
 
     if r == 0
@@ -361,8 +364,7 @@ end
     if ordT == lnull
         if isinteger(r)
             # TODO: get rid of allocations here
-            aux1 = deepcopy(res[ordT])
-            power_by_squaring!(res[ordT], a[l0], aux1, round(Int,r))
+            power_by_squaring!(res[ordT], a[l0], aux, round(Int,r))
             return nothing
         end
 
@@ -372,7 +374,7 @@ end
             in order to expand `^` around 0."""))
 
         for ordQ in eachindex(a[l0])
-            pow!(res[ordT], a[l0], r, ordQ)
+            pow!(res[ordT], a[l0], aux, r, ordQ)
         end
 
         return nothing
diff --git a/test/manyvariables.jl b/test/manyvariables.jl
index 1b7f243b..b765d7bf 100644
--- a/test/manyvariables.jl
+++ b/test/manyvariables.jl
@@ -507,22 +507,22 @@ end
     @test xx[1] == zero(xx[end])
     TS.div!(xx, 1.0+xT, 1.0+xT, 0)
     @test xx[0] == HomogeneousPolynomial([1.0])
-    TS.pow!(xx, 1.0+xT, 0.5, 1)
+    TS.pow!(xx, 1.0+xT, xx, 0.5, 1)
     @test xx[1] == HomogeneousPolynomial([0.5,0.0])
     xx = 1.0*zeroT
-    TS.pow!(xx, 1.0+xT, 1.5, 0)
+    TS.pow!(xx, 1.0+xT, xx, 1.5, 0)
     @test xx[0] == HomogeneousPolynomial([1.0])
-    TS.pow!(xx, 1.0+xT, 1.5, 1)
+    TS.pow!(xx, 1.0+xT, xx, 1.5, 1)
     @test xx[1] == HomogeneousPolynomial([1.5,0.0])
     xx = 1.0*zeroT
-    TS.pow!(xx, 1.0+xT, 0, 0)
+    TS.pow!(xx, 1.0+xT, xx, 0, 0)
     @test xx[0] == HomogeneousPolynomial([1.0])
-    TS.pow!(xx, 1.0+xT, 1, 1)
+    TS.pow!(xx, 1.0+xT, xx, 1, 1)
     @test xx[1] == HomogeneousPolynomial([1.0,0.0])
     xx = 1.0*zeroT
-    TS.pow!(xx, 1.0+xT, 2, 0)
+    TS.pow!(xx, 1.0+xT, xx, 2, 0)
     @test xx[0] == HomogeneousPolynomial([1.0])
-    TS.pow!(xx, 1.0+xT, 2, 1)
+    TS.pow!(xx, 1.0+xT, xx, 2, 1)
     @test xx[1] == HomogeneousPolynomial([2.0,0.0])
     xx = 1.0*zeroT
     TS.sqrt!(xx, 1.0+xT, 0)
diff --git a/test/mutatingfuncts.jl b/test/mutatingfuncts.jl
index 6599391f..db74a5d6 100644
--- a/test/mutatingfuncts.jl
+++ b/test/mutatingfuncts.jl
@@ -30,7 +30,8 @@ using Test
     # Some examples
     t1 = Taylor1(5)
     t2 = zero(t1)
-    TS.pow!(t2, t1, 2, 2)
+    t10 = deepcopy(t1[0])
+    TS.pow!(t2, t1, t10, 2, 2)
     @test t2[2] == 1.0
     #
     res = zero(t1)
diff --git a/test/onevariable.jl b/test/onevariable.jl
index 120a5a9f..8de3a602 100644
--- a/test/onevariable.jl
+++ b/test/onevariable.jl
@@ -407,26 +407,28 @@ Base.iszero(::SymbNumber) = false
     @test tt[0] == 1.0
     TS.div!(tt, 1, 1+ut, 0)
     @test tt[0] == 1.0
-    TS.pow!(tt, 1.0+t, 1.5, 0)
+    aux = deepcopy(tt[0])
+    TS.pow!(tt, 1.0+t, aux, 1.5, 0)
     @test tt[0] == 1.0
-    TS.pow!(tt, 0.0*t, 1.5, 0)
+    TS.pow!(tt, 0.0*t, aux, 1.5, 0)
     @test tt[0] == 0.0
-    TS.pow!(tt, 0.0+t, 18, 0)
+    TS.pow!(tt, 0.0+t, aux, 18, 0)
     @test tt[0] == 0.0
-    TS.pow!(tt, 1.0+t, 1.5, 0)
+    TS.pow!(tt, 1.0+t, aux, 1.5, 0)
     @test tt[0] == 1.0
-    TS.pow!(tt, 1.0+t, 0.5, 1)
+    TS.pow!(tt, 1.0+t, aux, 0.5, 1)
     @test tt[1] == 0.5
-    TS.pow!(tt, 1.0+t, 0, 0)
+    TS.pow!(tt, 1.0+t, aux, 0, 0)
     @test tt[0] == 1.0
-    TS.pow!(tt, 1.0+t, 1, 1)
+    TS.pow!(tt, 1.0+t, aux, 1, 1)
     @test tt[1] == 1.0
     tt = zero(ut)
-    TS.pow!(tt, 1.0+t, 2, 0)
+    aux = deepcopy(tt[0])
+    TS.pow!(tt, 1.0+t, aux, 2, 0)
     @test tt[0] == 1.0
-    TS.pow!(tt, 1.0+t, 2, 1)
+    TS.pow!(tt, 1.0+t, aux, 2, 1)
     @test tt[1] == 2.0
-    TS.pow!(tt, 1.0+t, 2, 2)
+    TS.pow!(tt, 1.0+t, aux, 2, 2)
     @test tt[2] == 1.0
     TS.sqrt!(tt, 1.0+t, 0, 0)
     @test tt[0] == 1.0

From 27c3d0e65d8c4bf45c9561eae4993791269a30b7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Thu, 11 Jul 2024 14:10:33 +0200
Subject: [PATCH 12/26] Update TaylorSeries IA extension

---
 ext/TaylorSeriesIAExt.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/ext/TaylorSeriesIAExt.jl b/ext/TaylorSeriesIAExt.jl
index 30ecc7c4..72c1d987 100644
--- a/ext/TaylorSeriesIAExt.jl
+++ b/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[0], 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

From dc6a68cfa6e3ed86199436bbca05fdcd9636bdc0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Thu, 11 Jul 2024 22:03:28 +0200
Subject: [PATCH 13/26] Update pow! and add fix suggested by @lbenet

---
 ext/TaylorSeriesIAExt.jl |  2 +-
 src/dictmutfunct.jl      |  2 +-
 src/power.jl             | 14 +++++++-------
 test/mutatingfuncts.jl   |  4 ++--
 test/onevariable.jl      |  4 ++--
 5 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/ext/TaylorSeriesIAExt.jl b/ext/TaylorSeriesIAExt.jl
index 72c1d987..d526bd16 100644
--- a/ext/TaylorSeriesIAExt.jl
+++ b/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, c[0], r, k)
+        TS.pow!(c, aa, c, r, k)
     end
 
     return c
diff --git a/src/dictmutfunct.jl b/src/dictmutfunct.jl
index 87f3b63f..7d1870ab 100644
--- a/src/dictmutfunct.jl
+++ b/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, :_arg3, :_k), :(_res = _arg1 ^ float(_arg3))),
+    :^ => (:pow!, (:_res, :_arg1, :_aux, :_arg2, :_k), :(_res = _arg1 ^ float(_arg2)), :(_aux = zero(_arg1))),
 );
 
 """
diff --git a/src/power.jl b/src/power.jl
index 3f1403e0..afd131d8 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -168,7 +168,7 @@ function ^(a::Taylor1{T}, r::S) where {T<:Number, S<:Real}
 
     c_order = l0 == 0 ? a.order : min(a.order, trunc(Int,r*a.order))
     c = Taylor1(zero(aux), c_order)
-    aux0 = deepcopy(c[0])
+    aux0 = deepcopy(c)
     for k in eachindex(c)
         pow!(c, aa, aux0, r, k)
     end
@@ -221,7 +221,7 @@ function ^(a::Taylor1{TaylorN{T}}, r::S) where {T<:NumberNotSeries, S<:Real}
 
     c_order = l0 == 0 ? a.order : min(a.order, trunc(Int,r*a.order))
     c = Taylor1(zero(aux), c_order)
-    aux0 = deepcopy(c[0])
+    aux0 = deepcopy(c)
     for k in eachindex(c)
         pow!(c, aa, aux0, r, k)
     end
@@ -248,7 +248,7 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`.
 
 """ pow!
 
-@inline function pow!(c::Taylor1{T}, a::Taylor1{T}, ::T, r::S, k::Int) where
+@inline function pow!(c::Taylor1{T}, a::Taylor1{T}, ::Taylor1{T}, r::S, k::Int) where
         {T<:Number, S<:Real}
 
     if r == 0
@@ -330,7 +330,7 @@ end
     return nothing
 end
 
-@inline function pow!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, aux::TaylorN{T}, r::S,
+@inline function pow!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, aux::Taylor1{TaylorN{T}}, r::S,
         ordT::Int) where {T<:NumberNotSeries, S<:Real}
 
     if r == 0
@@ -364,7 +364,7 @@ end
     if ordT == lnull
         if isinteger(r)
             # TODO: get rid of allocations here
-            power_by_squaring!(res[ordT], a[l0], aux, round(Int,r))
+            power_by_squaring!(res[ordT], a[l0], aux[0], round(Int,r))
             return nothing
         end
 
@@ -374,7 +374,7 @@ end
             in order to expand `^` around 0."""))
 
         for ordQ in eachindex(a[l0])
-            pow!(res[ordT], a[l0], aux, r, ordQ)
+            pow!(res[ordT], a[l0], aux[0], r, ordQ)
         end
 
         return nothing
@@ -528,7 +528,7 @@ for T = (:Taylor1, :TaylorN)
                     c[k] += c[i] * c[k-i]
                 end
                 @inbounds c[k] = 2 * c[k]
-                (kodd == 0) && ( @inbounds c[k] = c[k >> 1]^2 )
+                (kodd == 0) && ( @inbounds c[k] += c[k >> 1]^2 )
             else
                 (kend ≥ 0) && ( @inbounds mul!(c, c[0][1], c, k) )
                 @inbounds for i = 1:kend
diff --git a/test/mutatingfuncts.jl b/test/mutatingfuncts.jl
index db74a5d6..fb9be7ef 100644
--- a/test/mutatingfuncts.jl
+++ b/test/mutatingfuncts.jl
@@ -30,8 +30,8 @@ using Test
     # Some examples
     t1 = Taylor1(5)
     t2 = zero(t1)
-    t10 = deepcopy(t1[0])
-    TS.pow!(t2, t1, t10, 2, 2)
+    t1aux = deepcopy(t1)
+    TS.pow!(t2, t1, t1aux, 2, 2)
     @test t2[2] == 1.0
     #
     res = zero(t1)
diff --git a/test/onevariable.jl b/test/onevariable.jl
index 8de3a602..40bb13a0 100644
--- a/test/onevariable.jl
+++ b/test/onevariable.jl
@@ -407,7 +407,7 @@ Base.iszero(::SymbNumber) = false
     @test tt[0] == 1.0
     TS.div!(tt, 1, 1+ut, 0)
     @test tt[0] == 1.0
-    aux = deepcopy(tt[0])
+    aux = tt
     TS.pow!(tt, 1.0+t, aux, 1.5, 0)
     @test tt[0] == 1.0
     TS.pow!(tt, 0.0*t, aux, 1.5, 0)
@@ -423,7 +423,7 @@ Base.iszero(::SymbNumber) = false
     TS.pow!(tt, 1.0+t, aux, 1, 1)
     @test tt[1] == 1.0
     tt = zero(ut)
-    aux = deepcopy(tt[0])
+    aux = tt
     TS.pow!(tt, 1.0+t, aux, 2, 0)
     @test tt[0] == 1.0
     TS.pow!(tt, 1.0+t, aux, 2, 1)

From 364faf7f51478c2e67c821345f76ab75301a5c53 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Fri, 12 Jul 2024 13:59:04 +0200
Subject: [PATCH 14/26] Middle-of-the-road approach to suggestion by @lbenet

---
 src/power.jl | 54 +++++++++++++++++++++++++++++-----------------------
 1 file changed, 30 insertions(+), 24 deletions(-)

diff --git a/src/power.jl b/src/power.jl
index afd131d8..68aba8c1 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -99,35 +99,41 @@ end
 # power_by_squaring; slightly modified from base/intfuncs.jl
 # Licensed under MIT "Expat"
 for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
-    @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`
+    @eval begin
+        power_by_squaring(x::$T, p::Integer) = power_by_squaring(x, Val(p))
+        power_by_squaring(x::$T, ::Val{0}) = one(x)
+        power_by_squaring(x::$T, ::Val{1}) = copy(x)
+        power_by_squaring(x::$T, ::Val{2}) = square(x)
         if $T != HomogeneousPolynomial
-            y = zero(x)
-            aux = zero(x)
-            power_by_squaring!(y, x, aux, p)
-        else
-            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
+            function power_by_squaring(x::$T, ::Val{P}) where P
+                p = P # copy static parameter `P` into local variable `p`
+                y = zero(x)
+                aux = zero(x)
+                power_by_squaring!(y, x, aux, p)
+                return y
             end
         end
-        return y
     end
 end
+function power_by_squaring(x::HomogeneousPolynomial, ::Val{P}) where P
+    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
+
 
 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)

From 99e3061fd03cca29e8cabff1716f36f4ce7f8940 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Fri, 12 Jul 2024 16:01:14 +0200
Subject: [PATCH 15/26] Handle pow! cases by dispatch-by-value

---
 src/power.jl | 90 +++++++++++++++++++++++++++++++---------------------
 1 file changed, 54 insertions(+), 36 deletions(-)

diff --git a/src/power.jl b/src/power.jl
index 68aba8c1..91d947ed 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -254,18 +254,24 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`.
 
 """ pow!
 
-@inline function pow!(c::Taylor1{T}, a::Taylor1{T}, ::Taylor1{T}, r::S, k::Int) where
-        {T<:Number, S<:Real}
-
-    if r == 0
-        return one!(c, a, k)
-    elseif r == 1
-        return identity!(c, a, k)
-    elseif r == 2
-        return sqr!(c, a, k)
-    elseif r == 0.5
-        return sqrt!(c, a, k)
-    end
+# pow! main dispatcher
+@inline pow!(c::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}, r::S, k::Int) where
+        {T<:Number, S <: Real} = pow!(c, a, b, Val(r), k)
+# pow! dispatches by value
+@inline pow!(c::Taylor1{T}, a::Taylor1{T}, ::Taylor1{T}, ::Val{0}, k::Int) where
+        {T<:Number} = one!(c, a, k)
+@inline pow!(c::Taylor1{T}, a::Taylor1{T}, ::Taylor1{T}, ::Val{1}, k::Int) where
+        {T<:Number} = identity!(c, a, k)
+@inline pow!(c::Taylor1{T}, a::Taylor1{T}, ::Taylor1{T}, ::Val{2}, k::Int) where
+        {T<:Number} = sqr!(c, a, k)
+@inline pow!(c::Taylor1{T}, a::Taylor1{T}, ::Taylor1{T}, ::Val{0.5}, k::Int) where
+        {T<:Number} = sqrt!(c, a, k)
+# fallback pow! method
+@inline function pow!(c::Taylor1{T}, a::Taylor1{T}, ::Taylor1{T}, ::Val{R}, k::Int) where
+        {T<:Number, R}
+
+    # copy static parameter `R` into local variable `r`
+    r = R
 
     # Sanity
     zero!(c, k)
@@ -303,18 +309,24 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`.
     return nothing
 end
 
-@inline function pow!(c::TaylorN{T}, a::TaylorN{T}, ::TaylorN{T}, r::S, k::Int) where
-        {T<:NumberNotSeriesN, S<:Real}
-
-    if r == 0
-        return one!(c, a, k)
-    elseif r == 1
-        return identity!(c, a, k)
-    elseif r == 2
-        return sqr!(c, a, k)
-    elseif r == 0.5
-        return sqrt!(c, a, k)
-    end
+# pow! main dispatcher
+@inline pow!(c::TaylorN{T}, a::TaylorN{T}, b::TaylorN{T}, r::S, k::Int) where
+        {T<:NumberNotSeriesN, S <: Real} = pow!(c, a, b, Val(r), k)
+# pow! dispatches by value
+@inline pow!(c::TaylorN{T}, a::TaylorN{T}, ::TaylorN{T}, ::Val{0}, k::Int) where
+        {T<:NumberNotSeriesN} = one!(c, a, k)
+@inline pow!(c::TaylorN{T}, a::TaylorN{T}, ::TaylorN{T}, ::Val{1}, k::Int) where
+        {T<:NumberNotSeriesN} = identity!(c, a, k)
+@inline pow!(c::TaylorN{T}, a::TaylorN{T}, ::TaylorN{T}, ::Val{2}, k::Int) where
+        {T<:NumberNotSeriesN} = sqr!(c, a, k)
+@inline pow!(c::TaylorN{T}, a::TaylorN{T}, ::TaylorN{T}, ::Val{0.5}, k::Int) where
+        {T<:NumberNotSeriesN} = sqrt!(c, a, k)
+# fallback pow! method
+@inline function pow!(c::TaylorN{T}, a::TaylorN{T}, ::TaylorN{T}, ::Val{R}, k::Int) where
+        {T<:NumberNotSeriesN, R}
+
+    # copy static parameter `R` into local variable `r`
+    r = R
 
     if k == 0
         @inbounds c[0][1] = ( constant_term(a) )^r
@@ -336,18 +348,24 @@ end
     return nothing
 end
 
-@inline function pow!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, aux::Taylor1{TaylorN{T}}, r::S,
-        ordT::Int) where {T<:NumberNotSeries, S<:Real}
-
-    if r == 0
-        return one!(res, ordT)
-    elseif r == 1
-        return identity!(res, a, ordT)
-    elseif r == 2
-        return sqr!(res, a, ordT)
-    elseif r == 0.5
-        return sqrt!(res, a, ordT)
-    end
+# pow! main dispatcher for `Taylor1{TaylorN{T}}`
+@inline pow!(c::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}}, r::S, k::Int) where
+        {T<:NumberNotSeries, S <: Real} = pow!(c, a, b, Val(r), k)
+# pow! dispatches by value for `Taylor1{TaylorN{T}}`
+@inline pow!(c::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, ::Taylor1{TaylorN{T}}, ::Val{0}, k::Int) where
+        {T<:NumberNotSeries} = one!(c, a, k)
+@inline pow!(c::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, ::Taylor1{TaylorN{T}}, ::Val{1}, k::Int) where
+        {T<:NumberNotSeries} = identity!(c, a, k)
+@inline pow!(c::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, ::Taylor1{TaylorN{T}}, ::Val{2}, k::Int) where
+        {T<:NumberNotSeries} = sqr!(c, a, k)
+@inline pow!(c::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, ::Taylor1{TaylorN{T}}, ::Val{0.5}, k::Int) where
+        {T<:NumberNotSeries} = sqrt!(c, a, k)
+# fallback pow! method for `Taylor1{TaylorN{T}}`
+@inline function pow!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, aux::Taylor1{TaylorN{T}}, ::Val{R},
+        ordT::Int) where {T<:NumberNotSeries, R}
+
+    # copy static parameter `R` into local variable `r`
+    r = R
 
     # Sanity
     zero!(res, ordT)

From efa02a18c1051e127853a3c1d188bb18a27bc790 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Fri, 12 Jul 2024 16:29:34 +0200
Subject: [PATCH 16/26] Remove unneeded deepcopy in setindex! method

---
 src/auxiliary.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/auxiliary.jl b/src/auxiliary.jl
index 938120c7..4a5e855a 100644
--- a/src/auxiliary.jl
+++ b/src/auxiliary.jl
@@ -82,7 +82,7 @@ getindex(a::Taylor1{T}, u::StepRange{Int,Int}) where {T<:Number} =
     view(a.coeffs, u[:] .+ 1)
 
 setindex!(a::Taylor1{T}, x::T, n::Int) where {T<:Number} = a.coeffs[n+1] = x
-setindex!(a::Taylor1{T}, x::T, n::Int) where {T<:AbstractSeries} = setindex!(a.coeffs, deepcopy(x), n+1)
+setindex!(a::Taylor1{T}, x::T, n::Int) where {T<:AbstractSeries} = setindex!(a.coeffs, x, n+1)
 setindex!(a::Taylor1{T}, x::T, u::UnitRange{Int}) where {T<:Number} =
     a.coeffs[u .+ 1] .= x
 function setindex!(a::Taylor1{T}, x::Array{T,1}, u::UnitRange{Int}) where {T<:Number}

From d9af8ca9e1bd75d50916729337efedf1c26aadbc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Sun, 14 Jul 2024 12:42:11 +0200
Subject: [PATCH 17/26] Add `power_by_squaring(x, ::Val{3})` methods and add
 tests

---
 src/power.jl          |  2 ++
 test/manyvariables.jl | 20 ++++++++++++++++++++
 test/onevariable.jl   |  8 ++++++++
 3 files changed, 30 insertions(+)

diff --git a/src/power.jl b/src/power.jl
index 91d947ed..f95f66d6 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -104,6 +104,7 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
         power_by_squaring(x::$T, ::Val{0}) = one(x)
         power_by_squaring(x::$T, ::Val{1}) = copy(x)
         power_by_squaring(x::$T, ::Val{2}) = square(x)
+        power_by_squaring(x::$T, ::Val{3}) = x*square(x)
         if $T != HomogeneousPolynomial
             function power_by_squaring(x::$T, ::Val{P}) where P
                 p = P # copy static parameter `P` into local variable `p`
@@ -138,6 +139,7 @@ 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)
+power_by_squaring(x::TaylorN{Taylor1{T}}, ::Val{3}) where {T<:NumberNotSeries} = x*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
diff --git a/test/manyvariables.jl b/test/manyvariables.jl
index b765d7bf..32b4e88b 100644
--- a/test/manyvariables.jl
+++ b/test/manyvariables.jl
@@ -117,6 +117,26 @@ end
     @test findfirst(hpol_v) == -1
     @test findlast(hpol_v) == -1
 
+    xv = TaylorSeries.get_variables()
+    @test (xv[1] - 1.0)^3 == -1 + 3xv[1] - 3xv[1]^2 + xv[1]^3
+    @test (xv[1] - 1.0)^4 == 1 - 4xv[1] + 6xv[1]^2 - 4xv[1]^3 + xv[1]^4
+    @test (xv[2] - 1.0)^3 == -1 + 3xv[2] - 3xv[2]^2 + xv[2]^3
+    @test (xv[2] - 1.0)^4 == 1 - 4xv[2] + 6xv[2]^2 - 4xv[2]^3 + xv[2]^4
+    xN = 5.0 + 3xv[1] - 4.5xv[2] + 6.125xv[1]^2 - 7.25xv[1]*xv[2] - 5.5xv[2]^2
+    xNcopy = deepcopy(xN)
+    for k in reverse(eachindex(xNcopy))
+        TaylorSeries.sqr!(xNcopy, k)
+    end
+    @test xNcopy == xN^2
+    @test xN^2 == Base.power_by_squaring(xN, 2)
+    @test xN*xN*xN == xN*xNcopy
+    @test xN^3 == Base.power_by_squaring(xN, 3)
+    for k in reverse(eachindex(xNcopy))
+        TaylorSeries.sqr!(xNcopy, k)
+    end
+    @test xNcopy == xN^4
+    @test xN^4 == Base.power_by_squaring(xN, 4)
+
     tn_v = TaylorN(HomogeneousPolynomial(zeros(Int, 3)))
     tn_v[0] = 1
     @test tn_v == 1
diff --git a/test/onevariable.jl b/test/onevariable.jl
index 40bb13a0..847ff8b3 100644
--- a/test/onevariable.jl
+++ b/test/onevariable.jl
@@ -618,6 +618,14 @@ Base.iszero(::SymbNumber) = false
     #     @test a[i]*(180/pi) == b[i]
     # end
 
+    x = Taylor1([5.0,-1.5,3.0,-2.0,-20.0])
+    @test x*x == x^2
+    @test x*x*x == x*(x^2) == TaylorSeries.power_by_squaring(x, 3)
+    @test x*x*x*x == (x^2)*(x^2) == TaylorSeries.power_by_squaring(x, 4)
+    @test (x - 1.0)^2 == 1 - 2x + x^2
+    @test (x - 1.0)^3 == -1 + 3x - 3x^2 + x^3
+    @test (x - 1.0)^4 == 1 - 4x + 6x^2 - 4x^3 + x^4
+
     # Test additional Taylor1 constructors
     @test Taylor1{Float64}(true) == Taylor1([1.0])
     @test Taylor1{Float64}(false) == Taylor1([0.0])

From 68bf845b6d9f8dec474f9138ec10b872ef3f60b9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Sun, 14 Jul 2024 12:46:39 +0200
Subject: [PATCH 18/26] Update comments

---
 src/power.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/power.jl b/src/power.jl
index f95f66d6..b53a3d53 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -60,6 +60,7 @@ end
 
 # in-place form of power_by_squaring
 # this method assumes `y`, `x` and `aux` are of same order
+# TODO: add power_by_squaring! method for HomogeneousPolynomial
 for T in (:Taylor1, :TaylorN)
     @eval function power_by_squaring!(y::$T{T}, x::$T{T}, aux::$T{T},
             p::Integer) where {T<:NumberNotSeries}
@@ -389,7 +390,6 @@ end
 
     if ordT == lnull
         if isinteger(r)
-            # TODO: get rid of allocations here
             power_by_squaring!(res[ordT], a[l0], aux[0], round(Int,r))
             return nothing
         end

From 8610e211500c6574f7a667955f54504ec5a1972d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Sun, 14 Jul 2024 12:48:45 +0200
Subject: [PATCH 19/26] Bump patch version

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index e25d96af..a6dd8c1d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "TaylorSeries"
 uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
-version = "0.17.8"
+version = "0.17.9"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

From ba30e7a8596d02783b9de0fb402f34dd9350ce63 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Wed, 17 Jul 2024 16:02:20 +0200
Subject: [PATCH 20/26] Switch back from dispatch-by-value to if/else for pow!

---
 src/power.jl | 130 +++++++++++++++++++++------------------------------
 1 file changed, 52 insertions(+), 78 deletions(-)

diff --git a/src/power.jl b/src/power.jl
index b53a3d53..eb5e87b7 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -99,26 +99,23 @@ end
 
 # power_by_squaring; slightly modified from base/intfuncs.jl
 # Licensed under MIT "Expat"
-for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
-    @eval begin
-        power_by_squaring(x::$T, p::Integer) = power_by_squaring(x, Val(p))
-        power_by_squaring(x::$T, ::Val{0}) = one(x)
-        power_by_squaring(x::$T, ::Val{1}) = copy(x)
-        power_by_squaring(x::$T, ::Val{2}) = square(x)
-        power_by_squaring(x::$T, ::Val{3}) = x*square(x)
-        if $T != HomogeneousPolynomial
-            function power_by_squaring(x::$T, ::Val{P}) where P
-                p = P # copy static parameter `P` into local variable `p`
-                y = zero(x)
-                aux = zero(x)
-                power_by_squaring!(y, x, aux, p)
-                return y
-            end
-        end
+for T in (:Taylor1, :TaylorN)
+    @eval function power_by_squaring(x::$T, p::Integer)
+        (p == 0) && return one(x)
+        (p == 1) && return copy(x)
+        (p == 2) && return square(x)
+        (p == 3) && return x*square(x)
+        y = zero(x)
+        aux = zero(x)
+        power_by_squaring!(y, x, aux, p)
+        return y
     end
 end
-function power_by_squaring(x::HomogeneousPolynomial, ::Val{P}) where P
-    p = P # copy static parameter `P` into local variable `p`
+function power_by_squaring(x::HomogeneousPolynomial, p::Integer)
+    (p == 0) && return one(x)
+    (p == 1) && return copy(x)
+    (p == 2) && return square(x)
+    (p == 3) && return x*square(x)
     t = trailing_zeros(p) + 1
     p >>= t
     while (t -= 1) > 0
@@ -137,12 +134,7 @@ function power_by_squaring(x::HomogeneousPolynomial, ::Val{P}) where P
 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)
-power_by_squaring(x::TaylorN{Taylor1{T}}, ::Val{3}) where {T<:NumberNotSeries} = x*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`
+function power_by_squaring(x::TaylorN{Taylor1{T}}, p::Integer) where {T<:NumberNotSeries}
     t = trailing_zeros(p) + 1
     p >>= t
     while (t -= 1) > 0
@@ -257,24 +249,18 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`.
 
 """ pow!
 
-# pow! main dispatcher
-@inline pow!(c::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}, r::S, k::Int) where
-        {T<:Number, S <: Real} = pow!(c, a, b, Val(r), k)
-# pow! dispatches by value
-@inline pow!(c::Taylor1{T}, a::Taylor1{T}, ::Taylor1{T}, ::Val{0}, k::Int) where
-        {T<:Number} = one!(c, a, k)
-@inline pow!(c::Taylor1{T}, a::Taylor1{T}, ::Taylor1{T}, ::Val{1}, k::Int) where
-        {T<:Number} = identity!(c, a, k)
-@inline pow!(c::Taylor1{T}, a::Taylor1{T}, ::Taylor1{T}, ::Val{2}, k::Int) where
-        {T<:Number} = sqr!(c, a, k)
-@inline pow!(c::Taylor1{T}, a::Taylor1{T}, ::Taylor1{T}, ::Val{0.5}, k::Int) where
-        {T<:Number} = sqrt!(c, a, k)
-# fallback pow! method
-@inline function pow!(c::Taylor1{T}, a::Taylor1{T}, ::Taylor1{T}, ::Val{R}, k::Int) where
-        {T<:Number, R}
-
-    # copy static parameter `R` into local variable `r`
-    r = R
+@inline function pow!(c::Taylor1{T}, a::Taylor1{T}, ::Taylor1{T}, r::S, k::Int) where
+        {T<:Number, S <: Real}
+
+    if r == 0
+        return one!(c, a, k)
+    elseif r == 1
+        return identity!(c, a, k)
+    elseif r == 2
+        return sqr!(c, a, k)
+    elseif r == 0.5
+        return sqrt!(c, a, k)
+    end
 
     # Sanity
     zero!(c, k)
@@ -312,24 +298,18 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`.
     return nothing
 end
 
-# pow! main dispatcher
-@inline pow!(c::TaylorN{T}, a::TaylorN{T}, b::TaylorN{T}, r::S, k::Int) where
-        {T<:NumberNotSeriesN, S <: Real} = pow!(c, a, b, Val(r), k)
-# pow! dispatches by value
-@inline pow!(c::TaylorN{T}, a::TaylorN{T}, ::TaylorN{T}, ::Val{0}, k::Int) where
-        {T<:NumberNotSeriesN} = one!(c, a, k)
-@inline pow!(c::TaylorN{T}, a::TaylorN{T}, ::TaylorN{T}, ::Val{1}, k::Int) where
-        {T<:NumberNotSeriesN} = identity!(c, a, k)
-@inline pow!(c::TaylorN{T}, a::TaylorN{T}, ::TaylorN{T}, ::Val{2}, k::Int) where
-        {T<:NumberNotSeriesN} = sqr!(c, a, k)
-@inline pow!(c::TaylorN{T}, a::TaylorN{T}, ::TaylorN{T}, ::Val{0.5}, k::Int) where
-        {T<:NumberNotSeriesN} = sqrt!(c, a, k)
-# fallback pow! method
-@inline function pow!(c::TaylorN{T}, a::TaylorN{T}, ::TaylorN{T}, ::Val{R}, k::Int) where
-        {T<:NumberNotSeriesN, R}
-
-    # copy static parameter `R` into local variable `r`
-    r = R
+@inline function pow!(c::TaylorN{T}, a::TaylorN{T}, ::TaylorN{T}, r::S, k::Int) where
+        {T<:NumberNotSeriesN, S<:Real}
+
+        if r == 0
+            return one!(c, a, k)
+        elseif r == 1
+            return identity!(c, a, k)
+        elseif r == 2
+            return sqr!(c, a, k)
+        elseif r == 0.5
+            return sqrt!(c, a, k)
+        end
 
     if k == 0
         @inbounds c[0][1] = ( constant_term(a) )^r
@@ -351,24 +331,18 @@ end
     return nothing
 end
 
-# pow! main dispatcher for `Taylor1{TaylorN{T}}`
-@inline pow!(c::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}}, r::S, k::Int) where
-        {T<:NumberNotSeries, S <: Real} = pow!(c, a, b, Val(r), k)
-# pow! dispatches by value for `Taylor1{TaylorN{T}}`
-@inline pow!(c::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, ::Taylor1{TaylorN{T}}, ::Val{0}, k::Int) where
-        {T<:NumberNotSeries} = one!(c, a, k)
-@inline pow!(c::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, ::Taylor1{TaylorN{T}}, ::Val{1}, k::Int) where
-        {T<:NumberNotSeries} = identity!(c, a, k)
-@inline pow!(c::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, ::Taylor1{TaylorN{T}}, ::Val{2}, k::Int) where
-        {T<:NumberNotSeries} = sqr!(c, a, k)
-@inline pow!(c::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, ::Taylor1{TaylorN{T}}, ::Val{0.5}, k::Int) where
-        {T<:NumberNotSeries} = sqrt!(c, a, k)
-# fallback pow! method for `Taylor1{TaylorN{T}}`
-@inline function pow!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, aux::Taylor1{TaylorN{T}}, ::Val{R},
-        ordT::Int) where {T<:NumberNotSeries, R}
-
-    # copy static parameter `R` into local variable `r`
-    r = R
+@inline function pow!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, aux::Taylor1{TaylorN{T}}, r::S,
+        ordT::Int) where {T<:NumberNotSeries, S<:Real}
+
+    if r == 0
+        return one!(c, a, k)
+    elseif r == 1
+        return identity!(c, a, k)
+    elseif r == 2
+        return sqr!(c, a, k)
+    elseif r == 0.5
+        return sqrt!(c, a, k)
+    end
 
     # Sanity
     zero!(res, ordT)

From ea15252ffacaedd8b313efdc14d36caa3a014a8f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Wed, 17 Jul 2024 16:28:09 +0200
Subject: [PATCH 21/26] Revert change in setindex! overload for Taylor1

---
 src/auxiliary.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/auxiliary.jl b/src/auxiliary.jl
index 4a5e855a..938120c7 100644
--- a/src/auxiliary.jl
+++ b/src/auxiliary.jl
@@ -82,7 +82,7 @@ getindex(a::Taylor1{T}, u::StepRange{Int,Int}) where {T<:Number} =
     view(a.coeffs, u[:] .+ 1)
 
 setindex!(a::Taylor1{T}, x::T, n::Int) where {T<:Number} = a.coeffs[n+1] = x
-setindex!(a::Taylor1{T}, x::T, n::Int) where {T<:AbstractSeries} = setindex!(a.coeffs, x, n+1)
+setindex!(a::Taylor1{T}, x::T, n::Int) where {T<:AbstractSeries} = setindex!(a.coeffs, deepcopy(x), n+1)
 setindex!(a::Taylor1{T}, x::T, u::UnitRange{Int}) where {T<:Number} =
     a.coeffs[u .+ 1] .= x
 function setindex!(a::Taylor1{T}, x::Array{T,1}, u::UnitRange{Int}) where {T<:Number}

From 21c876b4a806e3a3ce1b56185f88beca5536a868 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Wed, 17 Jul 2024 20:33:30 +0200
Subject: [PATCH 22/26] Add test for nested Taylor1s setindex! method

---
 test/mixtures.jl | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/test/mixtures.jl b/test/mixtures.jl
index 0639893f..9c414bf7 100644
--- a/test/mixtures.jl
+++ b/test/mixtures.jl
@@ -574,4 +574,14 @@ end
     @test cos(to)(0.0) == cos(to[0])
     @test to(ti) == to[0] + ti
     @test evaluate(to*ti, ti) == to[0]*ti + ti^2
+
+    @testset "Test setindex! method for nested Taylor1s" begin
+        t = Taylor1(2)
+        y = one(t)
+        x = Taylor1([t,2t,t^2,0t,t^3])
+        x[3] = y
+        @test x[3] !== y
+        y[2] = -5.0
+        @test x[3][2] == 0.0
+    end
 end

From 702fbac004d28c4765c7330133cacad8e3ca29c9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Wed, 17 Jul 2024 20:54:51 +0200
Subject: [PATCH 23/26] De-bump patch version

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index a6dd8c1d..e25d96af 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "TaylorSeries"
 uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
-version = "0.17.9"
+version = "0.17.8"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

From 6264eae73889a312612ba3848294855c98a2c294 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Wed, 17 Jul 2024 21:06:45 +0200
Subject: [PATCH 24/26] Another approach to suggestion by @lbenet

---
 src/auxiliary.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/auxiliary.jl b/src/auxiliary.jl
index 938120c7..534b1054 100644
--- a/src/auxiliary.jl
+++ b/src/auxiliary.jl
@@ -17,9 +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
         @simd for ord in lencoef+1:order+1
-            @inbounds coeffs[ord] = zero(coeffs[1])
+            @inbounds coeffs[ord] = zero(c1)
         end
     end
     return nothing
@@ -38,8 +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)
+    c1 = coeffs[1]
     @simd for ord in lencoef+1:num_coeffs
-        @inbounds coeffs[ord] = zero(coeffs[1])
+        @inbounds coeffs[ord] = zero(c1)
     end
     return nothing
 end

From 4f613d63313151b5e9282a4d3c3faf6d455c9193 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Wed, 17 Jul 2024 21:26:30 +0200
Subject: [PATCH 25/26] Fix typo

---
 src/power.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/power.jl b/src/power.jl
index eb5e87b7..74986b73 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -335,13 +335,13 @@ end
         ordT::Int) where {T<:NumberNotSeries, S<:Real}
 
     if r == 0
-        return one!(c, a, k)
+        return one!(res, a, k)
     elseif r == 1
-        return identity!(c, a, k)
+        return identity!(res, a, k)
     elseif r == 2
-        return sqr!(c, a, k)
+        return sqr!(res, a, k)
     elseif r == 0.5
-        return sqrt!(c, a, k)
+        return sqrt!(res, a, k)
     end
 
     # Sanity

From ad1805730673f68cc8d96d5e428b2d26bacdb8b0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Wed, 17 Jul 2024 21:28:29 +0200
Subject: [PATCH 26/26] Fix another typo

---
 src/power.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/power.jl b/src/power.jl
index 74986b73..08ec1f6e 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -335,13 +335,13 @@ end
         ordT::Int) where {T<:NumberNotSeries, S<:Real}
 
     if r == 0
-        return one!(res, a, k)
+        return one!(res, a, ordT)
     elseif r == 1
-        return identity!(res, a, k)
+        return identity!(res, a, ordT)
     elseif r == 2
-        return sqr!(res, a, k)
+        return sqr!(res, a, ordT)
     elseif r == 0.5
-        return sqrt!(res, a, k)
+        return sqrt!(res, a, ordT)
     end
 
     # Sanity