diff --git a/ext/TaylorSeriesIAExt.jl b/ext/TaylorSeriesIAExt.jl
index d526bd16..36ad3ad8 100644
--- a/ext/TaylorSeriesIAExt.jl
+++ b/ext/TaylorSeriesIAExt.jl
@@ -8,67 +8,50 @@ import TaylorSeries: evaluate, _evaluate, normalize_taylor, square
 
 isdefined(Base, :get_extension) ? (using IntervalArithmetic) : (using ..IntervalArithmetic)
 
-# Method used for Taylor1{Interval{T}}^n
+
+# TS._pow
 for T in (:Taylor1, :TaylorN)
     @eval begin
-        function ^(a::$T{Interval{S}}, n::Integer) where {S<:Real}
-            n == 0 && return one(a)
-            n == 1 && return copy(a)
-            n == 2 && return square(a)
-            n < 0 && return a^float(n)
-            return power_by_squaring(a, n)
+        # Use power_by_squaring!
+        function TS._pow(a::$T{Interval{S}}, n::Integer) where {S<:Real}
+            aux = one(constant_term(a))^n
+            c = $T( zero(aux), a.order)
+            caux = zero(c)
+            TS.power_by_squaring!(c, a, caux, n)
+            return c
         end
-        ^(a::$T{Interval{S}}, r::Rational) where {S<:Real} = a^float(r)
     end
 end
 
-function ^(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real}
+function TS._pow(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real}
+    isinteger(r) && r >= 0 && return a^Int(r)
     a0 = constant_term(a) ∩ Interval(zero(T), T(Inf))
     aux = one(a0)^r
-
-    iszero(r) && return Taylor1(aux, a.order)
-    aa = one(aux) * a
-    aa[0] = one(aux) * a0
-    r == 1 && return aa
-    r == 2 && return square(aa)
-    r == 1/2 && return sqrt(aa)
-
+    a[0] = aux * a0
+    #
     l0 = findfirst(a)
     lnull = trunc(Int, r*l0 )
-    if (a.order-lnull < 0) || (lnull > a.order)
-        return Taylor1( zero(aux), a.order)
-    end
-    c_order = l0 == 0 ? a.order : min(a.order, trunc(Int,r*a.order))
+    (lnull > a.order) && return Taylor1( zero(aux), a.order)
+    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, r, k)
+    aux0 = deepcopy(c)
+    for k in eachindex(c)
+        TS.pow!(c, a, aux0, r, k)
     end
-
     return c
 end
-function ^(a::TaylorN{Interval{T}}, r::S) where {T<:Real, S<:Real}
+
+function TS._pow(a::TaylorN{Interval{T}}, r::S) where {T<:Real, S<:Real}
+    isinteger(r) && r >= 0 && return a^Int(r)
     a0 = constant_term(a) ∩ Interval(zero(T), T(Inf))
-    a0r = a0^r
-    aux = one(a0r)
-
-    iszero(r) && return TaylorN(aux, a.order)
-    aa = aux * a
-    aa[0] = aux * a0
-    r == 1 && return aa
-    r == 2 && return square(aa)
-    r == 1/2 && return sqrt(aa)
-    isinteger(r) && return aa^round(Int,r)
-
-    # @assert !iszero(a0)
-    iszero(a0) && throw(DomainError(a,
-        """The 0-th order TaylorN coefficient must be non-zero
-        in order to expand `^` around 0."""))
-
-    c = TaylorN( a0r, a.order)
-    for ord in 1:a.order
-        TS.pow!(c, aa, c, r, ord)
+    aux = one(a0^r)
+    a[0] = aux * a0
+    #
+    c = TaylorN( zero(aux), a.order)
+    aux0 = zero(c)
+    for ord in eachindex(c)
+        TS.pow!(c, a, aux0, r, ord)
     end
-
     return c
 end
 
diff --git a/test/intervals.jl b/test/intervals.jl
index d163acfa..fc66506f 100644
--- a/test/intervals.jl
+++ b/test/intervals.jl
@@ -9,6 +9,7 @@ using Test
 @testset "Tests Taylor1 and TaylorN expansions over Intervals" begin
     a = 1..2
     b = -1 .. 1
+    p3(x, a) = x^3 + 3*a*x^2 + 3*a^2*x + a^3
     p4(x, a) = x^4 + 4*a*x^3 + 6*a^2*x^2 + 4*a^3*x + a^4
     p5(x, a) = x^5 + 5*a*x^4 + 10*a^2*x^3 + 10*a^3*x^2 + 5*a^4*x + a^5
 
@@ -24,22 +25,34 @@ using Test
     @test normalize_taylor(ti) == ti
     @test normalize_taylor(x) == x
 
-    @test p4(ti,-a) == (ti-a)^4
-    @test p5(ti,-a) == (ti-a)^5
-    @test p4(ti,-b) == (ti-b)^4
+    @test p3(ti,-a) == (ti-a)^3 == (ti-a)^3.0
+    @test p4(ti,-a) == (ti-a)^4 == (ti-a)^4.0
+    @test p5(ti,-a) == (ti-a)^5 == (ti-a)^5.0
+    @test (ti-b)^3 == (ti-b)^3.0
+    @test all((p3(ti,-b)).coeffs .⊆ ((ti-b)^3).coeffs)
+    @test p4(ti,-b) == (ti-b)^4 == (ti-b)^4.0
+    @test (ti-b)^5 == (ti-b)^5.0
     @test all((p5(ti,-b)).coeffs .⊆ ((ti-b)^5).coeffs)
 
 
-    @test p4(x,-y) == (x-y)^4
-    @test p5(x,-y) == (x-y)^5
-    @test p4(x,-a) == (x-a)^4
-    @test p5(x,-a) == (x-a)^5
+    @test p3(x,-y) == (x-y)^3 == (x-y)^3.0
+    @test p4(x,-y) == (x-y)^4 == (x-y)^4.0
+    @test p5(x,-y) == (x-y)^5 == (x-y)^5.0
+    @test p3(x,-a) == (x-a)^3 == (x-a)^3.0
+    @test p4(x,-a) == (x-a)^4 == (x-a)^4.0
+    @test p5(x,-a) == (x-a)^5 == (x-a)^5.0
+    @test (x-b)^3 == (x-b)^3.0
+    for ind in eachindex(p3(x,-b))
+        @test all((p3(x,-b)[ind]).coeffs .⊆ (((x-b)^3)[ind]).coeffs)
+    end
     @test p4(x,-b) == (x-b)^4
+    @test (x-b)^5 == (x-b)^5.0
     for ind in eachindex(p5(x,-b))
         @test all((p5(x,-b)[ind]).coeffs .⊆ (((x-b)^5)[ind]).coeffs)
     end
 
     # Tests `evaluate`
+    @test evaluate(p3(x,y), IntervalBox(a,-b)) == p3(a, -b)
     @test evaluate(p4(x,y), IntervalBox(a,-b)) == p4(a, -b)
     @test (p5(x,y))(IntervalBox(a,b)) == p5(a, b)
     @test (a-b)^4 ⊆ ((x-y)^4)(a × b)