diff --git a/ext/TaylorSeriesIAExt.jl b/ext/TaylorSeriesIAExt.jl
index 30ecc7c4..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, 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
diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 2af3e8eb..4c8c0a9a 100644
--- a/src/arithmetic.jl
+++ b/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
diff --git a/src/auxiliary.jl b/src/auxiliary.jl
index 22c1255c..534b1054 100644
--- a/src/auxiliary.jl
+++ b/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
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/dictmutfunct.jl b/src/dictmutfunct.jl
index d5ec11cf..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, :_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
-
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/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..08ec1f6e 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -58,37 +58,99 @@ 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)
-        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}
         t = trailing_zeros(p) + 1
         p >>= t
+        # aux = x
+        for k in eachindex(aux)
+            identity!(aux, x, k)
+        end
         while (t -= 1) > 0
-            x = square(x)
+            # aux = square(aux)
+            for k in reverse(eachindex(aux))
+                sqr!(aux, k)
+            end
+        end
+        # y = aux
+        for k in eachindex(y)
+            identity!(y, aux, k)
         end
-        y = x
         while p > 0
             t = trailing_zeros(p) + 1
             p >>= t
             while (t -= 1) ≥ 0
-                x = square(x)
+                # aux = square(aux)
+                for k in reverse(eachindex(aux))
+                    sqr!(aux, k)
+                end
             end
-            y *= x
+            # y = y * aux
+            mul!(y, aux)
         end
+        return nothing
+    end
+end
+
+
+# power_by_squaring; slightly modified from base/intfuncs.jl
+# Licensed under MIT "Expat"
+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, 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
+        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
+
 
+function power_by_squaring(x::TaylorN{Taylor1{T}}, p::Integer) where {T<:NumberNotSeries}
+    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}
@@ -107,8 +169,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)
     for k in eachindex(c)
-        pow!(c, aa, r, k)
+        pow!(c, aa, aux0, r, k)
     end
 
     return c
@@ -132,8 +195,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
@@ -158,8 +222,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)
     for k in eachindex(c)
-        pow!(c, aa, r, k)
+        pow!(c, aa, aux0, r, k)
     end
 
     return c
@@ -184,8 +249,8 @@ 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
-        {T<:Number, S<:Real}
+@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)
@@ -233,18 +298,18 @@ 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
-        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 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
@@ -266,11 +331,11 @@ 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::Taylor1{TaylorN{T}}, r::S,
         ordT::Int) where {T<:NumberNotSeries, S<:Real}
 
     if r == 0
-        return one!(res, ordT)
+        return one!(res, a, ordT)
     elseif r == 1
         return identity!(res, a, ordT)
     elseif r == 2
@@ -299,8 +364,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
+            power_by_squaring!(res[ordT], a[l0], aux[0], round(Int,r))
             return nothing
         end
 
@@ -310,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[0], r, ordQ)
         end
 
         return nothing
@@ -339,7 +403,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 +511,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
+                (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) && ( @inbounds 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
 
diff --git a/test/manyvariables.jl b/test/manyvariables.jl
index 1b7f243b..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
@@ -507,22 +527,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/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
diff --git a/test/mutatingfuncts.jl b/test/mutatingfuncts.jl
index 6599391f..fb9be7ef 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)
+    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 120a5a9f..847ff8b3 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 = tt
+    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 = tt
+    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
@@ -616,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])