Separate methods for ^ to use power_by_squaring or pow!

src/power.jl

@@ -15,15 +15,35 @@ function ^(a::HomogeneousPolynomial, n::Integer)
     return power_by_squaring(a, n)
 end
 
+#= The following three methods are coded like that, to use
+preferentially a^float(n), but for cases like Taylor1{Interval{T}}^n
+power_by_squaring is used. The latter is important when the
+0-th order coefficient is/contains zero.
+=#
 function ^(a::Taylor1{T}, n::Integer) where {T<:Number}
     n == 0 && return one(a)
     n == 1 && return copy(a)
     n == 2 && return square(a)
-    # Return a^float(n) instead of power_by_squaring(a,n)
-    # because it yields better accuracy
     return a^float(n)
 end
 
+# Method used for Taylor1{Interval{T}}^n
+function ^(a::Taylor1{T}, n::Integer) where {T<: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)
+end
+
+function ^(a::Taylor1{T}, n::Integer) where {T<:AbstractFloat}
+    n == 0 && return one(a)
+    n == 1 && return copy(a)
+    n == 2 && return square(a)
+    return a^float(n)
+end
+
+
 function ^(a::TaylorN{T}, n::Integer) where {T<:Number}
     n == 0 && return one(a)
     n == 1 && return copy(a)
@@ -328,7 +348,7 @@ Return `c = a*a` with no allocation; all parameters are `HomogeneousPolynomial`.
         iszero(ca) && continue
         inda = idxTb[na]
         pos = posTb[2*inda]
-        c[pos] += ca * ca
+        c[pos] += ca^2
         @inbounds for nb = na+1:num_coeffs_a
             cb = a[nb]
             iszero(cb) && continue

test/REQUIRE

@@ -0,0 +1 @@
+IntervalArithmetic 0.13.0

test/intervals.jl

@@ -0,0 +1,41 @@
+# This file is part of TaylorSeries.jl, MIT licensed
+#
+
+using TaylorSeries, IntervalArithmetic
+
+if VERSION < v"0.7.0-DEV.2004"
+    using Base.Test
+    eeuler = Base.e
+else
+    using Test
+    eeuler = Base.MathConstants.e
+end
+
+@testset "Tests Taylor1 and TaylorN expansions over Intervals" begin
+    a = 1..2
+    b = -1 .. 1
+    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
+
+    ti = Taylor1(Interval{Float64}, 10)
+    x, y = set_variables(Interval{Float64}, "x y")
+
+    @test eltype(ti) == Interval{Float64}
+    @test eltype(x) == Interval{Float64}
+
+    @test p4(ti,-a) == (ti-a)^4
+    @test p5(ti,-a) == (ti-a)^5
+    @test p4(ti,-b) == (ti-b)^4
+    @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 p4(x,-b) == (x-b)^4
+    for ind in eachindex((p5(x,-b)).coeffs)
+        @test all(((p5(x,-b)).coeffs[ind]).coeffs .⊆ (((x-b)^5).coeffs[ind]).coeffs)
+    end
+
+end

test/onevariable.jl

@@ -145,7 +145,7 @@ end
     @test (Rational(1,2)*tsquare)[2] == 1//2
     @test t^2/tsquare == ot
     @test ((1+t)^(1/3))[2]+1/9 ≤ tol1
-    @test 1-tsquare == (1+t)-t*(1+t)
+    @test (1.0-tsquare)^3 == (1.0-t)^3*(1.0+t)^3
     @test (1-tsquare)^2 == (1+t)^2.0 * (1-t)^2.0
     @test (sqrt(1+t))[2] == -1/8
     @test ((1-tsquare)^(1//2))^2 == 1-tsquare

test/runtests.jl

@@ -8,7 +8,8 @@ testfiles = (
     "mixtures.jl",
     "mutatingfuncts.jl",
     "identities_Euler.jl",
-    "fateman40.jl"
+    "fateman40.jl",
+    "intervals.jl"
     )
 
 for file in testfiles