From d463b4c9812aaafabaa34eef4a4f7c8f234ee35b Mon Sep 17 00:00:00 2001
From: Luis Benet <benet@icf.unam.mx>
Date: Thu, 22 Aug 2024 17:35:02 -0600
Subject: [PATCH 1/8] Fixes in evaluate, and rename one method evaluate! ->
 _evaluate!

---
 src/evaluate.jl       | 107 +++++++++++++++++++-----------------------
 test/manyvariables.jl |   4 +-
 2 files changed, 50 insertions(+), 61 deletions(-)

diff --git a/src/evaluate.jl b/src/evaluate.jl
index a496c607..6e1f1ccd 100644
--- a/src/evaluate.jl
+++ b/src/evaluate.jl
@@ -254,39 +254,32 @@ Note that the syntax `a(vals)` is equivalent to
 `evaluate(a)`; use a(b::Bool, x) corresponds to
 evaluate(a, x, sorting=b).
 """
-function evaluate(a::TaylorN, vals::NTuple{N,<:Number};
-        sorting::Bool=true) where {N}
+function evaluate(a::TaylorN, vals::NTuple{N,<:Number}; sorting::Bool=true) where {N}
     @assert get_numvars() == N
     return _evaluate(a, vals, Val(sorting))
 end
 
-function evaluate(a::TaylorN, vals::NTuple{N,<:AbstractSeries};
-        sorting::Bool=false) where {N}
+function evaluate(a::TaylorN, vals::NTuple{N,<:AbstractSeries}; sorting::Bool=false) where {N}
     @assert get_numvars() == N
     return _evaluate(a, vals, Val(sorting))
 end
 
 evaluate(a::TaylorN{T}, vals::AbstractVector{<:Number};
-        sorting::Bool=true) where {T<:NumberNotSeries} =
-    evaluate(a, (vals...,); sorting=sorting)
+    sorting::Bool=true) where {T<:NumberNotSeries} = evaluate(a, (vals...,); sorting=sorting)
 
 evaluate(a::TaylorN{T}, vals::AbstractVector{<:AbstractSeries};
-        sorting::Bool=false) where {T<:NumberNotSeries} =
-    evaluate(a, (vals...,); sorting=sorting)
+    sorting::Bool=false) where {T<:NumberNotSeries} = evaluate(a, (vals...,); sorting=sorting)
 
 evaluate(a::TaylorN{Taylor1{T}}, vals::AbstractVector{S};
-        sorting::Bool=false) where {T, S} =
-    evaluate(a, (vals...,); sorting=sorting)
+    sorting::Bool=false) where {T, S} = evaluate(a, (vals...,); sorting=sorting)
 
-function evaluate(a::TaylorN{T}, s::Symbol, val::S) where
-        {T<:Number, S<:NumberNotSeriesN}
+function evaluate(a::TaylorN{T}, s::Symbol, val::S) where {T<:Number, S<:NumberNotSeriesN}
     ind = lookupvar(s)
     @assert (1 ≤ ind ≤ get_numvars()) "Symbol is not a TaylorN variable; see `get_variable_names()`"
     return evaluate(a, ind, val)
 end
 
-function evaluate(a::TaylorN{T}, ind::Int, val::S) where
-        {T<:Number, S<:NumberNotSeriesN}
+function evaluate(a::TaylorN{T}, ind::Int, val::S) where {T<:Number, S<:NumberNotSeriesN}
     @assert (1 ≤ ind ≤ get_numvars()) "Invalid `ind`; it must be between 1 and `get_numvars()`"
     R = promote_type(T,S)
     return _evaluate(convert(TaylorN{R}, a), ind, convert(R, val))
@@ -305,8 +298,7 @@ function evaluate(a::TaylorN{T}, ind::Int, val::TaylorN) where {T<:Number}
     return _evaluate(a, ind, val)
 end
 
-evaluate(a::TaylorN{T}, x::Pair{Symbol,S}) where {T, S} =
-    evaluate(a, first(x), last(x))
+evaluate(a::TaylorN{T}, x::Pair{Symbol,S}) where {T, S} = evaluate(a, first(x), last(x))
 
 evaluate(a::TaylorN{T}) where {T<:Number} = constant_term(a)
 
@@ -338,14 +330,6 @@ function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:Number}) where {N,T<:Number}
     return suma
 end
 
-function _evaluate!(res::Vector{TaylorN{T}}, a::TaylorN{T}, vals::NTuple{N,<:TaylorN},
-        valscache::Vector{TaylorN{T}}, aux::TaylorN{T}) where {N,T<:Number}
-    @inbounds for homPol in eachindex(a)
-        _evaluate!(res[homPol+1], a[homPol], vals, valscache, aux)
-    end
-    return nothing
-end
-
 function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:TaylorN}) where {N,T<:Number}
     R = promote_type(T, TS.numtype(vals[1]))
     suma = [TaylorN(zero(R), vals[1].order) for _ in eachindex(a)]
@@ -375,27 +359,37 @@ function _evaluate(a::TaylorN{T}, ind::Int, val::TaylorN{T}) where {T<:NumberNot
     return suma
 end
 
+function _evaluate!(res::Vector{TaylorN{T}}, a::TaylorN{T}, vals::NTuple{N,<:TaylorN},
+        valscache::Vector{TaylorN{T}}, aux::TaylorN{T}) where {N,T<:Number}
+    @inbounds for homPol in eachindex(a)
+        _evaluate!(res[homPol+1], a[homPol], vals, valscache, aux)
+    end
+    return nothing
+end
+
 function _evaluate!(suma::TaylorN{T}, a::HomogeneousPolynomial{T}, ind::Int, val::T) where
         {T<:NumberNotSeriesN}
     order = a.order
+    orderTN = get_order()
     if order == 0
-        suma[0] = a[1]*one(val)
+        suma[0][1] = a[1]*one(val)
         return nothing
     end
     vv = val .^ (0:order)
-    # ct = @isonethread coeff_table[order+1]
-    ct = deepcopy(coeff_table[order+1])
+    vct = zero(coeff_table[order+1][1])
+    zct = zero(coeff_table[order+1][1])
     for (i, a_coeff) in enumerate(a.coeffs)
         iszero(a_coeff) && continue
-        if ct[i][ind] == 0
+        vpow = coeff_table[order+1][i][ind]
+        if vpow == 0
             suma[order][i] += a_coeff
             continue
         end
-        vpow = ct[i][ind]
+        vct .= coeff_table[order+1][i]
+        zct[ind] = vpow
         red_order = order - vpow
-        ct[i][ind] -= vpow
-        kdic = in_base(get_order(), ct[i])
-        ct[i][ind] += vpow
+        kdic = in_base(orderTN, vct - zct)
+        zct[ind] = 0
         pos = pos_table[red_order+1][kdic]
         suma[red_order][pos] += a_coeff * vv[vpow+1]
     end
@@ -406,32 +400,33 @@ function _evaluate!(suma::TaylorN{T}, a::HomogeneousPolynomial{T}, ind::Int,
         val::TaylorN{T}, aux::TaylorN{T}) where {T<:NumberNotSeriesN}
     order = a.order
     if order == 0
-        suma[0] = a[1]
+        suma[0][1] = a[1]
         return nothing
     end
     vv = zero(suma)
-    ct = coeff_table[order+1]
+    vvaux = zero(vv)
     za = zero(a)
     for (i, a_coeff) in enumerate(a.coeffs)
         iszero(a_coeff) && continue
-        if ct[i][ind] == 0
+        vpow = coeff_table[order+1][i][ind]
+        if vpow == 0
             suma[order][i] += a_coeff
             continue
         end
-        za[i] = a_coeff
-        zero!(aux)
-        _evaluate!(aux, za, ind, one(T))
-        za[i] = zero(T)
-        vpow = ct[i][ind]
         # vv = val ^ vpow
         if constant_term(val) == 0
-            vv = val ^ vpow
+            zero!(vvaux)
+            power_by_squaring!(vv, val, vvaux, vpow)
         else
             for ordQ in eachindex(val)
                 zero!(vv, ordQ)
-                pow!(vv, val, vv, vpow, ordQ)
+                pow!(vv, val, vvaux, vpow, ordQ)
             end
         end
+        za[i] = a_coeff
+        zero!(aux)
+        _evaluate!(aux, za, ind, one(T))
+        za[i] = zero(a_coeff)
         for ordQ in eachindex(suma)
             mul!(suma, vv, aux, ordQ)
         end
@@ -456,7 +451,7 @@ evaluate(A::AbstractArray{TaylorN{T}}) where {T<:Number} = evaluate.(A)
 #function-like behavior for TaylorN
 (p::TaylorN)(x) = evaluate(p, x)
 (p::TaylorN)() = evaluate(p)
-(p::TaylorN)(s::S, x) where {S<:Union{Symbol, Int}}= evaluate(p, s, x)
+(p::TaylorN)(s::S, x) where {S<:Union{Symbol, Int}} = evaluate(p, s, x)
 (p::TaylorN)(x::Pair) = evaluate(p, first(x), last(x))
 (p::TaylorN)(x, v::Vararg{T}) where {T} = evaluate(p, (x, v...,))
 (p::TaylorN)(b::Bool, x) = evaluate(p, x, sorting=b)
@@ -470,39 +465,33 @@ evaluate(A::AbstractArray{TaylorN{T}}) where {T<:Number} = evaluate.(A)
 
 
 """
-    evaluate!(x, δt, x0)
+    evaluate!(x, δt, dest)
 
 Evaluates each element of `x::AbstractArray{Taylor1{T}}`,
 representing the Taylor expansion for the dependent variables
-of an ODE at *time* `δt`. It updates the vector `x0` with the
+of an ODE at *time* `δt`. It updates the vector `dest` with the
 computed values.
 """
 function evaluate!(x::AbstractArray{Taylor1{T}}, δt::S,
-        x0::AbstractArray{T}) where {T<:Number, S<:Number}
-    x0 .= evaluate.( x, δt )
+        dest::AbstractArray{T}) where {T<:Number, S<:Number}
+    dest .= evaluate.( x, δt )
     return nothing
 end
-# function evaluate!(x::AbstractArray{Taylor1{Taylor1{T}}}, δt::Taylor1{T},
-#         x0::AbstractArray{Taylor1{T}}) where {T<:Number}
-#     x0 .= evaluate.( x, Ref(δt) )
-#     # x0 .= evaluate.( x, δt )
-#     return nothing
-# end
 
 ## In place evaluation of multivariable arrays
 function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{T,1},
-        x0::AbstractArray{T}) where {T<:Number}
-    x0 .= evaluate.( x, Ref(δx) )
+        dest::AbstractArray{T}) where {T<:Number}
+    dest .= evaluate.( x, Ref(δx) )
     return nothing
 end
 
 function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{TaylorN{T},1},
-        x0::AbstractArray{TaylorN{T}}; sorting::Bool=true) where {T<:NumberNotSeriesN}
-    x0 .= evaluate.( x, Ref(δx), sorting = sorting)
+        dest::AbstractArray{TaylorN{T}}; sorting::Bool=true) where {T<:NumberNotSeriesN}
+    dest .= evaluate.( x, Ref(δx), sorting = sorting)
     return nothing
 end
 
-function evaluate!(a::TaylorN{T}, vals::NTuple{N,TaylorN{T}}, dest::TaylorN{T},
+function _evaluate!(a::TaylorN{T}, vals::NTuple{N,TaylorN{T}}, dest::TaylorN{T},
         valscache::Vector{TaylorN{T}}, aux::TaylorN{T}) where {N,T<:Number}
     @inbounds for homPol in eachindex(a)
         _evaluate!(dest, a[homPol], vals, valscache, aux)
@@ -518,7 +507,7 @@ function evaluate!(a::AbstractArray{TaylorN{T}}, vals::NTuple{N,TaylorN{T}},
     # loop over elements of `a`
     for i in eachindex(a)
         (!iszero(dest[i])) && zero!(dest[i])
-        evaluate!(a[i], vals, dest[i], valscache, aux)
+        _evaluate!(a[i], vals, dest[i], valscache, aux)
     end
     return nothing
 end
diff --git a/test/manyvariables.jl b/test/manyvariables.jl
index 8e49a1c0..1a09958c 100644
--- a/test/manyvariables.jl
+++ b/test/manyvariables.jl
@@ -853,9 +853,9 @@ end
         radntn!.(v)
         x1 = randn(4) .+ x
         # warmup
-        TaylorSeries.evaluate!(v, (x1...,), r)
+        evaluate!(v, (x1...,), r)
         # call twice to make sure `r` is reset on second call
-        TaylorSeries.evaluate!(v, (x1...,), r)
+        evaluate!(v, (x1...,), r)
         r2 = evaluate.(v, Ref(x1))
         @test r == r2
         @test iszero(norm(r-r2, Inf))

From 9bf3b393f4359968f09a365327a4e7eb39d7f25a Mon Sep 17 00:00:00 2001
From: Luis Benet <benet@icf.unam.mx>
Date: Fri, 23 Aug 2024 10:54:14 -0600
Subject: [PATCH 2/8] Add function barrier _pow

---
 src/power.jl          | 273 ++++++++++++++++++------------------------
 test/manyvariables.jl |   2 +-
 2 files changed, 117 insertions(+), 158 deletions(-)

diff --git a/src/power.jl b/src/power.jl
index 67ca6726..ba9982c2 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -15,62 +15,103 @@ function ^(a::HomogeneousPolynomial, n::Integer)
     return power_by_squaring(a, n)
 end
 
-#= 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!`.
-=#
-^(a::Taylor1, n::Integer) = a^float(n)
-
-function ^(a::TaylorN{T}, n::Integer) where {T<:Number}
-    n == 0 && return one(a)
-    n == 1 && return copy(a)
-    n == 2 && return square(a)
-    n < 0 && return inv( a^(-n) )
-    return power_by_squaring(a, n)
-end
-
 
 for T in (:Taylor1, :TaylorN)
-    @eval function ^(a::$T{T}, n::Integer) where {T<:Integer}
+    @eval function ^(a::$T, n::Integer)
         n == 0 && return one(a)
         n == 1 && return copy(a)
         n == 2 && return square(a)
+        return _pow(a, n)
+    end
+    @eval ^(a::$T, r::S) where {S<:Rational} = a^float(r)
+    @eval ^(a::$T, b::$T) = exp( b*log(a) )
+    @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}}, r::Rational) where {T<:NumberNotSeries} =
+    a^float(r)
+
+
+## Real power ##
+function ^(a::Taylor1{T}, r::S) where {T<:Number, S<:Real}
+    a0 = constant_term(a)
+    aux = one(a0)^r
+    iszero(r) && return Taylor1(aux, a.order)
+    aa = aux*a
+    r == 1 && return aa
+    r == 2 && return square(aa)
+    r == 0.5 && return sqrt(aa)
+    return _pow(aa, r)
+end
+
+function ^(a::TaylorN{T}, r::S) where {T<:Number, S<:Real}
+    a0 = constant_term(a)
+    aux = one(a0^r)
+    aa = aux*a
+    isinteger(r) && r ≥ 0 && return _pow(aa, round(Int, r))
+    iszero(a0) && throw(DomainError(a,
+        """The 0-th order TaylorN coefficient must be non-zero
+        in order to expand `^` around 0."""))
+    r == 0.5 && return sqrt(aa)
+    return _pow(aa, r)
+end
+
+
+# _pow
+_pow(a::Taylor1, n::Integer) = a^float(n)
+
+_pow(a::TaylorN, n::Integer) = power_by_squaring(a, n)
+
+for T in (:Taylor1, :TaylorN)
+    @eval function _pow(a::$T{T}, n::Integer) where {T<:Integer}
         n < 0 && throw(DomainError())
         return power_by_squaring(a, n)
     end
-
-    @eval function ^(a::$T{Rational{T}}, n::Integer) where {T<:Integer}
-        n == 0 && return one(a)
-        n == 1 && return copy(a)
-        n == 2 && return square(a)
+    @eval function _pow(a::$T{Rational{T}}, n::Integer) where {T<:Integer}
         n < 0 && return inv( a^(-n) )
         return power_by_squaring(a, n)
     end
+end
 
-    @eval ^(a::$T, x::S) where {S<:Rational} = a^(x.num/x.den)
-
-    @eval ^(a::$T, b::$T) = exp( b*log(a) )
-
-    @eval ^(a::$T, x::T) where {T<:Complex} = exp( x*log(a) )
+function _pow(a::Taylor1{T}, r::S) where {T<:Number, S<:Real}
+    aux = one(constant_term(a)^r)
+    l0 = findfirst(a)
+    lnull = trunc(Int, r*l0 )
+    (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)
+    for k in eachindex(c)
+        pow!(c, a, aux0, r, k)
+    end
+    return c
 end
 
-^(a::Taylor1{TaylorN{T}}, n::Integer) where {T<:NumberNotSeries} = a^float(n)
+function _pow(a::TaylorN{T}, r::S) where {T<:Number, S<:Real}
+    aux = one(constant_term(a)^r)
+    c = TaylorN( zero(aux), a.order)
+    aux0 = zero(c)
+    for ord in eachindex(a)
+        pow!(c, a, aux0, r, ord)
+    end
+    return c
+end
 
-^(a::Taylor1{TaylorN{T}}, r::Rational) where {T<:NumberNotSeries} = a^(r.num/r.den)
 
 # 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 and mixtures
 for T in (:Taylor1, :TaylorN)
-    @eval @inline function power_by_squaring_0!(y::$T{T}, x::$T{T}) where {T<:NumberNotSeries}
-        for k in eachindex(y)
-            one!(y, x, k)
+    @eval function power_by_squaring!(y::$T, x::$T, aux::$T, p::Integer)
+        if p == 0
+            for k in eachindex(y)
+                one!(y, x, k)
+            end
+            return nothing
         end
-        return nothing
-    end
-    @eval @inline function power_by_squaring!(y::$T{T}, x::$T{T}, aux::$T{T},
-            p::Integer) where {T<:NumberNotSeries}
-        (p == 0) && return power_by_squaring_0!(y, x)
         t = trailing_zeros(p) + 1
         p >>= t
         # aux = x
@@ -108,15 +149,10 @@ end
 # Licensed under MIT "Expat"
 for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
     @eval function power_by_squaring(x::$T, p::Integer)
-        if p == 0
-            return one(x)
-        elseif p == 1
-            return copy(x)
-        elseif p == 2
-            return square(x)
-        elseif p == 3
-            return x*square(x)
-        end
+        (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
@@ -135,10 +171,10 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
     end
 end
 
-# power_by_squaring specializations for non-mixed Taylor1 and TaylorN
+# power_by_squaring specializations for non-mixtures of Taylor1 and TaylorN;
 # uses internally mutating method `power_by_squaring!`
 for T in (:Taylor1, :TaylorN)
-    @eval function power_by_squaring(x::$T{T}, p::Integer) where {T <: NumberNotSeries}
+    @eval function power_by_squaring(x::$T{T}, p::Integer) where {T<:NumberNotSeries}
         (p == 0) && return one(x)
         (p == 1) && return copy(x)
         (p == 2) && return square(x)
@@ -150,83 +186,6 @@ for T in (:Taylor1, :TaylorN)
     end
 end
 
-## Real power ##
-function ^(a::Taylor1{T}, r::S) where {T<:Number, S<:Real}
-    a0 = constant_term(a)
-    aux = one(a0)^r
-
-    iszero(r) && return Taylor1(aux, a.order)
-    aa = aux*a
-    r == 1 && return aa
-    r == 2 && return square(aa)
-    r == 0.5 && return sqrt(aa)
-
-    l0 = findfirst(a)
-    lnull = trunc(Int, r*l0 )
-    (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)
-    for k in eachindex(c)
-        pow!(c, aa, aux0, r, k)
-    end
-
-    return c
-end
-
-## Real power ##
-# TODO: get rid of allocations
-function ^(a::TaylorN, r::S) where {S<:Real}
-    a0 = constant_term(a)
-    aux = one(a0^r)
-
-    iszero(r) && return TaylorN(aux, a.order)
-    aa = aux*a
-    r == 1 && return aa
-    r == 2 && return square(aa)
-    r == 0.5 && return sqrt(aa)
-    isinteger(r) && return aa^round(Int,r) # uses power_by_squaring
-
-    iszero(a0) && throw(DomainError(a,
-        """The 0-th order TaylorN coefficient must be non-zero
-        in order to expand `^` around 0."""))
-
-    c = TaylorN( zero(aux), a.order)
-    aux = deepcopy(c)
-    for ord in eachindex(a)
-        pow!(c, aa, aux, r, ord)
-    end
-
-    return c
-end
-
-function ^(a::Taylor1{TaylorN{T}}, r::S) where {T<:NumberNotSeries, S<:Real}
-    a0 = constant_term(a)
-    aux = one(a0)^r
-
-    iszero(r) && return Taylor1(aux, a.order)
-    aa = aux*a
-    r == 1 && return aa
-    r == 2 && return square(aa)
-    r == 0.5 && return sqrt(aa)
-    # Is the following needed?
-    # isinteger(r) && r > 0 && iszero(constant_term(a[0])) &&
-    #     return power_by_squaring(aa, round(Int,r))
-
-    l0 = findfirst(a)
-    lnull = trunc(Int, r*l0 )
-    (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)
-    for k in eachindex(c)
-        pow!(c, aa, aux0, r, k)
-    end
-
-    return c
-end
 
 
 # Homogeneous coefficients for real power
@@ -247,18 +206,13 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`.
 
 """ pow!
 
-@inline function pow!(c::Taylor1{T}, a::Taylor1{T}, ::Taylor1{T}, r::S, k::Int) where
-        {T<:Number, S <: Real}
+@inline function pow!(c::Taylor1{T}, a::Taylor1{T}, aux::Taylor1{T},
+                      r::S, k::Int) where {T<:Number, 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
+    (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)
 
     # Sanity
     zero!(c, k)
@@ -296,18 +250,13 @@ 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}, ::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
+@inline function pow!(c::TaylorN{T}, a::TaylorN{T}, aux::TaylorN{T},
+                      r::S, k::Int) where {T<:NumberNotSeriesN, 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)
 
     if k == 0
         @inbounds c[0][1] = ( constant_term(a) )^r
@@ -329,18 +278,13 @@ end
     return nothing
 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}
+@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, a, ordT)
-    elseif r == 1
-        return identity!(res, a, ordT)
-    elseif r == 2
-        return sqr!(res, a, ordT)
-    elseif r == 0.5
-        return sqrt!(res, a, ordT)
-    end
+    (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)
 
     # Sanity
     zero!(res, ordT)
@@ -391,6 +335,21 @@ end
     return nothing
 end
 
+for T in (:Taylor1, :TaylorN)
+    @eval begin
+        @inline function pow!(res::$T{T}, a::$T{T}, aux::$T{T},
+                              r::Integer, k::Int) where {T<:NumberNotSeries}
+            (r == 0) && return one!(res, a, k)
+            (r == 1) && return identity!(res, a, k)
+            (r == 2) && return sqr!(res, a, k)
+            # Sanity
+            zero!(res, k)
+            power_by_squaring!(res, a, aux, r)
+            return nothing
+        end
+    end
+end
+
 
 ## Square ##
 """
diff --git a/test/manyvariables.jl b/test/manyvariables.jl
index 1a09958c..ebb4b24d 100644
--- a/test/manyvariables.jl
+++ b/test/manyvariables.jl
@@ -452,7 +452,7 @@ end
     @test_throws AssertionError q[end-2:2:end] = pol.coeffs[end-1:2:end]
 
     @test_throws AssertionError yT^(-2)
-    @test_throws AssertionError yT^(-2.0)
+    @test_throws DomainError yT^(-2.0)
     @test (1+xT)^(3//2) == ((1+xT)^0.5)^3
     @test real(xH) == xH
     @test imag(xH) == zero(xH)

From 05a9a837a6ccf3b2f53ea8216c1e81f4e8620c6b Mon Sep 17 00:00:00 2001
From: Luis Benet <benet@icf.unam.mx>
Date: Fri, 23 Aug 2024 13:10:01 -0600
Subject: [PATCH 3/8] Fixes in pow!

---
 src/power.jl | 67 ++++++++++++++++++----------------------------------
 1 file changed, 23 insertions(+), 44 deletions(-)

diff --git a/src/power.jl b/src/power.jl
index ba9982c2..f0f82763 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -23,8 +23,11 @@ for T in (:Taylor1, :TaylorN)
         n == 2 && return square(a)
         return _pow(a, n)
     end
+
     @eval ^(a::$T, r::S) where {S<:Rational} = a^float(r)
+
     @eval ^(a::$T, b::$T) = exp( b*log(a) )
+
     @eval ^(a::$T, z::T) where {T<:Complex} = exp( z*log(a) )
 end
 
@@ -70,6 +73,7 @@ for T in (:Taylor1, :TaylorN)
         n < 0 && throw(DomainError())
         return power_by_squaring(a, n)
     end
+
     @eval function _pow(a::$T{Rational{T}}, n::Integer) where {T<:Integer}
         n < 0 && return inv( a^(-n) )
         return power_by_squaring(a, n)
@@ -175,6 +179,7 @@ end
 # uses internally mutating method `power_by_squaring!`
 for T in (:Taylor1, :TaylorN)
     @eval function power_by_squaring(x::$T{T}, p::Integer) where {T<:NumberNotSeries}
+        @assert p > 0
         (p == 0) && return one(x)
         (p == 1) && return copy(x)
         (p == 2) && return square(x)
@@ -208,19 +213,15 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`.
 
 @inline function pow!(c::Taylor1{T}, a::Taylor1{T}, aux::Taylor1{T},
                       r::S, k::Int) where {T<:Number, 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)
-
     # Sanity
     zero!(c, k)
-
     # First non-zero coefficient
     l0 = findfirst(a)
     l0 < 0 && return nothing
-
     # The first non-zero coefficient of the result; must be integer
     !isinteger(r*l0) && throw(DomainError(a,
         """The 0-th order Taylor1 coefficient must be non-zero
@@ -228,15 +229,13 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`.
     lnull = trunc(Int, r*l0 )
     kprime = k-lnull
     (kprime < 0 || lnull > a.order) && return nothing
-
     # Relevant for positive integer r, to avoid round-off errors
     isinteger(r) && r > 0 && (k > r*findlast(a)) && return nothing
-
+    # First non-zero coeff
     if k == lnull
         @inbounds c[k] = (a[l0])^float(r)
         return nothing
     end
-
     # The recursion formula
     if l0+kprime ≤ a.order
         @inbounds c[k] = r * kprime * c[lnull] * a[l0+kprime]
@@ -252,20 +251,15 @@ end
 
 @inline function pow!(c::TaylorN{T}, a::TaylorN{T}, aux::TaylorN{T},
                       r::S, k::Int) where {T<:NumberNotSeriesN, S<:Real}
-
-    (r == 0) && return one!(c, a, k)
-    (r == 1) && return identity!(c, a, k)
-    (r == 2) && return sqr!(c, a, k)
+    isinteger(r) && r > 0 && return pow!(c, a, aux, Int(r), k)
     (r == 0.5) && return sqrt!(c, a, k)
-
+    # 0-th order coeff
     if k == 0
         @inbounds c[0][1] = ( constant_term(a) )^r
         return nothing
     end
-
     # Sanity
     zero!(c, k)
-
     # The recursion formula
     @inbounds for i = 0:k-1
         aux = r*(k-i) - i
@@ -274,25 +268,31 @@ end
     end
     # c[k] <- c[k]/(k * constant_term(a))
     @inbounds div!(c[k], c[k], k * constant_term(a))
-
     return nothing
 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}
+# 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
+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)
-
     # Sanity
     zero!(res, ordT)
-
     # First non-zero coefficient
     l0 = findfirst(a)
     l0 < 0 && return nothing
-
     # The first non-zero coefficient of the result; must be integer
     !isinteger(r*l0) && throw(DomainError(a,
         """The 0-th order Taylor1 coefficient must be non-zero
@@ -300,28 +300,23 @@ end
     lnull = trunc(Int, r*l0 )
     kprime = ordT-lnull
     (kprime < 0 || lnull > a.order) && return nothing
-
     # Relevant for positive integer r, to avoid round-off errors
     isinteger(r) && r > 0 && (ordT > r*findlast(a)) && return nothing
-
     if ordT == lnull
-        if isinteger(r)
-            power_by_squaring!(res[ordT], a[l0], aux[0], round(Int,r))
+        if isinteger(r) && r > 0
+            pow!(res[ordT], a[l0], aux[0], round(Int, r), 1)
             return nothing
         end
-
         a0 = constant_term(a[l0])
         iszero(a0) && throw(DomainError(a[l0],
             """The 0-th order TaylorN coefficient must be non-zero
             in order to expand `^` around 0."""))
-
+        # Recursion formula
         for ordQ in eachindex(a[l0])
             pow!(res[ordT], a[l0], aux[0], r, ordQ)
         end
-
         return nothing
     end
-
     # The recursion formula
     for i = 0:ordT-lnull-1
         ((i+lnull) > a.order || (l0+kprime-i > a.order)) && continue
@@ -331,25 +326,9 @@ end
     end
     # res[ordT] /= a[l0]*kprime
     @inbounds div_scalar!(res[ordT], 1/kprime, a[l0])
-
     return nothing
 end
 
-for T in (:Taylor1, :TaylorN)
-    @eval begin
-        @inline function pow!(res::$T{T}, a::$T{T}, aux::$T{T},
-                              r::Integer, k::Int) where {T<:NumberNotSeries}
-            (r == 0) && return one!(res, a, k)
-            (r == 1) && return identity!(res, a, k)
-            (r == 2) && return sqr!(res, a, k)
-            # Sanity
-            zero!(res, k)
-            power_by_squaring!(res, a, aux, r)
-            return nothing
-        end
-    end
-end
-
 
 ## Square ##
 """

From f0c823aa205e97599471cd3b5428f8684f3dcb27 Mon Sep 17 00:00:00 2001
From: Luis Benet <benet@icf.unam.mx>
Date: Sat, 24 Aug 2024 17:45:38 -0600
Subject: [PATCH 4/8] Fixes in _pow for Interval coeffs, and tests

---
 ext/TaylorSeriesIAExt.jl | 73 +++++++++++++++-------------------------
 test/intervals.jl        | 27 +++++++++++----
 2 files changed, 48 insertions(+), 52 deletions(-)

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)

From b9dcbf31014d7555d7db194e816b14517c8fd380 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 5/8] 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

From ee69417afc64109b6bbd714fbfb9f9b891831973 Mon Sep 17 00:00:00 2001
From: Luis Benet <benet@icf.unam.mx>
Date: Fri, 4 Oct 2024 17:42:46 -0600
Subject: [PATCH 6/8] Metaprogramming and fixes in pow!

for integer power, pow! uses power_by_squaring, is defined for TaylorN, or when coeffs are intervals
---
 ext/TaylorSeriesIAExt.jl | 13 +++++++-
 src/power.jl             | 65 +++++++++++++++++++---------------------
 2 files changed, 43 insertions(+), 35 deletions(-)

diff --git a/ext/TaylorSeriesIAExt.jl b/ext/TaylorSeriesIAExt.jl
index a6b541f3..2e8ab91a 100644
--- a/ext/TaylorSeriesIAExt.jl
+++ b/ext/TaylorSeriesIAExt.jl
@@ -12,7 +12,7 @@ isdefined(Base, :get_extension) ? (using IntervalArithmetic) : (using ..Interval
 # _pow
 for T in (:Taylor1, :TaylorN)
     @eval begin
-        # Uses TS.power_by_squaring! through `pow!`
+        # Uses TS.power_by_squaring! through TS.pow!
         function TS._pow(a::$T{Interval{S}}, n::Integer) where {S<:Real}
             a0 = constant_term(a)
             aux = one(a0)^n
@@ -58,6 +58,17 @@ function TS._pow(a::TaylorN{Interval{T}}, r::S) where {T<:Real, S<:Real}
     return c
 end
 
+for T in (:Taylor1, :TaylorN)
+    @eval @inline function TS.pow!(res::$T{Interval{T}}, a::$T{Interval{T}}, aux::$T{Interval{T}},
+            r::S, k::Int) where {T<:Real, S<:Integer}
+        (r == 0) && return TS.one!(res, a, k)
+        (r == 1) && return TS.identity!(res, a, k)
+        (r == 2) && return TS.sqr!(res, a, k)
+        TS.power_by_squaring!(res, a, aux, r)
+        return nothing
+    end
+end
+
 
 for T in (:Taylor1, :TaylorN)
     @eval function log(a::$T{Interval{S}}) where {S<:Real}
diff --git a/src/power.jl b/src/power.jl
index 901bc59b..6c31f85f 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -37,27 +37,25 @@ end
 
 
 ## Real power ##
-function ^(a::Taylor1{T}, r::S) where {T<:Number, S<:Real}
-    a0 = constant_term(a)
-    aux = one(a0)^r
-    iszero(r) && return Taylor1(aux, a.order)
-    aa = aux*a
-    r == 1 && return aa
-    r == 2 && return square(aa)
-    r == 0.5 && return sqrt(aa)
-    return _pow(aa, r)
-end
-
-function ^(a::TaylorN{T}, r::S) where {T<:Number, S<:Real}
-    a0 = constant_term(a)
-    aux = one(a0^r)
-    aa = aux*a
-    isinteger(r) && r ≥ 0 && return _pow(aa, round(Int, r))
-    iszero(a0) && throw(DomainError(a,
-        """The 0-th order TaylorN coefficient must be non-zero
-        in order to expand `^` around 0."""))
-    r == 0.5 && return sqrt(aa)
-    return _pow(aa, r)
+for T in (:Taylor1, :TaylorN)
+    @eval function ^(a::$T{T}, r::S) where {T<:Number, S<:Real}
+        a0 = constant_term(a)
+        aux = one(a0)^r
+        iszero(r) && return $T(aux, a.order)
+        aa = aux*a
+        r == 1 && return aa
+        r == 2 && return square(aa)
+        r == 0.5 && return sqrt(aa)
+        if $T == TaylorN
+            if iszero(aa[0])
+                isinteger(r) && r ≥ 0 && return _pow(aa, round(Int, r))
+                throw(DomainError(a,
+                """The 0-th order TaylorN coefficient must be non-zero
+                in order to expand `^` around 0."""))
+            end
+        end
+        return _pow(aa, r)
+    end
 end
 
 
@@ -79,7 +77,7 @@ for T in (:Taylor1, :TaylorN)
 end
 
 function _pow(a::Taylor1{T}, r::S) where {T<:Number, S<:Real}
-    aux = one(constant_term(a)^r)
+    aux = one(constant_term(a))^r
     l0 = findfirst(a)
     lnull = trunc(Int, r*l0 )
     (lnull > a.order) && return Taylor1( zero(aux), a.order)
@@ -93,7 +91,7 @@ function _pow(a::Taylor1{T}, r::S) where {T<:Number, S<:Real}
 end
 
 function _pow(a::TaylorN{T}, r::S) where {T<:Number, S<:Real}
-    aux = one(constant_term(a)^r)
+    aux = one(constant_term(a))^r
     c = TaylorN( zero(aux), a.order)
     aux0 = zero(c)
     for ord in eachindex(a)
@@ -151,6 +149,7 @@ end
 # Licensed under MIT "Expat"
 for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
     @eval function power_by_squaring(x::$T, p::Integer)
+        @assert p ≥ 0
         (p == 0) && return one(x)
         (p == 1) && return copy(x)
         (p == 2) && return square(x)
@@ -177,7 +176,7 @@ end
 # uses internally mutating method `power_by_squaring!`
 for T in (:Taylor1, :TaylorN)
     @eval function power_by_squaring(x::$T{T}, p::Integer) where {T<:NumberNotSeries}
-        @assert p > 0
+        @assert p ≥ 0
         (p == 0) && return one(x)
         (p == 1) && return copy(x)
         (p == 2) && return square(x)
@@ -249,7 +248,7 @@ end
 
 @inline function pow!(c::TaylorN{T}, a::TaylorN{T}, aux::TaylorN{T},
                       r::S, k::Int) where {T<:NumberNotSeriesN, S<:Real}
-    isinteger(r) && r > 0 && return pow!(c, a, aux, Int(r), k)
+    # isinteger(r) && r > 0 && return pow!(c, a, aux, Int(r), k)
     (r == 0.5) && return sqrt!(c, a, k)
     # 0-th order coeff
     if k == 0
@@ -270,15 +269,13 @@ end
 end
 
 # Uses power_by_squaring!
-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
+@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
 end
 
 @inline function pow!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},

From 774e83199744b0e2c8791ef7cfb98ad1ca866a08 Mon Sep 17 00:00:00 2001
From: Luis Benet <benet@icf.unam.mx>
Date: Wed, 9 Oct 2024 13:22:48 -0600
Subject: [PATCH 7/8] Further minor fixes

---
 ext/TaylorSeriesIAExt.jl | 114 +++++++++++++++++++++------------------
 src/power.jl             |  28 +++++-----
 2 files changed, 76 insertions(+), 66 deletions(-)

diff --git a/ext/TaylorSeriesIAExt.jl b/ext/TaylorSeriesIAExt.jl
index 2e8ab91a..f9dd1227 100644
--- a/ext/TaylorSeriesIAExt.jl
+++ b/ext/TaylorSeriesIAExt.jl
@@ -8,64 +8,76 @@ import TaylorSeries: evaluate, _evaluate, normalize_taylor, square
 
 isdefined(Base, :get_extension) ? (using IntervalArithmetic) : (using ..IntervalArithmetic)
 
-
-# _pow
 for T in (:Taylor1, :TaylorN)
     @eval begin
-        # Uses TS.power_by_squaring! through TS.pow!
-        function TS._pow(a::$T{Interval{S}}, n::Integer) where {S<:Real}
-            a0 = constant_term(a)
-            aux = one(a0)^n
+        @eval ^(a::$T{Interval{T}}, n::Integer) where {T<:Real} = TS.power_by_squaring(a, n)
+
+        @eval ^(a::$T{Interval{T}}, r::Rational) where {T<:Real} = a^float(r)
+
+        @eval function ^(a::$T{Interval{T}}, r::S) where {T<:Real, S<:Real}
+            isinteger(r) && r ≥ 0 && return TS.power_by_squaring(a, Integer(r))
+            a0 = constant_term(a) ∩ Interval(zero(T), T(Inf))
+            @assert !isempty(a0)
+            aux = one(a0^r)
             a[0] = aux * a0
-            c = $T( zero(aux), a.order)
-            caux = zero(c)
-            # TS.power_by_squaring!(c, a, caux, n)
-            TS.pow!(c, a, caux, n, 1)
-            return c
+            r == 0.5 && return sqrt(a)
+            if $T == TaylorN
+                if iszero(a0)
+                    throw(DomainError(a,
+                    """The 0-th order TaylorN coefficient must be non-zero
+                    in order to expand `^` around 0."""))
+                end
+            end
+            return TS._pow(a, r)
         end
-    end
-end
 
-function TS._pow(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real}
-    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
-    #
-    l0 = findfirst(a)
-    lnull = trunc(Int, r*l0 )
-    (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 = zero(c)
-    for k in eachindex(c)
-        TS.pow!(c, a, aux0, r, k)
-    end
-    return c
-end
+        # _pow
+        function TS._pow(a::$T{Interval{S}}, n::Integer) where {S<:Real}
+            n < 0 && return TS._pow(a, float(n))
+            return TS.power_by_squaring(a, n)
+        end
 
-function TS._pow(a::TaylorN{Interval{T}}, r::S) where {T<:Real, S<:Real}
-    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
-    #
-    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
+        function TS._pow(a::$T{Interval{T}}, r::S) where {T<:Real, S<:Real}
+            isinteger(r) && r ≥ 0 && return TS.power_by_squaring(a, Integer(r))
+            a0 = constant_term(a) ∩ Interval(zero(T), T(Inf))
+            @assert !isempty(a0)
+            aux = one(a0^r)
+            a[0] = aux * a0
+            r == 0.5 && return sqrt(a)
+            if $T == Taylor1
+                l0 = findfirst(a)
+                # Index of first non-zero coefficient of the result; must be integer
+                !isinteger(r*l0) && throw(DomainError(a,
+                    """The 0-th order Taylor1 coefficient must be non-zero
+                    to raise the Taylor1 polynomial to a non-integer exponent."""))
+                lnull = trunc(Int, r*l0 )
+                (lnull > a.order) && return $T( zero(aux), a.order)
+                c_order = l0 == 0 ? a.order : min(a.order, trunc(Int, r*a.order))
+            else
+                if iszero(a0)
+                    throw(DomainError(a,
+                    """The 0-th order TaylorN coefficient must be non-zero
+                    in order to expand `^` around 0."""))
+                end
+                c_order = a.order
+            end
+            #
+            c = $T(zero(aux), c_order)
+            aux0 = zero(c)
+            for k in eachindex(c)
+                TS.pow!(c, a, aux0, r, k)
+            end
+            return c
+        end
 
-for T in (:Taylor1, :TaylorN)
-    @eval @inline function TS.pow!(res::$T{Interval{T}}, a::$T{Interval{T}}, aux::$T{Interval{T}},
-            r::S, k::Int) where {T<:Real, S<:Integer}
-        (r == 0) && return TS.one!(res, a, k)
-        (r == 1) && return TS.identity!(res, a, k)
-        (r == 2) && return TS.sqr!(res, a, k)
-        TS.power_by_squaring!(res, a, aux, r)
-        return nothing
+        function TS.pow!(res::$T{Interval{T}}, a::$T{Interval{T}}, aux::$T{Interval{T}},
+                r::S, k::Int) where {T<:Real, S<:Integer}
+            (r == 0) && return TS.one!(res, a, k)
+            (r == 1) && return TS.identity!(res, a, k)
+            (r == 2) && return TS.sqr!(res, a, k)
+            TS.power_by_squaring!(res, a, aux, r)
+            return nothing
+        end
     end
 end
 
diff --git a/src/power.jl b/src/power.jl
index 6c31f85f..7e90bb85 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -31,24 +31,22 @@ 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}}, r::Rational) where {T<:NumberNotSeries} = a^float(r)
-
 
 ## Real power ##
 for T in (:Taylor1, :TaylorN)
     @eval function ^(a::$T{T}, r::S) where {T<:Number, S<:Real}
         a0 = constant_term(a)
-        aux = one(a0)^r
+        aux = a0^zero(r)
         iszero(r) && return $T(aux, a.order)
         aa = aux*a
         r == 1 && return aa
         r == 2 && return square(aa)
+        if $T == TaylorN
+            isinteger(r) && r ≥ 0 && return TS.power_by_squaring(a, Integer(r))
+        end
         r == 0.5 && return sqrt(aa)
         if $T == TaylorN
-            if iszero(aa[0])
-                isinteger(r) && r ≥ 0 && return _pow(aa, round(Int, r))
+            if iszero(a0)
                 throw(DomainError(a,
                 """The 0-th order TaylorN coefficient must be non-zero
                 in order to expand `^` around 0."""))
@@ -65,10 +63,7 @@ _pow(a::Taylor1, n::Integer) = a^float(n)
 _pow(a::TaylorN, n::Integer) = power_by_squaring(a, n)
 
 for T in (:Taylor1, :TaylorN)
-    @eval function _pow(a::$T{T}, n::Integer) where {T<:Integer}
-        n < 0 && throw(DomainError())
-        return power_by_squaring(a, n)
-    end
+    @eval _pow(a::$T{T}, n::Integer) where {T<:Integer} = power_by_squaring(a, n)
 
     @eval function _pow(a::$T{Rational{T}}, n::Integer) where {T<:Integer}
         n < 0 && return inv( a^(-n) )
@@ -78,6 +73,7 @@ end
 
 function _pow(a::Taylor1{T}, r::S) where {T<:Number, S<:Real}
     aux = one(constant_term(a))^r
+    iszero(r) && return Taylor1(aux, a.order)
     l0 = findfirst(a)
     lnull = trunc(Int, r*l0 )
     (lnull > a.order) && return Taylor1( zero(aux), a.order)
@@ -91,6 +87,7 @@ function _pow(a::Taylor1{T}, r::S) where {T<:Number, S<:Real}
 end
 
 function _pow(a::TaylorN{T}, r::S) where {T<:Number, S<:Real}
+    isinteger(r) && r ≥ 0 && return power_by_squaring(a, Integer(r))
     aux = one(constant_term(a))^r
     c = TaylorN( zero(aux), a.order)
     aux0 = zero(c)
@@ -219,7 +216,7 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`.
     # First non-zero coefficient
     l0 = findfirst(a)
     l0 < 0 && return nothing
-    # The first non-zero coefficient of the result; must be integer
+    # Index of first non-zero coefficient of the result; must be integer
     !isinteger(r*l0) && throw(DomainError(a,
         """The 0-th order Taylor1 coefficient must be non-zero
         to raise the Taylor1 polynomial to a non-integer exponent."""))
@@ -248,7 +245,7 @@ end
 
 @inline function pow!(c::TaylorN{T}, a::TaylorN{T}, aux::TaylorN{T},
                       r::S, k::Int) where {T<:NumberNotSeriesN, S<:Real}
-    # isinteger(r) && r > 0 && return pow!(c, a, aux, Int(r), k)
+    isinteger(r) && r > 0 && return pow!(c, a, aux, Integer(r), k)
     (r == 0.5) && return sqrt!(c, a, k)
     # 0-th order coeff
     if k == 0
@@ -300,11 +297,12 @@ end
     # Relevant for positive integer r, to avoid round-off errors
     isinteger(r) && r > 0 && (ordT > r*findlast(a)) && return nothing
     if ordT == lnull
+        a0 = constant_term(a[l0])
         if isinteger(r) && r > 0
-            pow!(res[ordT], a[l0], aux[0], round(Int, r), 1)
+            pow!(res[ordT], a[l0], aux[0], round(Integer, r), 1)
+            # power_by_squaring!(res[ordT], a[l0], aux[0], round(Integer, r))
             return nothing
         end
-        a0 = constant_term(a[l0])
         iszero(a0) && throw(DomainError(a[l0],
             """The 0-th order TaylorN coefficient must be non-zero
             in order to expand `^` around 0."""))

From 81ff6204d4482a50a94f5865c3b096b51ad2c2ea Mon Sep 17 00:00:00 2001
From: Luis Benet <benet@icf.unam.mx>
Date: Wed, 9 Oct 2024 15:22:50 -0600
Subject: [PATCH 8/8] Another minor fix

---
 src/power.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/power.jl b/src/power.jl
index 7e90bb85..d5f38b57 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -299,8 +299,8 @@ end
     if ordT == lnull
         a0 = constant_term(a[l0])
         if isinteger(r) && r > 0
-            pow!(res[ordT], a[l0], aux[0], round(Integer, r), 1)
-            # power_by_squaring!(res[ordT], a[l0], aux[0], round(Integer, r))
+            # pow!(res[ordT], a[l0], aux[0], round(Integer, r), 1)
+            power_by_squaring!(res[ordT], a[l0], aux[0], round(Integer, r))
             return nothing
         end
         iszero(a0) && throw(DomainError(a[l0],