From 0937596692f79f84cb55b00123117632e913112d Mon Sep 17 00:00:00 2001
From: Luis Benet <lbenet@users.noreply.github.com>
Date: Fri, 2 Oct 2020 12:22:20 -0500
Subject: [PATCH] For ^(a::Taylor1, n::Integer) fall back to a^float(n) (#250)

* Move a method oF ^ to intervals.jl

* Fall back ^(a::Taylor1, n::Integer) to a^float(n) and fix a test

This solves a subtle type unstability.

* Bump a patch version
---
 Project.toml         |  2 +-
 src/dictmutfunct.jl  |  2 +-
 src/intervals.jl     | 11 +++++++++++
 src/power.jl         | 31 ++++---------------------------
 test/broadcasting.jl | 10 +++++++---
 5 files changed, 24 insertions(+), 32 deletions(-)

diff --git a/Project.toml b/Project.toml
index 59fcd22b..8aabf93d 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.10.7"
+version = "0.10.8"
 
 [deps]
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
diff --git a/src/dictmutfunct.jl b/src/dictmutfunct.jl
index aade056d..5b5b4b8f 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 ^ _arg2)],
+    :^ => [:pow!, (:_res, :_arg1, :_arg2, :_k), :(_res = _arg1 ^ float(_arg2))],
 );
 
 """
diff --git a/src/intervals.jl b/src/intervals.jl
index 6a8c3f7d..9a28afa8 100644
--- a/src/intervals.jl
+++ b/src/intervals.jl
@@ -1,5 +1,16 @@
 using .IntervalArithmetic
 
+# Method used for Taylor1{Interval{T}}^n
+function ^(a::Taylor1{T}, n::Integer) where {T<:Interval}
+    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 evaluate(a::Taylor1, dx::Interval)
     order = a.order
     uno = one(dx)
diff --git a/src/power.jl b/src/power.jl
index 52a6636b..2f7591dd 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -15,34 +15,11 @@ 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.
+#= The following method computes `a^float(n)` (except for cases like 
+Taylor1{Interval{T}}^n, where `power_by_squaring` is used), to 
+use internally `pow!`.
 =#
-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)
-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
-
+^(a::Taylor1, n::Integer) = a^float(n)
 
 function ^(a::TaylorN{T}, n::Integer) where {T<:Number}
     n == 0 && return one(a)
diff --git a/test/broadcasting.jl b/test/broadcasting.jl
index b956d39a..b559a1d1 100644
--- a/test/broadcasting.jl
+++ b/test/broadcasting.jl
@@ -42,7 +42,9 @@ using Test
     @. ts = 3 * t^2 - 1
     @test ts == 3 * t^2 - 1
 
-    tt = Taylor1([zero(t), one(t)], 2)
+    #  `tt` has to be `Taylor1{Taylor1{Float64}}` (instead of `Taylor1{Taylor1{Int}}`)
+    # since the method a^n (n integer) is equivalent to `a^float(n).`
+    tt = Taylor1([zero(1.0*t), one(t)], 2)
     tts = zero(tt)
     @test tt .== tt
     @. tts = 3 * tt^2 - 1
@@ -51,8 +53,10 @@ using Test
     ttt = Taylor1([zero(tt), one(tt)])
     ttts = zero(ttt)
     @test ttt .≈ ttt
-    @. ttts = 3 * ttt^2 - 1
-    @test ttts == -1
+    @. ttts = 3 * ttt^1 - 1
+    @test ttts == 3 * ttt^1 - 1
+    @. ttts = 3 * ttt^3 - 1
+    @test ttts == - 1.0
 end
 
 @testset "Broadcasting with HomogeneousPolynomial and TaylorN" begin