From 7879ceba0bc22854f0d26cd0ba5eaf8e9982e1e2 Mon Sep 17 00:00:00 2001
From: Luis Benet <benet@icf.unam.mx>
Date: Wed, 28 Aug 2024 16:47:44 -0600
Subject: [PATCH] More fixes

---
 ext/TaylorSeriesIAExt.jl | 19 +++++++++++--------
 src/power.jl             | 32 ++++++++++++++++----------------
 2 files changed, 27 insertions(+), 24 deletions(-)

diff --git a/ext/TaylorSeriesIAExt.jl b/ext/TaylorSeriesIAExt.jl
index 36ad3ad8..a6b541f3 100644
--- a/ext/TaylorSeriesIAExt.jl
+++ b/ext/TaylorSeriesIAExt.jl
@@ -9,22 +9,25 @@ import TaylorSeries: evaluate, _evaluate, normalize_taylor, square
 isdefined(Base, :get_extension) ? (using IntervalArithmetic) : (using ..IntervalArithmetic)
 
 
-# TS._pow
+# _pow
 for T in (:Taylor1, :TaylorN)
     @eval begin
-        # Use power_by_squaring!
+        # Uses TS.power_by_squaring! through `pow!`
         function TS._pow(a::$T{Interval{S}}, n::Integer) where {S<:Real}
-            aux = one(constant_term(a))^n
+            a0 = constant_term(a)
+            aux = one(a0)^n
+            a[0] = aux * a0
             c = $T( zero(aux), a.order)
             caux = zero(c)
-            TS.power_by_squaring!(c, a, caux, n)
+            # TS.power_by_squaring!(c, a, caux, n)
+            TS.pow!(c, a, caux, n, 1)
             return c
         end
     end
 end
 
 function TS._pow(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real}
-    isinteger(r) && r >= 0 && return a^Int(r)
+    isinteger(r) && r >= 0 && return TS._pow(a, Int(r))
     a0 = constant_term(a) ∩ Interval(zero(T), T(Inf))
     aux = one(a0)^r
     a[0] = aux * a0
@@ -34,7 +37,7 @@ function TS._pow(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real}
     (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)
-    aux0 = deepcopy(c)
+    aux0 = zero(c)
     for k in eachindex(c)
         TS.pow!(c, a, aux0, r, k)
     end
@@ -42,7 +45,7 @@ function TS._pow(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real}
 end
 
 function TS._pow(a::TaylorN{Interval{T}}, r::S) where {T<:Real, S<:Real}
-    isinteger(r) && r >= 0 && return a^Int(r)
+    isinteger(r) && r >= 0 && return TS._pow(a, Int(r))
     a0 = constant_term(a) ∩ Interval(zero(T), T(Inf))
     aux = one(a0^r)
     a[0] = aux * a0
@@ -59,7 +62,7 @@ end
 for T in (:Taylor1, :TaylorN)
     @eval function log(a::$T{Interval{S}}) where {S<:Real}
         iszero(constant_term(a)) && throw(DomainError(a,
-                """The 0-th order coefficient must be non-zero in order to expand `log` around 0."""))
+            """The 0-th order coefficient must be non-zero in order to expand `log` around 0."""))
         a0 = constant_term(a) ∩ Interval(S(0.0), S(Inf))
         order = a.order
         aux = log(a0)
diff --git a/src/power.jl b/src/power.jl
index f0f82763..901bc59b 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -31,11 +31,9 @@ for T in (:Taylor1, :TaylorN)
     @eval ^(a::$T, z::T) where {T<:Complex} = exp( z*log(a) )
 end
 
-^(a::Taylor1{TaylorN{T}}, n::Integer) where {T<:NumberNotSeries} =
-    a^float(n)
+^(a::Taylor1{TaylorN{T}}, n::Integer) where {T<:NumberNotSeries} = a^float(n)
 
-^(a::Taylor1{TaylorN{T}}, r::Rational) where {T<:NumberNotSeries} =
-    a^float(r)
+^(a::Taylor1{TaylorN{T}}, r::Rational) where {T<:NumberNotSeries} = a^float(r)
 
 
 ## Real power ##
@@ -87,7 +85,7 @@ function _pow(a::Taylor1{T}, r::S) where {T<:Number, S<:Real}
     (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)
-    aux0 = deepcopy(c)
+    aux0 = zero(c)
     for k in eachindex(c)
         pow!(c, a, aux0, r, k)
     end
@@ -272,22 +270,24 @@ end
 end
 
 # Uses power_by_squaring!
-@inline function pow!(res::TaylorN{T}, a::TaylorN{T}, aux::TaylorN{T},
-        r::S, k::Int) where {T<:NumberNotSeriesN, S<:Integer}
-    (r == 0) && return one!(res, a, k)
-    (r == 1) && return identity!(res, a, k)
-    (r == 2) && return sqr!(res, a, k)
-    power_by_squaring!(res, a, aux, r)
-    return nothing
+for T in (:Taylor1, :TaylorN)
+    @eval @inline function pow!(res::$T{T}, a::$T{T}, aux::$T{T},
+            r::S, k::Int) where {T<:NumberNotSeriesN, S<:Integer}
+        (r == 0) && return one!(res, a, k)
+        (r == 1) && return identity!(res, a, k)
+        (r == 2) && return sqr!(res, a, k)
+        power_by_squaring!(res, a, aux, r)
+        return nothing
+    end
 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}
-    (r == 0) && return one!(c, a, k)
-    (r == 1) && return identity!(c, a, k)
-    (r == 2) && return sqr!(c, a, k)
-    (r == 0.5) && return sqrt!(c, a, k)
+    (r == 0) && return one!(res, a, ordT)
+    (r == 1) && return identity!(res, a, ordT)
+    (r == 2) && return sqr!(res, a, ordT)
+    (r == 0.5) && return sqrt!(res, a, ordT)
     # Sanity
     zero!(res, ordT)
     # First non-zero coefficient