From 8951e74c6ffc4f18b8898aec7b31a533cc7927dc Mon Sep 17 00:00:00 2001
From: Jorge Perez-Hernandez <jorge.perez_hernandez@cfa.harvard.edu>
Date: Fri, 29 May 2026 17:44:38 -0400
Subject: [PATCH 01/20] Some performance fixes for Taylor1; add nowarn kwarg
 for variables! and use it in tests

---
 src/arithmetic.jl        | 105 +++++++++++++++++++++++++++++++++----
 src/hash_tables.jl       |   7 ++-
 src/identity.jl          |   5 +-
 src/parameters.jl        |  56 ++++++++++++--------
 src/power.jl             | 110 ++++++++++++++++++++++++++-------------
 test/broadcasting.jl     |   6 +--
 test/fateman40.jl        |   2 +-
 test/identities_Euler.jl |   3 +-
 test/intervals.jl        |   4 +-
 test/jld2.jl             |   2 +-
 test/manyvariables.jl    |  42 +++++++--------
 test/mixtures.jl         |  10 ++--
 test/mutatingfuncts.jl   |  10 ++++
 test/runtests.jl         |   2 +-
 test/staticarrays.jl     |   2 +-
 15 files changed, 259 insertions(+), 107 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 2798d496..aa093827 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -317,6 +317,46 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
     end
 end
 
+function add!(v::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}) where
+        {T<:NumberNotSeries}
+    v_coeffs = v.coeffs
+    a_coeffs = a.coeffs
+    b_coeffs = b.coeffs
+    @inbounds for i in eachindex(v_coeffs)
+        v_coeffs[i] = a_coeffs[i] + b_coeffs[i]
+    end
+    return nothing
+end
+
+function subst!(v::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}) where
+        {T<:NumberNotSeries}
+    v_coeffs = v.coeffs
+    a_coeffs = a.coeffs
+    b_coeffs = b.coeffs
+    @inbounds for i in eachindex(v_coeffs)
+        v_coeffs[i] = a_coeffs[i] - b_coeffs[i]
+    end
+    return nothing
+end
+
+function +(a::Taylor1{T}, b::Taylor1{T}) where {T<:NumberNotSeries}
+    if order(a) != order(b)
+        a, b = fixorder(a, b)
+    end
+    c = zero(a)
+    add!(c, a, b)
+    return c
+end
+
+function -(a::Taylor1{T}, b::Taylor1{T}) where {T<:NumberNotSeries}
+    if order(a) != order(b)
+        a, b = fixorder(a, b)
+    end
+    c = zero(a)
+    subst!(c, a, b)
+    return c
+end
+
 for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
     @eval begin
         function ($f)(a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}}) where
@@ -527,13 +567,20 @@ end
 # Internal multiplication functions
 function mul!(c::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}, k::Int) where
         {T<:NumberNotSeries}
-    zero!(c, k)
-    muladd!(c, a, b, k)
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    b_coeffs = b.coeffs
+    kk = k+1
+    @inbounds acc = zero(c_coeffs[kk])
+    @inbounds for i = 1:kk
+        acc += a_coeffs[i] * b_coeffs[kk-i+1]
+    end
+    @inbounds c_coeffs[kk] = acc
     return nothing
 end
 function mul!(v::Taylor1{T}, a::Taylor1{S}, b::NumberNotSeries, k::Int) where
         {T<:NumberNotSeries, S<:NumberNotSeries}
-    @inbounds v[k] = a[k] * b
+    @inbounds v.coeffs[k+1] = a.coeffs[k+1] * b
     return nothing
 end
 mul!(v::Taylor1{T}, a::NumberNotSeries, b::Taylor1{S}, k::Int) where
@@ -541,16 +588,29 @@ mul!(v::Taylor1{T}, a::NumberNotSeries, b::Taylor1{S}, k::Int) where
 #
 function muladd!(c::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}, k::Int) where
         {T<:NumberNotSeries}
-    @inbounds for i = 0:k
-        c[k] += a[i] * b[k-i]
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    b_coeffs = b.coeffs
+    kk = k+1
+    @inbounds acc = c_coeffs[kk]
+    @inbounds for i = 1:kk
+        acc += a_coeffs[i] * b_coeffs[kk-i+1]
     end
+    @inbounds c_coeffs[kk] = acc
     return nothing
 end
 
 function muladd!(c::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}) where
         {T<:NumberNotSeries}
-    for k in eachindex(c)
-        muladd!(c, a, b, k)
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    b_coeffs = b.coeffs
+    @inbounds for kk in eachindex(c_coeffs)
+        acc = c_coeffs[kk]
+        for i = 1:kk
+            acc += a_coeffs[i] * b_coeffs[kk-i+1]
+        end
+        c_coeffs[kk] = acc
     end
     return nothing
 end
@@ -571,8 +631,15 @@ end
 # Implements c[k] = scalar \sum_i a[i] b[k-i]
 function mul_scalar!(c::Taylor1{T}, scalar::NumberNotSeries, a::Taylor1{T},
         b::Taylor1{T}, k::Int) where {T<:NumberNotSeries}
-    mul!(c, a, b, k)
-    mul!(c, scalar, c, k)
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    b_coeffs = b.coeffs
+    kk = k+1
+    @inbounds acc = zero(c_coeffs[kk])
+    @inbounds for i = 1:kk
+        acc += a_coeffs[i] * b_coeffs[kk-i+1]
+    end
+    @inbounds c_coeffs[kk] = scalar * acc
     return nothing
 end
 
@@ -771,8 +838,24 @@ mul!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, b::NumberNotSeries,
     k::Int) where {T<:NumberNotSeries} = mul!(res, b, a, k)
 
 
-# in-place product (assumes equal order among TaylorNs)
-# NOTE: the result of the product is *accumulated* in c[k]
+# in-place product (assumes equal order)
+function mul!(c::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}) where
+        {T<:NumberNotSeries}
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    b_coeffs = b.coeffs
+    @inbounds for kk in eachindex(c_coeffs)
+        acc = zero(c_coeffs[kk])
+        for i = 1:kk
+            acc += a_coeffs[i] * b_coeffs[kk-i+1]
+        end
+        c_coeffs[kk] = acc
+    end
+    return nothing
+end
+
+# Fallback for nested coefficient types; scalar Taylor1 has a coefficient-vector
+# specialization above.
 function mul!(c::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}) where {T<:Number}
     for k in eachindex(c)
         mul!(c, a, b, k)
diff --git a/src/hash_tables.jl b/src/hash_tables.jl
index cb5480f0..29976339 100644
--- a/src/hash_tables.jl
+++ b/src/hash_tables.jl
@@ -242,8 +242,11 @@ mutable. Module loading initializes `default_space[]`; later compatibility calls
 such as `variables!` replace `default_space[]` with `space`. Existing `TaylorN`
 and `HomogeneousPolynomial` objects keep their original spaces, while future
 default-space constructors use the new default algebra.
+
+Use `nowarn=true` to suppress the warning emitted when replacing the default
+space.
 """
-function set_default_space!(space::JetSpace)
+function set_default_space!(space::JetSpace; nowarn::Bool=false)
     old_space = default_space[]
     old_space === space && return space
 
@@ -254,7 +257,7 @@ function set_default_space!(space::JetSpace)
     old_numvars = get_numvars(old_space)
     new_order = order(space)
     new_numvars = get_numvars(space)
-    @warn msg old_order old_numvars new_order new_numvars
+    nowarn || @warn msg old_order old_numvars new_order new_numvars
 
     # Replace the active default space. Do not mutate the old object in place:
     # existing TaylorN/HomogeneousPolynomial objects may still depend on it.
diff --git a/src/identity.jl b/src/identity.jl
index 5469cf0b..7197ce20 100644
--- a/src/identity.jl
+++ b/src/identity.jl
@@ -17,8 +17,11 @@ for T in (:Taylor1, :TaylorN)
         identity!(c[k], a[k])
 end
 
-identity!(c::Taylor1{T}, a::Taylor1{T}, k::Int) where {T<:NumberNotSeries} =
+function identity!(c::Taylor1{T}, a::Taylor1{T}, k::Int) where
+        {T<:NumberNotSeries}
     c[k] = a[k]
+    return nothing
+end
 
 @inline function identity!(c::Taylor1{T}, a::T, k::Int) where {T<:NumberNotSeries}
     zero!(c[k])
diff --git a/src/parameters.jl b/src/parameters.jl
index 5a5c107e..11dc9b69 100644
--- a/src/parameters.jl
+++ b/src/parameters.jl
@@ -230,14 +230,15 @@ variables(space::JetSpace; order::Int=TS.order(space)) =
     variables(Float64, space, order=order)
 
 """
-    variables!([T::Type], names::String; [order=order(), numvars=-1])
+    variables!([T::Type], names::String; [order=order(), numvars=-1, nowarn=false])
 
 Return a `TaylorN{T}` vector with each entry representing an
 independent variable. `names` defines the output for each variable
 (separated by a space). The default type `T` is `Float64`,
 and the default for `order` is the current default-space order.
 
-This function updates the compatibility default `JetSpace`.
+This function updates the compatibility default `JetSpace`. Set `nowarn=true`
+to suppress the warning emitted when the default space changes.
 
 If `numvars` is not specified, it is inferred from `names`. If only
 one variable name is defined and `numvars>1`, it uses this name with
@@ -261,34 +262,47 @@ julia> variables!("x", order=6, numvars=2)
   1.0 x₂ + 𝒪(‖x‖⁷)
 ```
 """
-function variables!(::Type{R}, names::Vector{T}; order=TS.order()) where
+function variables!(::Type{R}, names::Vector{T}; order=TS.order(),
+        nowarn::Bool=false) where
         {R, T<:AbstractString}
 
-    space = set_default_space!(JetSpace(order, _variable_names_from(names)))
+    space = set_default_space!(JetSpace(order, _variable_names_from(names));
+        nowarn=nowarn)
 
     # return a list of the new variables
     return variables(R, space)
 end
-variables!(::Type{R}, symbs::Vector{T}; order=TS.order()) where
-    {R,T<:Symbol} = variables!(R, string.(symbs), order=order)
-
-variables!(names::Vector{T}; order=TS.order()) where {T<:AbstractString} =
-    variables!(Float64, names, order=order)
-variables!(symbs::Vector{T}; order=TS.order()) where {T<:Symbol} =
-    variables!(Float64, symbs, order=order)
-
-function variables!(::Type{R}, names::T; order=TS.order(), numvars=-1) where
+variables!(::Type{R}, symbs::Vector{T}; order=TS.order(),
+    nowarn::Bool=false) where {R,T<:Symbol} =
+        variables!(R, string.(symbs), order=order, nowarn=nowarn)
+
+variables!(names::Vector{T}; order=TS.order(),
+    nowarn::Bool=false) where {T<:AbstractString} =
+        variables!(Float64, names, order=order, nowarn=nowarn)
+variables!(symbs::Vector{T}; order=TS.order(),
+    nowarn::Bool=false) where {T<:Symbol} =
+        variables!(Float64, symbs, order=order, nowarn=nowarn)
+
+function variables!(::Type{R}, names::T; order=TS.order(), numvars=-1,
+        nowarn::Bool=false) where
         {R,T<:AbstractString}
 
-    return variables!(R, _variable_names_from(names, numvars=numvars), order=order)
+    return variables!(R, _variable_names_from(names, numvars=numvars),
+        order=order, nowarn=nowarn)
 end
-variables!(::Type{R}, symbs::Symbol; order=TS.order(), numvars=-1) where {R} =
-    variables!(R, string(symbs), order=order, numvars=numvars)
-
-variables!(names::T; order=TS.order(), numvars=-1) where {T<:AbstractString} =
-    variables!(Float64, names, order=order, numvars=numvars)
-variables!(symbs::Symbol; order=TS.order(), numvars=-1) =
-    variables!(Float64, string(symbs), order=order, numvars=numvars)
+variables!(::Type{R}, symbs::Symbol; order=TS.order(), numvars=-1,
+    nowarn::Bool=false) where {R} =
+        variables!(R, string(symbs), order=order, numvars=numvars,
+            nowarn=nowarn)
+
+variables!(names::T; order=TS.order(), numvars=-1,
+    nowarn::Bool=false) where {T<:AbstractString} =
+        variables!(Float64, names, order=order, numvars=numvars,
+            nowarn=nowarn)
+variables!(symbs::Symbol; order=TS.order(), numvars=-1,
+    nowarn::Bool=false) =
+        variables!(Float64, string(symbs), order=order, numvars=numvars,
+            nowarn=nowarn)
 
 
 """
diff --git a/src/power.jl b/src/power.jl
index f73c23de..d0934dac 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -69,6 +69,7 @@ end
 
 # _pow
 _pow(a::Taylor1, n::Integer) = a^float(n)
+_pow(a::Taylor1{T}, n::Integer) where {T<:NumberNotSeries} = a^float(n)
 
 _pow(a::TaylorN, n::Integer) = power_by_squaring(a, n)
 
@@ -81,6 +82,22 @@ for T in (:Taylor1, :TaylorN)
     end
 end
 
+function _pow(a::Taylor1{T}, r::S) where {T<:NumberNotSeries, S<:Real}
+    aux = one(constant_term(a))^r
+    iszero(r) && return Taylor1(aux, order(a))
+    l0 = findfirst(a)
+    lnull = trunc(Int, r*l0 )
+    (lnull > order(a)) && return Taylor1( zero(aux), order(a))
+    c_order = l0 == 0 ? order(a) : min(order(a), trunc(Int, r*order(a)))
+    c = Taylor1(zero(aux), c_order)
+    order_a = order(a)
+    lastnz = isinteger(r) && r > 0 ? findlast(a) : order_a
+    for k in eachindex(c)
+        _pow_taylor1_cached!(c, a, r, k, l0, lnull, lastnz, order_a)
+    end
+    return c
+end
+
 function _pow(a::Taylor1{T}, r::S) where {T<:Number, S<:Real}
     aux = one(constant_term(a))^r
     iszero(r) && return Taylor1(aux, order(a))
@@ -214,44 +231,57 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`.
 
 """ pow!
 
-function pow!(c::Taylor1{T}, a::Taylor1{T}, aux::Taylor1{T},
-                r::S, k::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, constant_term(aux), k)
-    (r == 0.5) && return sqrt!(c, a, aux, k)
-    # Sanity
+@inline function _pow_taylor1_cached!(c::Taylor1{T}, a::Taylor1{T},
+        r::S, k::Int, l0::Int, lnull::Int, lastnz::Int,
+        order_a::Int) where {T<:NumberNotSeries, S<:Real}
     zero!(c, k)
-    # First non-zero coefficient
-    l0 = findfirst(a)
     l0 < 0 && return nothing
-    # 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 )
     kprime = k-lnull
-    (kprime < 0 || lnull > order(a)) && 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
+    (kprime < 0 || lnull > order_a) && return nothing
+    isinteger(r) && r > 0 && (k > r*lastnz) && return nothing
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    kk = k+1
+    @inbounds a_l0 = a_coeffs[l0+1]
     if k == lnull
-        @inbounds c[k] = (a[l0])^float(r)
+        @inbounds c_coeffs[kk] = a_l0^float(r)
         return nothing
     end
-    # The recursion formula
-    if l0+kprime ≤ order(a)
-        @inbounds c[k] = r * kprime * c[lnull] * a[l0+kprime]
+
+    @inbounds acc = zero(c_coeffs[kk])
+    if l0+kprime ≤ order_a
+        @inbounds acc = r * kprime * c_coeffs[lnull+1] * a_coeffs[l0+kprime+1]
     end
-    for i = 1:k-lnull-1
-        ((i+lnull) > order(a) || (l0+kprime-i > order(a))) && continue
+    ilo = max(1, l0+kprime-order_a)
+    ihi = min(k-lnull-1, order_a-lnull)
+    @inbounds for i = ilo:ihi
         aaux = r*(kprime-i) - i
-        @inbounds c[k] += aaux * c[i+lnull] * a[l0+kprime-i]
+        acc += aaux * c_coeffs[i+lnull+1] * a_coeffs[l0+kprime-i+1]
     end
-    @inbounds c[k] = c[k] / (kprime * a[l0])
+    @inbounds c_coeffs[kk] = acc / (kprime * a_l0)
     return nothing
 end
 
+function pow!(c::Taylor1{T}, a::Taylor1{T}, aux::Taylor1{T},
+                r::S, k::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, constant_term(aux), k)
+    (r == 0.5) && return sqrt!(c, a, aux, k)
+    l0 = findfirst(a)
+    if l0 < 0
+        zero!(c, k)
+        return nothing
+    end
+    lnull = trunc(Int, r*l0)
+    order_a = order(a)
+    lastnz = isinteger(r) && r > 0 ? findlast(a) : order_a
+    return _pow_taylor1_cached!(c, a, r, k, l0, lnull, lastnz, order_a)
+end
+
 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, Integer(r), k)
@@ -487,17 +517,21 @@ function sqr!(c::Taylor1{T}, a::Taylor1{T}, ::T, k::Int) where {T<:Number}
         sqr_orderzero!(c, a)
         return nothing
     end
-    # Sanity
-    zero!(c, k)
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    kk = k+1
+    @inbounds acc = zero(c_coeffs[kk])
     # Recursion formula
     kodd = k%2
     kend = (k - 2 + kodd) >> 1
-    @inbounds for i = 0:kend
-        c[k] += a[i] * a[k-i]
+    @inbounds for i = 1:kend+1
+        acc += a_coeffs[i] * a_coeffs[kk-i+1]
     end
-    @inbounds c[k] = 2 * c[k]
-    kodd == 1 && return nothing
-    @inbounds c[k] += a[k >> 1]^2
+    acc = 2 * acc
+    if kodd == 0
+        @inbounds acc += a_coeffs[(k >> 1)+1]^2
+    end
+    @inbounds c_coeffs[kk] = acc
     return nothing
 end
 
@@ -526,15 +560,19 @@ function sqr!(c::Taylor1{T}, k::Int) where {T<:NumberNotSeries}
         sqr_orderzero!(c, c)
         return nothing
     end
+    c_coeffs = c.coeffs
+    kk = k+1
     # Recursion formula
     kodd = k%2
     kend = (k - 2 + kodd) >> 1
-    (kend >= 0) && ( @inbounds c[k] = c[0] * c[k] )
-    @inbounds for i = 1:kend
-        c[k] += c[i] * c[k-i]
-    end
-    @inbounds c[k] = 2 * c[k]
-    (kodd == 0) && ( @inbounds c[k] += c[k >> 1]^2 )
+    @inbounds acc = zero(c_coeffs[kk])
+    (kend >= 0) && ( @inbounds acc = c_coeffs[1] * c_coeffs[kk] )
+    @inbounds for i = 2:kend+1
+        acc += c_coeffs[i] * c_coeffs[kk-i+1]
+    end
+    acc = 2 * acc
+    (kodd == 0) && ( @inbounds acc += c_coeffs[(k >> 1)+1]^2 )
+    @inbounds c_coeffs[kk] = acc
     return nothing
 end
 
diff --git a/test/broadcasting.jl b/test/broadcasting.jl
index 4d9a67c5..2b582069 100644
--- a/test/broadcasting.jl
+++ b/test/broadcasting.jl
@@ -60,7 +60,7 @@ using Test
 end
 
 @testset "Broadcasting with HomogeneousPolynomial and TaylorN" begin
-    x, y = variables!("x y", order=3)
+    x, y = variables!("x y", order=3, nowarn=true)
     xH = x[1]
     yH = y[1]
 
@@ -106,7 +106,7 @@ end
 end
 
 @testset "Broadcasting with mixtures Taylor1{TaylorN{T}}" begin
-    x, y = variables!("x", numvars=2, order=6)
+    x, y = variables!("x", numvars=2, order=6, nowarn=true)
     tN = Taylor1(TaylorN{Float64}, 3)
 
     @test tN .== tN
@@ -133,7 +133,7 @@ end
 end
 
 @testset "Broadcasting with mixtures TaylorN{Taylor1{T}}" begin
-    variables!("x", numvars=2, order=6)
+    variables!("x", numvars=2, order=6, nowarn=true)
     t = Taylor1(3)
     xHt = HomogeneousPolynomial([one(t), zero(t)])
     yHt = HomogeneousPolynomial([zero(t), t])
diff --git a/test/fateman40.jl b/test/fateman40.jl
index 3d9c243f..23a7eb3c 100644
--- a/test/fateman40.jl
+++ b/test/fateman40.jl
@@ -5,7 +5,7 @@ using TaylorSeries
 using Test
 
 @testset "Test inspired by Fateman (takes a few seconds)" begin
-    x, y, z, w = variables!(Int128, "x", numvars=4, order=40)
+    x, y, z, w = variables!(Int128, "x", numvars=4, order=40, nowarn=true)
 
     function fateman2(degree::Int)
         T = Int128
diff --git a/test/identities_Euler.jl b/test/identities_Euler.jl
index 0536a4ed..2bdb44ee 100644
--- a/test/identities_Euler.jl
+++ b/test/identities_Euler.jl
@@ -10,7 +10,8 @@ using Test
     variable_names = String[make_variable("α", i) for i in 1:4]
     append!(variable_names, [make_variable("β", i) for i in 1:4])
 
-    a1, a2, a3, a4, b1, b2, b3, b4 = variables!(variable_names, order=4)
+    a1, a2, a3, a4, b1, b2, b3, b4 =
+        variables!(variable_names, order=4, nowarn=true)
 
     lhs1 = a1^2 + a2^2 + a3^2 + a4^2
     lhs2 = b1^2 + b2^2 + b3^2 + b4^2
diff --git a/test/intervals.jl b/test/intervals.jl
index 8c8e2c6b..ee40dbf6 100644
--- a/test/intervals.jl
+++ b/test/intervals.jl
@@ -22,7 +22,7 @@ setdisplay(:full)
     p5(x, a) = x^5 + ifive*a*x^4 + iten*a^2*x^3 + iten*a^3*x^2 + ifive*a^4*x + a^5
 
     ti = Taylor1(Interval{Float64}, 10)
-    x, y = variables!(Interval{Float64}, "x y", order=6)
+    x, y = variables!(Interval{Float64}, "x y", order=6, nowarn=true)
 
     @test eltype(ti) == Taylor1{Interval{Float64}}
     @test eltype(x) == TaylorN{Interval{Float64}}
@@ -217,7 +217,7 @@ setdisplay(:full)
     # Iss 405 (TM iss 158)
     for TT in (:Float32, :Float64, :BigFloat)
         @eval begin
-            x₁, x₂, x₃ = variables!($TT, ["x₁", "x₂", "x₃"]; order=5)
+            x₁, x₂, x₃ = variables!($TT, ["x₁", "x₂", "x₃"]; order=5, nowarn=true)
             Dx₁ = interval(($TT)(1.0), ($TT)(3.0))
             Dx₂ = interval(($TT)(-1.0), ($TT)(1.0))
             Dx₃ = interval(($TT)(-1.0), ($TT)(0.0))
diff --git a/test/jld2.jl b/test/jld2.jl
index fc8294ae..8bf11126 100644
--- a/test/jld2.jl
+++ b/test/jld2.jl
@@ -13,7 +13,7 @@ function same_taylorn_coefficients(a::TaylorN, b::TaylorN)
 end
 
 @testset "Test TaylorSeries JLD2 extension" begin
-    dq = variables!("q", order=4, numvars=6)
+    dq = variables!("q", order=4, numvars=6, nowarn=true)
     random_TaylorN = [exp(sum(rand(6) .* dq)) for _ in 1:10_000]
     jldsave("test.jld2"; random_TaylorN = random_TaylorN)
     recovered_taylorN = JLD2.load("test.jld2", "random_TaylorN")
diff --git a/test/manyvariables.jl b/test/manyvariables.jl
index 4a3ffe41..b9d16085 100644
--- a/test/manyvariables.jl
+++ b/test/manyvariables.jl
@@ -13,7 +13,7 @@ using Test
     @test b isa TaylorN{Complex{Float64}}
     @test_throws MethodError TaylorN(-4//7im)
     # Issue #85 is solved!
-    variables!("x", numvars=66, order=1)
+    variables!("x", numvars=66, order=1, nowarn=true)
     sp = TS.default_space[]
     @test order(sp) == order() == 1
     @test get_numvars(sp) == get_numvars() == 66
@@ -24,7 +24,7 @@ using Test
     @test sp.index_table[2][1] == 3^65
     @test sp.pos_table[2][3^64] == 2
     #
-    variables!("x", numvars=66, order=2)
+    variables!("x", numvars=66, order=2, nowarn=true)
     sp = TS.default_space[]
     @test order(sp) == order() == 2
     @test get_numvars(sp) == get_numvars() == 66
@@ -33,17 +33,17 @@ using Test
     @test sp.index_table[2][1] == 3^65
     @test sp.pos_table[2][3^64] == 2
 
-    @test eltype(variables!("a", order=1, numvars=2))  == TaylorN{Float64}
+    @test eltype(variables!("a", order=1, numvars=2, nowarn=true))  == TaylorN{Float64}
     @test get_variable_names()[1] == "a₁"
     @test get_variable_symbols()[2] == :a₂
-    @test eltype(variables!(Int, "x", numvars=2, order=6)) == TaylorN{Int}
-    @test eltype(variables!("x", numvars=2, order=6)) == TaylorN{Float64}
-    @test eltype(variables!(BigInt, "x y", order=6)) == TaylorN{BigInt}
-    @test eltype(variables!("x y", order=6)) == TaylorN{Float64}
-    @test eltype(variables!(Int, :x, numvars=2, order=6)) == TaylorN{Int}
-    @test eltype(variables!(:x, numvars=2, order=6)) == TaylorN{Float64}
-    @test eltype(variables!(BigInt, [:x, :y], order=6)) == TaylorN{BigInt}
-    @test eltype(variables!([:x, :y], order=6)) == TaylorN{Float64}
+    @test eltype(variables!(Int, "x", numvars=2, order=6, nowarn=true)) == TaylorN{Int}
+    @test eltype(variables!("x", numvars=2, order=6, nowarn=true)) == TaylorN{Float64}
+    @test eltype(variables!(BigInt, "x y", order=6, nowarn=true)) == TaylorN{BigInt}
+    @test eltype(variables!("x y", order=6, nowarn=true)) == TaylorN{Float64}
+    @test eltype(variables!(Int, :x, numvars=2, order=6, nowarn=true)) == TaylorN{Int}
+    @test eltype(variables!(:x, numvars=2, order=6, nowarn=true)) == TaylorN{Float64}
+    @test eltype(variables!(BigInt, [:x, :y], order=6, nowarn=true)) == TaylorN{BigInt}
+    @test eltype(variables!([:x, :y], order=6, nowarn=true)) == TaylorN{Float64}
     @test typeof(show_params_TaylorN()) == Nothing
     @test typeof(show_monomials(2)) == Nothing
 
@@ -100,13 +100,13 @@ end
     old_default_space = TS.default_space[]
     dx, dy = variables()
     new_default_space = JetSpace(order=4, variables=[:δx, :δy, :δz])
-    @test TS.set_default_space!(new_default_space) === new_default_space
+    @test TS.set_default_space!(new_default_space; nowarn=true) === new_default_space
     @test TS.default_space[] === new_default_space
     @test dx.space === old_default_space
     @test dx.space !== TS.default_space[]
     @test sprint(show, dx) isa String
 
-    u, v = variables!("u v", order=2)
+    u, v = variables!("u v", order=2, nowarn=true)
     @test u.space === v.space
     @test u.space !== sx
     @test x*y == TaylorN(sx, HomogeneousPolynomial(sx, [0.0, 1.0, 0.0, 0.0, 0.0, 0.0], 2), 5)
@@ -125,7 +125,7 @@ end
     @test HomogeneousPolynomial{Int} <: AbstractSeries{Int}
     @test TaylorN{Float64} <: AbstractSeries{Float64}
 
-    variables!([:x, :y], order=6)
+    variables!([:x, :y], order=6, nowarn=true)
     @test order() == 6
     @test get_numvars() == 2
 
@@ -135,7 +135,7 @@ end
     @test variables(Int, 3)[1] == TaylorN(Int, 1, order=3)
     @test length(variables()) == get_numvars()
 
-    x, y = variables!("x y", order=6)
+    x, y = variables!("x y", order=6, nowarn=true)
     @test axes(x) == axes(y) == ()
     @test axes(x[1]) == axes(y[2]) == ()
     @test size(x) == (7,)
@@ -163,7 +163,7 @@ end
     @test y == HomogeneousPolynomial(2)
     @test !isnan(x)
 
-    variables!("x", numvars=2, order=17)
+    variables!("x", numvars=2, order=17, nowarn=true)
     v = [1,2]
     vv = TS._coeffsHP(v, 2)
     @test vv isa FixedSizeVectorDefault{Int}
@@ -863,7 +863,7 @@ end
     @test ( -sinh(xT+yT)^2 + cosh(xT+yT)^2 )(rand(2)) == 1
     @test ( -sinh(xT+yT)^2 + cosh(xT+yT)^2 )(zeros(2)) == 1
     #number of variables changed to 4...
-    dx = variables!("x", numvars=4, order=10)
+    dx = variables!("x", numvars=4, order=10, nowarn=true)
     P = sin.(dx)
     v = [1.0,2,3,4]
     for i in 1:4
@@ -885,7 +885,7 @@ end
     @test Q() == evaluate(Q)
     @test Q[1:3]() == evaluate(Q[1:3])
 
-    dx = variables!("x", numvars=4, order=10)
+    dx = variables!("x", numvars=4, order=10, nowarn=true)
     for i in 1:4
         @test deg2rad(180+dx[i]) == pi + deg2rad(1.0)dx[i]
         @test rad2deg(pi+dx[i]) == 180.0+rad2deg(1.0)dx[i]
@@ -924,7 +924,7 @@ end
     end
 
     @testset "Test evaluate! method for arrays of TaylorN" begin
-        x = variables!("x", order=2, numvars=4)
+        x = variables!("x", order=2, numvars=4, nowarn=true)
         function radntn!(y)
             for k in eachindex(y)
                 for l in eachindex(y[k])
@@ -954,7 +954,7 @@ end
 
 @testset "Integrate for several variables" begin
 
-    t, x, y = variables!("t x y")
+    t, x, y = variables!("t x y", nowarn=true)
 
     @test integrate(t, 1) == 0.5*t^2
     @test integrate(t, 2) == t * x
@@ -967,7 +967,7 @@ end
 
 @testset "Consistency of coeff_table" begin
     order = 20
-    x, y, z, w = variables!(Int128, "x y z w", numvars=4, order=2order)
+    x, y, z, w = variables!(Int128, "x y z w", numvars=4, order=2order, nowarn=true)
     sp = x.space
     ctab = deepcopy(sp.coeff_table);
 
diff --git a/test/mixtures.jl b/test/mixtures.jl
index 30e50eb3..89245a1a 100644
--- a/test/mixtures.jl
+++ b/test/mixtures.jl
@@ -11,7 +11,7 @@ using Test
     @test TS.NumberNotSeries == Union{Real,Complex}
     @test TS.NumberNotSeriesN == Union{Real,Complex,Taylor1}
 
-    variables!("x", numvars=2, order=6)
+    variables!("x", numvars=2, order=6, nowarn=true)
     xH = HomogeneousPolynomial(Int, 1)
     yH = HomogeneousPolynomial(Int, 2)
     tN = Taylor1(TaylorN{Float64}, 3)
@@ -231,7 +231,7 @@ using Test
     @test q1N(-xN^2) == 1.0 - xN^2*yN
 
     # See #92 and #94
-    δx, δy = variables!("δx δy")
+    δx, δy = variables!("δx δy", nowarn=true)
     xx = 1+Taylor1(δx, 5)
     yy = 1+Taylor1(δy, 5)
     tt = Taylor1([zero(δx), one(δx)], order(xx))
@@ -348,7 +348,7 @@ using Test
     p = cos(t)
     q = sin(t)
     a = [p,q]
-    dx = variables!("x", numvars=4, order=10)
+    dx = variables!("x", numvars=4, order=10, nowarn=true)
     P = sin.(dx)
     v = [1.0,2,3,4]
     F(x) = [sin(sin(x[4]+x[3])), sin(cos(x[3]-x[2])), cos(sin(x[1]^2+x[2]^2)), cos(cos(x[2]*x[3]))]
@@ -404,7 +404,7 @@ using Test
     diff_a11b11 = a11(t)-b11
     @test norm(diff_a11b11.coeffs, Inf) < 1E-19
 
-    X, Y = variables!(Taylor1{Float64}, "x y")
+    X, Y = variables!(Taylor1{Float64}, "x y", nowarn=true)
     @test typeof( norm(X) ) == Float64
     @test norm(X) > 0
     @test norm(X+Y) == sqrt(2)
@@ -455,7 +455,7 @@ using Test
     @test !isapprox(Q, P, atol=eps(), rtol=0)
     @test P ≉ Q^2
 
-    X, Y = variables!(BigFloat, "x y", numvars=2, order=6)
+    X, Y = variables!(BigFloat, "x y", numvars=2, order=6, nowarn=true)
     p1N = Taylor1([X^2,X*Y,Y+X,Y^2])
     q1N = Taylor1([X^2,(1.0+sqrt(eps(BigFloat)))*X*Y,Y+X,Y^2])
     @test p1N ≈ p1N
diff --git a/test/mutatingfuncts.jl b/test/mutatingfuncts.jl
index fb9be7ef..2a85c9d9 100644
--- a/test/mutatingfuncts.jl
+++ b/test/mutatingfuncts.jl
@@ -35,6 +35,16 @@ using Test
     @test t2[2] == 1.0
     #
     res = zero(t1)
+    TS.add!(res, t1, t2)
+    @test res == t1 + t2
+    TS.subst!(res, t1, t2)
+    @test res == t1 - t2
+    TS.mul!(res, t1, t2)
+    @test res == t1 * t2
+    TS.zero!(res)
+    TS.muladd!(res, t1, t2)
+    @test res == t1 * t2
+
     TS.add!(res, t1, t2, 3)
     @test res[3] == 0.0
     TS.add!(res, 1, t2, 3)
diff --git a/test/runtests.jl b/test/runtests.jl
index d9cd11cb..a84e13a3 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -23,4 +23,4 @@ for file in testfiles
 end
 
 # Back to TaylorN default parameters
-variables!("x", order=6, numvars=2);
+variables!("x", order=6, numvars=2, nowarn=true);
diff --git a/test/staticarrays.jl b/test/staticarrays.jl
index b9c0019d..79c2e97a 100644
--- a/test/staticarrays.jl
+++ b/test/staticarrays.jl
@@ -6,7 +6,7 @@ using TaylorSeries, StaticArrays
 using Test
 
 @testset "Tests TaylorSeries operations over StaticArrays types" begin
-    q = variables!("q", order=2, numvars=2)
+    q = variables!("q", order=2, numvars=2, nowarn=true)
     m = @SMatrix fill(Taylor1(rand(2).*q), 3, 3)
     mt = m'
     @test m isa SMatrix{3, 3, Taylor1{TaylorN{Float64}}, 9}

From eb911aa0fe3c44349f51bfdb4b1f06dfa1d18bc8 Mon Sep 17 00:00:00 2001
From: Jorge Perez-Hernandez <jorge.perez_hernandez@cfa.harvard.edu>
Date: Fri, 29 May 2026 17:45:00 -0400
Subject: [PATCH 02/20] Bump patch version

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 85d67f6d..a916bbd9 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "TaylorSeries"
 uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
-version = "0.22.0"
+version = "0.22.1"
 repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
 
 [deps]

From 85a0ec18f30e01ae66258f9356902c162ec41e70 Mon Sep 17 00:00:00 2001
From: Jorge Perez-Hernandez <jorge.perez_hernandez@cfa.harvard.edu>
Date: Sat, 30 May 2026 22:50:33 -0400
Subject: [PATCH 03/20] Improvements for nested Taylor1s

---
 src/arithmetic.jl      | 198 +++++++++++++++++++++++++++++++++++++++++
 src/power.jl           |  85 ++++++++++++++++++
 test/mutatingfuncts.jl |  22 +++++
 3 files changed, 305 insertions(+)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index aa093827..8aa6f639 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -339,6 +339,28 @@ function subst!(v::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}) where
     return nothing
 end
 
+function add!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
+        b::Taylor1{Taylor1{T}}) where {T<:NumberNotSeries}
+    v_coeffs = v.coeffs
+    a_coeffs = a.coeffs
+    b_coeffs = b.coeffs
+    @inbounds for i in eachindex(v_coeffs)
+        add!(v_coeffs[i], a_coeffs[i], b_coeffs[i])
+    end
+    return nothing
+end
+
+function subst!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
+        b::Taylor1{Taylor1{T}}) where {T<:NumberNotSeries}
+    v_coeffs = v.coeffs
+    a_coeffs = a.coeffs
+    b_coeffs = b.coeffs
+    @inbounds for i in eachindex(v_coeffs)
+        subst!(v_coeffs[i], a_coeffs[i], b_coeffs[i])
+    end
+    return nothing
+end
+
 function +(a::Taylor1{T}, b::Taylor1{T}) where {T<:NumberNotSeries}
     if order(a) != order(b)
         a, b = fixorder(a, b)
@@ -357,6 +379,26 @@ function -(a::Taylor1{T}, b::Taylor1{T}) where {T<:NumberNotSeries}
     return c
 end
 
+function +(a::Taylor1{Taylor1{T}}, b::Taylor1{Taylor1{T}}) where
+        {T<:NumberNotSeries}
+    if order(a) != order(b)
+        a, b = fixorder(a, b)
+    end
+    c = zero(a)
+    add!(c, a, b)
+    return c
+end
+
+function -(a::Taylor1{Taylor1{T}}, b::Taylor1{Taylor1{T}}) where
+        {T<:NumberNotSeries}
+    if order(a) != order(b)
+        a, b = fixorder(a, b)
+    end
+    c = zero(a)
+    subst!(c, a, b)
+    return c
+end
+
 for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
     @eval begin
         function ($f)(a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}}) where
@@ -585,6 +627,18 @@ function mul!(v::Taylor1{T}, a::Taylor1{S}, b::NumberNotSeries, k::Int) where
 end
 mul!(v::Taylor1{T}, a::NumberNotSeries, b::Taylor1{S}, k::Int) where
         {T<:NumberNotSeries, S<:NumberNotSeries} = mul!(v, b, a, k)
+
+function mul!(v::Taylor1{T}, a::Taylor1{S}, b::NumberNotSeries) where
+        {T<:NumberNotSeries, S<:NumberNotSeries}
+    v_coeffs = v.coeffs
+    a_coeffs = a.coeffs
+    @inbounds for i in eachindex(v_coeffs)
+        v_coeffs[i] = a_coeffs[i] * b
+    end
+    return nothing
+end
+mul!(v::Taylor1{T}, a::NumberNotSeries, b::Taylor1{S}) where
+        {T<:NumberNotSeries, S<:NumberNotSeries} = mul!(v, b, a)
 #
 function muladd!(c::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}, k::Int) where
         {T<:NumberNotSeries}
@@ -643,6 +697,21 @@ function mul_scalar!(c::Taylor1{T}, scalar::NumberNotSeries, a::Taylor1{T},
     return nothing
 end
 
+function mul_scalar!(c::Taylor1{T}, scalar::NumberNotSeries, a::Taylor1{T},
+        b::Taylor1{T}) where {T<:NumberNotSeries}
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    b_coeffs = b.coeffs
+    @inbounds for kk in eachindex(c_coeffs)
+        acc = zero(c_coeffs[kk])
+        for i = 1:kk
+            acc += a_coeffs[i] * b_coeffs[kk-i+1]
+        end
+        c_coeffs[kk] = scalar * acc
+    end
+    return nothing
+end
+
 # NOTE: For TaylorN, `mul!` (`muladd!`) *accumulates* the result of a * b in c[k]
 mul!(c::TaylorN{T}, a::TaylorN{T}, b::TaylorN{T}, k::Int) where
         {T<:Number} = muladd!(c, a, b, k)
@@ -678,6 +747,32 @@ function mul_scalar!(c::TaylorN{T}, scalar::NumberNotSeries, a::TaylorN{T},
 end
 
 # Nested Taylor1s
+function add!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
+        b::Taylor1{Taylor1{T}}, k::Int) where {T<:NumberNotSeries}
+    @inbounds add!(v.coeffs[k+1], a.coeffs[k+1], b.coeffs[k+1])
+    return nothing
+end
+
+function subst!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
+        b::Taylor1{Taylor1{T}}, k::Int) where {T<:NumberNotSeries}
+    @inbounds subst!(v.coeffs[k+1], a.coeffs[k+1], b.coeffs[k+1])
+    return nothing
+end
+
+function mul!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, b::Taylor1{Taylor1{T}},
+        k::Int) where {T<:NumberNotSeries}
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    b_coeffs = b.coeffs
+    kk = k+1
+    @inbounds c_k = c_coeffs[kk]
+    zero!(c_k)
+    @inbounds for i = 1:kk
+        muladd!(c_k, a_coeffs[i], b_coeffs[kk-i+1])
+    end
+    return nothing
+end
+
 function mul!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, b::Taylor1{Taylor1{T}},
         k::Int) where {T<:NumberNotSeriesN}
     @inbounds for j in eachindex(c[k])
@@ -688,8 +783,21 @@ function mul!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, b::Taylor1{Taylor1
     end
     return nothing
 end
+
 mul!(v::Taylor1{Taylor1{T}}, a::Taylor1{T}, b::Taylor1{Taylor1{T}}, k::Int) where
         {T<:NumberNotSeriesN} = mul!(v, b, a, k)
+function mul!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, b::Taylor1{T}, k::Int) where
+        {T<:NumberNotSeries}
+    @inbounds v_k = v.coeffs[k+1]
+    @inbounds a_k = a.coeffs[k+1]
+    if v_k === a_k
+        mul!(v_k, b)
+    else
+        mul!(v_k, a_k, b)
+    end
+    return nothing
+end
+
 function mul!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, b::Taylor1{T}, k::Int) where
         {T<:NumberNotSeriesN}
     @inbounds for i in eachindex(v[k])
@@ -697,8 +805,15 @@ function mul!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, b::Taylor1{T}, k::
     end
     return nothing
 end
+
 mul!(v::Taylor1{Taylor1{T}}, a::NumberNotSeries, b::Taylor1{Taylor1{T}}, k::Int) where
         {T<:NumberNotSeriesN} = mul!(v, b, a, k)
+function mul!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, b::NumberNotSeries,
+        k::Int) where {T<:NumberNotSeries}
+    @inbounds mul!(v.coeffs[k+1], a.coeffs[k+1], b)
+    return nothing
+end
+
 function mul!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, b::NumberNotSeries,
         k::Int) where {T<:NumberNotSeriesN}
     @inbounds for i in eachindex(v[k])
@@ -706,6 +821,20 @@ function mul!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, b::NumberNotSeries
     end
     return nothing
 end
+
+function muladd!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
+        b::Taylor1{Taylor1{T}}, k::Int) where {T<:NumberNotSeries}
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    b_coeffs = b.coeffs
+    kk = k+1
+    @inbounds c_k = c_coeffs[kk]
+    @inbounds for i = 1:kk
+        muladd!(c_k, a_coeffs[i], b_coeffs[kk-i+1])
+    end
+    return nothing
+end
+
 function muladd!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
         b::Taylor1{Taylor1{T}}, k::Int) where {T<:NumberNotSeriesN}
     @inbounds for j in eachindex(c[k])
@@ -715,6 +844,7 @@ function muladd!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
     end
     return nothing
 end
+
 # muladd!(v::Taylor1{Taylor1{T}}, a::Taylor1{T}, b::Taylor1{Taylor1{T}}, k::Int) where
 #         {T<:NumberNotSeriesN} = muladd!(v, b, a, k)
 # function muladd!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, b::Taylor1{T}, k::Int) where
@@ -733,6 +863,25 @@ end
 # end
 # muladd!(v::Taylor1{Taylor1{T}}, a::T, b::Taylor1{Taylor1{T}}, k::Int) where
 #     {T<:NumberNotSeriesN} = muladd!(v, b, a, k)
+function mul_scalar!(c::Taylor1{Taylor1{T}}, scalar::NumberNotSeries,
+        a::Taylor1{Taylor1{T}}, b::Taylor1{Taylor1{T}}, k::Int) where
+        {T<:NumberNotSeries}
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    b_coeffs = b.coeffs
+    kk = k+1
+    @inbounds c_k = c_coeffs[kk]
+    zero!(c_k)
+    @inbounds for i = 1:kk
+        muladd!(c_k, a_coeffs[i], b_coeffs[kk-i+1])
+    end
+    c_k_coeffs = c_k.coeffs
+    @inbounds for i in eachindex(c_k_coeffs)
+        c_k_coeffs[i] *= scalar
+    end
+    return nothing
+end
+
 function mul_scalar!(c::Taylor1{Taylor1{T}}, scalar::NumberNotSeries,
         a::Taylor1{Taylor1{T}}, b::Taylor1{Taylor1{T}}, k::Int) where {T<:Number}
     mul!(c, a, b, k)
@@ -743,6 +892,26 @@ function mul_scalar!(c::Taylor1{Taylor1{T}}, scalar::NumberNotSeries,
     return nothing
 end
 
+function mul_scalar!(c::Taylor1{Taylor1{T}}, scalar::NumberNotSeries,
+        a::Taylor1{Taylor1{T}}, b::Taylor1{Taylor1{T}}) where
+        {T<:NumberNotSeries}
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    b_coeffs = b.coeffs
+    @inbounds for kk in eachindex(c_coeffs)
+        c_k = c_coeffs[kk]
+        zero!(c_k)
+        for i = 1:kk
+            muladd!(c_k, a_coeffs[i], b_coeffs[kk-i+1])
+        end
+        c_k_coeffs = c_k.coeffs
+        for i in eachindex(c_k_coeffs)
+            c_k_coeffs[i] *= scalar
+        end
+    end
+    return nothing
+end
+
 
 # for T in (:Taylor1, :TaylorN)
 #     @eval begin
@@ -792,6 +961,20 @@ function mul!(a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}}) where
     end
     return nothing
 end
+function mul!(a::Taylor1{Taylor1{T}}, b::Taylor1{Taylor1{T}}) where
+        {T<:NumberNotSeries}
+    a_coeffs = a.coeffs
+    b_coeffs = b.coeffs
+    @inbounds for kk in reverse(eachindex(a_coeffs))
+        a_k = a_coeffs[kk]
+        mul!(a_k, b_coeffs[1])
+        for l = 2:kk
+            muladd!(a_k, a_coeffs[kk-l+1], b_coeffs[l])
+        end
+    end
+    return nothing
+end
+
 function mul!(a::Taylor1{Taylor1{T}}, b::Taylor1{Taylor1{T}}) where
         {T<:NumberNotSeriesN}
     @inbounds for k in reverse(eachindex(a))
@@ -862,6 +1045,21 @@ function mul!(c::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}) where {T<:Number}
     end
 end
 
+function mul!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
+        b::Taylor1{Taylor1{T}}) where {T<:NumberNotSeries}
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    b_coeffs = b.coeffs
+    @inbounds for kk in eachindex(c_coeffs)
+        c_k = c_coeffs[kk]
+        zero!(c_k)
+        for i = 1:kk
+            muladd!(c_k, a_coeffs[i], b_coeffs[kk-i+1])
+        end
+    end
+    return nothing
+end
+
 function mul!(c::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
         b::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
     _check_same_space(c[0], a[0], b[0])
diff --git a/src/power.jl b/src/power.jl
index d0934dac..28a30364 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -362,6 +362,61 @@ function pow!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
     return nothing
 end
 
+function pow!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
+        aux::Taylor1{Taylor1{T}}, r::S, k::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, constant_term(aux), k)
+    (r == 0.5) && return sqrt!(c, a, aux, k)
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    aux_coeffs = aux.coeffs
+    kk = k+1
+    @inbounds c_k = c_coeffs[kk]
+    zero!(c_k)
+    # First non-zero coefficient
+    l0 = findfirst(a)
+    l0 < 0 && return nothing
+    # 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 )
+    kprime = k-lnull
+    order_a = order(a)
+    (kprime < 0 || lnull > order_a) && return nothing
+    # Relevant for positive integer r, to avoid round-off errors
+    lastnz = isinteger(r) && r > 0 ? findlast(a) : order_a
+    isinteger(r) && r > 0 && (k > r*lastnz) && return nothing
+    @inbounds a_l0 = a_coeffs[l0+1]
+    if k == lnull
+        @inbounds aux0 = aux_coeffs[1]
+        @inbounds for j in eachindex(a_l0)
+            pow!(c_k, a_l0, aux0, float(r), j)
+        end
+        return nothing
+    end
+    # The recursion formula
+    @inbounds aux_k = aux_coeffs[kk]
+    ilo = max(0, l0+kprime-order_a)
+    ihi = min(k-lnull-1, order_a-lnull)
+    @inbounds for i = ilo:ihi
+        rr = r*(kprime-i) - i
+        mul_scalar!(aux_k, rr, c_coeffs[i+lnull+1], a_coeffs[l0+kprime-i+1])
+        add!(c_k, c_k, aux_k)
+    end
+    # c[k] = c[k] / (kprime * a[l0])
+    @inbounds copyto!(aux_k.coeffs, c_k.coeffs)
+    @inbounds for j in eachindex(a_l0)
+        div!(c_k, aux_k, a_l0, j)
+    end
+    @inbounds for j in eachindex(a_l0)
+        div!(c_k, c_k, kprime, j)
+    end
+    return nothing
+end
+
 function pow!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
         aux::Taylor1{Taylor1{T}}, r::S, k::Int) where
         {T<:NumberNotSeriesN, S<:Real}
@@ -625,6 +680,36 @@ function sqr!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, aux::TaylorN{T},
     return nothing
 end
 
+function sqr!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, aux::Taylor1{T},
+        k::Int) where {T<:NumberNotSeries}
+    if k == 0
+        sqr_orderzero!(c, a)
+        return nothing
+    end
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    kk = k+1
+    @inbounds c_k = c_coeffs[kk]
+    zero!(c_k)
+    # Recursion formula
+    kodd = k%2
+    kend = (k - 2 + kodd) >> 1
+    @inbounds for i = 0:kend
+        mul_scalar!(aux, 2, a_coeffs[i+1], a_coeffs[kk-i])
+        add!(c_k, c_k, aux)
+    end
+    kodd == 1 && return nothing
+    # c[k] += a[k >> 1]^2
+    aaux = zero(aux[0])
+    zero!(aux)
+    @inbounds a_mid = a_coeffs[(k >> 1)+1]
+    @inbounds for j in eachindex(a_mid)
+        sqr!(aux, a_mid, aaux, j)
+    end
+    add!(c_k, c_k, aux)
+    return nothing
+end
+
 function sqr!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, aux::Taylor1{T},
         k::Int) where {T<:Number}
     if k == 0
diff --git a/test/mutatingfuncts.jl b/test/mutatingfuncts.jl
index 2a85c9d9..6c6eeb56 100644
--- a/test/mutatingfuncts.jl
+++ b/test/mutatingfuncts.jl
@@ -103,4 +103,26 @@ using Test
     TaylorSeries.zero!(t2)
     @test TaylorSeries.iszero(t2)
 
+    outer_order = 5
+    inner_order = 4
+    nt1 = Taylor1([Taylor1([0.5 + i + 0.1*j for j in 0:inner_order],
+        inner_order) for i in 0:outer_order], outer_order)
+    nt2 = Taylor1([Taylor1([1.0 + 0.2*i - 0.03*j for j in 0:inner_order],
+        inner_order) for i in 0:outer_order], outer_order)
+    nres = zero(nt1)
+    TS.add!(nres, nt1, nt2)
+    @test nres == nt1 + nt2
+    TS.subst!(nres, nt1, nt2)
+    @test nres == nt1 - nt2
+    TS.mul!(nres, nt1, nt2)
+    @test nres ≈ nt1 * nt2
+    TS.zero!(nres)
+    for k in eachindex(nres)
+        TS.muladd!(nres, nt1, nt2, k)
+    end
+    @test nres ≈ nt1 * nt2
+    ninplace = deepcopy(nt1)
+    TS.mul!(ninplace, nt2)
+    @test ninplace ≈ nt1 * nt2
+
 end

From d0de25176de995a5d666cf3f32e8173c554527ca Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Sat, 6 Jun 2026 16:14:57 -0400
Subject: [PATCH 04/20] Switch to direct coeffs access for mutating methods of:
 elementary functions, arithmetic helpers, zero, one, identity, scalar sqrt!,
 nested differentiation/integration and some evaluation methods (Horner)

---
 src/arithmetic.jl |  67 +++++++++++-----
 src/calculus.jl   |  37 +++++++--
 src/evaluate.jl   |  35 +++++++--
 src/functions.jl  | 195 +++++++++++++++++++++++++++-------------------
 src/identity.jl   |  11 ++-
 src/power.jl      |  16 ++--
 src/zero_one.jl   |  28 ++++++-
 7 files changed, 261 insertions(+), 128 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 8aa6f639..d08f3cf8 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -1489,53 +1489,80 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`.
 # @inline
 function div!(c::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}, k::Int) where
         {T<:NumberNotSeries}
-    zero!(c, k)
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    b_coeffs = b.coeffs
+    kk = k+1
+    @inbounds c_coeffs[kk] = zero(c_coeffs[kk])
     iszero(a) && !iszero(b) && return nothing
     # order and coefficient of first factorized term
     ordfact, cdivfact = divfactorization(a, b)
     if k == 0
-        @inbounds c[0] = cdivfact
+        @inbounds c_coeffs[1] = cdivfact
         return nothing
     end
     b_order = order(b)
     imin = max(0, k+ordfact-b_order)
-    @inbounds c[k] = c[imin] * b[k+ordfact-imin]
+    @inbounds acc = c_coeffs[imin+1] * b_coeffs[k+ordfact-imin+1]
     @inbounds for i = imin+1:k-1
-        c[k] += c[i] * b[k+ordfact-i]
+        acc += c_coeffs[i+1] * b_coeffs[k+ordfact-i+1]
     end
     if k+ordfact ≤ b_order
-        @inbounds c[k] = a[k+ordfact] - c[k]
+        @inbounds acc = a_coeffs[k+ordfact+1] - acc
     else
-        @inbounds c[k] = - c[k]
+        acc = -acc
     end
-    c[k] = c[k] / b[ordfact]
+    @inbounds c_coeffs[kk] = acc / b_coeffs[ordfact+1]
     return nothing
 end
 
 # @inline
 function div!(v::Taylor1{T}, a::Taylor1{T}, b::NumberNotSeries,
         k::Int) where {T<:Number}
-    @inbounds v[k] = a[k] / b
+    @inbounds v.coeffs[k+1] = a.coeffs[k+1] / b
+    return nothing
+end
+
+function div!(v::Taylor1{T}, a::Taylor1{S}, b::NumberNotSeries) where
+        {T<:NumberNotSeries, S<:NumberNotSeries}
+    v_coeffs = v.coeffs
+    a_coeffs = a.coeffs
+    @inbounds for i in eachindex(v_coeffs)
+        v_coeffs[i] = a_coeffs[i] / b
+    end
+    return nothing
+end
+
+function div!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
+        b::NumberNotSeries) where {T<:NumberNotSeries}
+    v_coeffs = v.coeffs
+    a_coeffs = a.coeffs
+    @inbounds for i in eachindex(v_coeffs)
+        div!(v_coeffs[i], a_coeffs[i], b)
+    end
     return nothing
 end
 
 # @inline
 function div!(c::Taylor1{T}, a::NumberNotSeries, b::Taylor1{T}, k::Int) where
         {T<:NumberNotSeries}
-    zero!(c, k)
+    c_coeffs = c.coeffs
+    b_coeffs = b.coeffs
+    kk = k+1
+    @inbounds c_coeffs[kk] = zero(c_coeffs[kk])
     iszero(a) && !iszero(b) && return nothing
     # order and coefficient of first factorized term
     # In this case, since a[k]=0 for k>0, we can simplify to:
     # ordfact, cdivfact = 0, a/b[0]
     if k == 0
-        @inbounds c[0] = a / b[0]
+        @inbounds c_coeffs[1] = a / b_coeffs[1]
         return nothing
     end
-    @inbounds c[k] = c[0] * b[k]
+    @inbounds acc = c_coeffs[1] * b_coeffs[kk]
     @inbounds for i = 1:k-1
-        c[k] += c[i] * b[k-i]
+        acc += c_coeffs[i+1] * b_coeffs[k-i+1]
     end
-    @inbounds c[k] = -c[k] / b[0]
+    @inbounds c_coeffs[kk] = -acc / b_coeffs[1]
     return nothing
 end
 
@@ -1723,18 +1750,20 @@ end
 
 @inline function div_scalar!(c::Taylor1{T}, scalar::NumberNotSeries,
         a::Taylor1{T}, k::Int) where {T <: NumberNotSeries}
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
     if k==0
-        @inbounds c[0] = scalar*constant_term(c) / constant_term(a)
+        @inbounds c_coeffs[1] = scalar*c_coeffs[1] / a_coeffs[1]
         return nothing
     end
 
-    aux = c[k] = scalar * c[k]
-    @inbounds zero!(c, k)
+    kk = k+1
+    @inbounds aux = scalar * c_coeffs[kk]
+    @inbounds acc = zero(c_coeffs[kk])
     @inbounds for i = 0:k-1
-        c[k] -= c[i] * a[k-i]
+        acc -= c_coeffs[i+1] * a_coeffs[k-i+1]
     end
-    c[k] += aux
-    @inbounds c[k] = c[k] / constant_term(a)
+    @inbounds c_coeffs[kk] = (acc + aux) / a_coeffs[1]
     return nothing
 end
 
diff --git a/src/calculus.jl b/src/calculus.jl
index 9412e230..728bab68 100644
--- a/src/calculus.jl
+++ b/src/calculus.jl
@@ -16,7 +16,7 @@ The order of the result is `order(a)-1`.
 The function `derivative` is an exact synonym of `differentiate`.
 """
 function differentiate(a::Taylor1)
-    res = Taylor1(zero(a[0]), order(a)-1)
+    res = Taylor1(zero(a.coeffs[1]), order(a)-1)
     for ord in eachindex(res)
         differentiate!(res, a, ord)
     end
@@ -60,7 +60,14 @@ p_k = (k+1) a_{k+1}.
 function differentiate!(p::Taylor1{T}, a::Taylor1{T}, k::Int) where
         {T<:NumberNotSeries}
     k >= order(a) && return nothing
-    @inbounds p[k] = (k+1)*a[k+1]
+    @inbounds p.coeffs[k+1] = (k+1)*a.coeffs[k+2]
+    return nothing
+end
+
+function differentiate!(p::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
+        k::Int) where {T<:NumberNotSeries}
+    k >= order(a) && return nothing
+    @inbounds mul!(p.coeffs[k+1], a.coeffs[k+2], k+1)
     return nothing
 end
 function differentiate!(p::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, k::Int) where
@@ -111,7 +118,7 @@ Return the value of the `n`-th differentiate of the polynomial `a`.
 """
 function differentiate(n::Int, a::Taylor1{T}) where {T<:Number}
     @assert order(a) ≥ n ≥ 0
-    return factorial( widen(n) ) * a[n] :: T
+    return factorial( widen(n) ) * a.coeffs[n+1] :: T
 end
 
 
@@ -125,24 +132,38 @@ Note that the order of the result is `order(a)+1`.
 """
 function integrate(a::Taylor1{T}, x::S) where {T<:Number, S<:Number}
     order = TS.order(a)
-    aa = zero(a[0])/1 + zero(x)
+    aa = zero(a.coeffs[1])/1 + zero(x)
     res = Taylor1(aa, order)
     integrate!(res, a, convert(typeof(aa), x))
     return res
 end
-integrate(a::Taylor1{T}) where {T<:Number} = integrate(a, zero(a[0]))
+integrate(a::Taylor1{T}) where {T<:Number} = integrate(a, zero(a.coeffs[1]))
 
 # In-place method
 function integrate!(res::Taylor1{T}, a::Taylor1{T}, x::T) where {T<:Number}
     order = TS.order(a)
+    res_coeffs = res.coeffs
+    a_coeffs = a.coeffs
     @inbounds for i = 1:order
-        res.coeffs[i+1] = a[i-1] / i
+        res_coeffs[i+1] = a_coeffs[i] / i
     end
-    @inbounds res.coeffs[1] = x
+    @inbounds res_coeffs[1] = x
     return nothing
 end
 integrate!(res::Taylor1{T}, a::Taylor1{T}) where {T<:Number} =
-    integrate!(res, a, zero(a[0]))
+    integrate!(res, a, zero(a.coeffs[1]))
+
+function integrate!(res::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
+        x::Taylor1{T}) where {T<:NumberNotSeries}
+    order = TS.order(a)
+    res_coeffs = res.coeffs
+    a_coeffs = a.coeffs
+    @inbounds for i = 1:order
+        div!(res_coeffs[i+1], a_coeffs[i], i)
+    end
+    @inbounds res_coeffs[1] = x
+    return nothing
+end
 
 
 
diff --git a/src/evaluate.jl b/src/evaluate.jl
index ebd88e2d..44676e80 100644
--- a/src/evaluate.jl
+++ b/src/evaluate.jl
@@ -15,23 +15,25 @@ omitted, its value is considered as zero. Note that the syntax `a(dx)` is
 equivalent to `evaluate(a,dx)`, and `a()` is equivalent to `evaluate(a)`.
 """
 function evaluate(a::Taylor1{T}, dx::S) where {T<:NumberNotSeries, S<:NumberNotSeries}
-    @inbounds suma = zero(a[end])
+    a_coeffs = a.coeffs
+    @inbounds suma = zero(a_coeffs[end])
     # suma = evalpoly(dx, view(a.coeffs,:) )
-    @inbounds for k in reverse(eachindex(a))
-        suma = suma * dx + a[k]
+    @inbounds for k in reverse(eachindex(a_coeffs))
+        suma = suma * dx + a_coeffs[k]
     end
     return suma
 end
 
 function evaluate(a::Taylor1{T}, dx::S) where {T<:Number, S<:Number}
-    suma = a[end]*zero(dx)
-    @inbounds for k in reverse(eachindex(a))
-        suma = suma * dx + a[k]
+    a_coeffs = a.coeffs
+    suma = a_coeffs[end]*zero(dx)
+    @inbounds for k in reverse(eachindex(a_coeffs))
+        suma = suma * dx + a_coeffs[k]
     end
     return suma
 end
 
-evaluate(a::Taylor1{T}) where {T<:Number} = a[0]
+evaluate(a::Taylor1{T}) where {T<:Number} = a.coeffs[1]
 
 
 """
@@ -57,11 +59,28 @@ Note that the syntax `a(x)` is equivalent to `evaluate(a, x)`.
 evaluate(a::Taylor1{T}, x::Taylor1{S}) where {T<:Number, S<:Number} =
     evaluate(promote(a, x)...)
 
+function evaluate(a::Taylor1{T}, x::Taylor1{T}) where {T<:NumberNotSeries}
+    if order(a) != order(x)
+        a, x = fixorder(a, x)
+    end
+    suma = zero(x)
+    aux = zero(x)
+    suma_coeffs = suma.coeffs
+    aux_coeffs = aux.coeffs
+    a_coeffs = a.coeffs
+    @inbounds for k in reverse(eachindex(a_coeffs))
+        mul!(aux, suma, x)
+        copyto!(suma_coeffs, aux_coeffs)
+        suma_coeffs[1] += a_coeffs[k]
+    end
+    return suma
+end
+
 function evaluate(a::Taylor1{T}, x::Taylor1{T}) where {T<:Number}
     if order(a) != order(x)
         a, x = fixorder(a, x)
     end
-    @inbounds suma = a[end]*zero(x)
+    @inbounds suma = a.coeffs[end]*zero(x)
     _horner!(suma, a, x, zero(suma))
     return suma
 end
diff --git a/src/functions.jl b/src/functions.jl
index 97871047..8f6ee70b 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -235,15 +235,18 @@ end
 # Taylor1 mutating functions
 function exp!(c::Taylor1{T}, a::Taylor1{T}, k::Int) where
         {T<:NumberNotSeries}
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
     if k == 0
-        @inbounds c[0] = exp(constant_term(a))
+        @inbounds c_coeffs[1] = exp(a_coeffs[1])
         return nothing
     end
-    zero!(c, k)
-    @inbounds for i = 0:k-1
-        c[k] += (k-i) * a[k-i] * c[i]
+    kk = k+1
+    @inbounds acc = zero(c_coeffs[kk])
+    @inbounds for i = 1:k
+        acc += (k-i+1) * a_coeffs[k-i+2] * c_coeffs[i]
     end
-    div!(c, c, k, k)
+    @inbounds c_coeffs[kk] = acc / k
     return nothing
 end
 function exp!(c::Taylor1{T}, a::Taylor1{T}, k::Int) where
@@ -308,13 +311,14 @@ end
 function expm1!(c::Taylor1{T}, a::Taylor1{T}, k::Int) where
         {T<:NumberNotSeries}
     if k == 0
-        @inbounds c[0] = expm1(constant_term(a))
+        @inbounds c.coeffs[1] = expm1(a.coeffs[1])
         return nothing
     end
-    cc = c[0]
-    c[0] += one(c[0])
+    c_coeffs = c.coeffs
+    @inbounds cc = c_coeffs[1]
+    @inbounds c_coeffs[1] += one(cc)
     exp!(c, a, k)
-    c[0] = cc
+    @inbounds c_coeffs[1] = cc
     return nothing
 end
 function expm1!(c::Taylor1{T}, a::Taylor1{T}, k::Int) where
@@ -334,18 +338,22 @@ end
 
 function log!(c::Taylor1{T}, a::Taylor1{T}, k::Int) where
         {T<:NumberNotSeries}
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    @inbounds a0 = a_coeffs[1]
     if k == 0
-        @inbounds c[0] = log(constant_term(a))
+        @inbounds c_coeffs[1] = log(a0)
         return nothing
     elseif k == 1
-        @inbounds c[1] = a[1] / constant_term(a)
+        @inbounds c_coeffs[2] = a_coeffs[2] / a0
         return nothing
     end
-    zero!(c, k)
+    kk = k+1
+    @inbounds acc = zero(c_coeffs[kk])
     @inbounds for i = 1:k-1
-        c[k] += (k-i) * a[i] * c[k-i]
+        acc += (k-i) * a_coeffs[i+1] * c_coeffs[k-i+1]
     end
-    @inbounds c[k] = (a[k] - c[k]/k) / constant_term(a)
+    @inbounds c_coeffs[kk] = (a_coeffs[kk] - acc/k) / a0
     return nothing
 end
 function log!(c::Taylor1{T}, a::Taylor1{T}, k::Int) where
@@ -388,20 +396,23 @@ end
 
 function log1p!(c::Taylor1{T}, a::Taylor1{T}, k::Int) where
         {T<:NumberNotSeries}
-    a0 = constant_term(a)
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    @inbounds a0 = a_coeffs[1]
     a0p1 = a0+one(a0)
     if k == 0
-        @inbounds c[0] = log1p(a0)
+        @inbounds c_coeffs[1] = log1p(a0)
         return nothing
     elseif k == 1
-        @inbounds c[1] = a[1] / a0p1
+        @inbounds c_coeffs[2] = a_coeffs[2] / a0p1
         return nothing
     end
-    zero!(c, k)
+    kk = k+1
+    @inbounds acc = zero(c_coeffs[kk])
     @inbounds for i = 1:k-1
-        c[k] += (k-i) * a[i] * c[k-i]
+        acc += (k-i) * a_coeffs[i+1] * c_coeffs[k-i+1]
     end
-    @inbounds c[k] = (a[k] - c[k]/k) / a0p1
+    @inbounds c_coeffs[kk] = (a_coeffs[kk] - acc/k) / a0p1
     return nothing
 end
 function log1p!(c::Taylor1{T}, a::Taylor1{T}, k::Int) where
@@ -443,19 +454,23 @@ end
 
 function sincos!(s::Taylor1{T}, c::Taylor1{T}, a::Taylor1{T}, k::Int) where
         {T<:NumberNotSeries}
+    s_coeffs = s.coeffs
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
     if k == 0
-        @inbounds s[0], c[0] = sincos( constant_term(a) )
+        @inbounds s_coeffs[1], c_coeffs[1] = sincos(a_coeffs[1])
         return nothing
     end
-    zero!(s, k)
-    zero!(c, k)
+    kk = k+1
+    @inbounds s_acc = zero(s_coeffs[kk])
+    @inbounds c_acc = zero(c_coeffs[kk])
     @inbounds for i = 1:k
-        x = i * a[i]
-        s[k] += x * c[k-i]
-        c[k] -= x * s[k-i]
+        x = i * a_coeffs[i+1]
+        s_acc += x * c_coeffs[k-i+1]
+        c_acc -= x * s_coeffs[k-i+1]
     end
-    s[k] = s[k] / k
-    c[k] = c[k] / k
+    @inbounds s_coeffs[kk] = s_acc / k
+    @inbounds c_coeffs[kk] = c_acc / k
     return nothing
 end
 function sincos!(s::Taylor1{T}, c::Taylor1{T}, a::Taylor1{T}, k::Int) where
@@ -489,20 +504,23 @@ end
 
 function sincospi!(s::Taylor1{T}, c::Taylor1{T}, a::Taylor1{T}, k::Int) where
         {T<:NumberNotSeries}
+    s_coeffs = s.coeffs
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
     if k == 0
-        a0 = constant_term(a)
-        @inbounds s[0], c[0] = sincospi( a0 )
+        @inbounds s_coeffs[1], c_coeffs[1] = sincospi(a_coeffs[1])
         return nothing
     end
-    zero!(s, k)
-    zero!(c, k)
+    kk = k+1
+    @inbounds s_acc = zero(s_coeffs[kk])
+    @inbounds c_acc = zero(c_coeffs[kk])
     @inbounds for i = 1:k
-        x = (i * pi) * a[i]
-        s[k] += x * c[k-i]
-        c[k] -= x * s[k-i]
+        x = (i * pi) * a_coeffs[i+1]
+        s_acc += x * c_coeffs[k-i+1]
+        c_acc -= x * s_coeffs[k-i+1]
     end
-    s[k] = s[k] / k
-    c[k] = c[k] / k
+    @inbounds s_coeffs[kk] = s_acc / k
+    @inbounds c_coeffs[kk] = c_acc / k
     return nothing
 end
 function sincospi!(s::Taylor1{T}, c::Taylor1{T}, a::Taylor1{T}, k::Int) where
@@ -537,21 +555,22 @@ end
 
 function tan!(c::Taylor1{T}, a::Taylor1{T}, c2::Taylor1{T}, k::Int) where
         {T<:NumberNotSeries}
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    c2_coeffs = c2.coeffs
     if k == 0
-        @inbounds aux = tan( constant_term(a) )
-        @inbounds c[0] = aux
-        @inbounds c2[0] = aux^2
+        @inbounds aux = tan(a_coeffs[1])
+        @inbounds c_coeffs[1] = aux
+        @inbounds c2_coeffs[1] = aux^2
         return nothing
     end
-    zero!(c, k)
-    @inbounds for i = 0:k-1
-        c[k] += (k-i) * a[k-i] * c2[i]
+    kk = k+1
+    @inbounds acc = zero(c_coeffs[kk])
+    @inbounds for i = 1:k
+        acc += (k-i+1) * a_coeffs[k-i+2] * c2_coeffs[i]
     end
-    # c[k] <- c[k]/k
-    div!(c, c, k, k)
-    # c[k] <- c[k] + a[k]
-    add!(c, a, c, k)
-    sqr!(c2, c, zero(c[0]), k)
+    @inbounds c_coeffs[kk] = acc / k + a_coeffs[kk]
+    @inbounds sqr!(c2, c, zero(c_coeffs[1]), k)
     return nothing
 end
 function tan!(c::Taylor1{T}, a::Taylor1{T}, c2::Taylor1{T}, k::Int) where
@@ -700,18 +719,22 @@ end
 
 function atan!(c::Taylor1{T}, a::Taylor1{T}, r::Taylor1{T}, k::Int) where
         {T<:NumberNotSeries}
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    r_coeffs = r.coeffs
     if k == 0
-        a0 = constant_term(a)
-        @inbounds c[0] = atan( a0 )
-        @inbounds r[0] = 1 + a0^2
+        @inbounds a0 = a_coeffs[1]
+        @inbounds c_coeffs[1] = atan(a0)
+        @inbounds r_coeffs[1] = 1 + a0^2
         return nothing
     end
-    zero!(c, k)
+    kk = k+1
+    @inbounds acc = zero(c_coeffs[kk])
     @inbounds for i in 1:k-1
-        c[k] += (k-i) * r[i] * c[k-i]
+        acc += (k-i) * r_coeffs[i+1] * c_coeffs[k-i+1]
     end
-    sqr!(r, a, zero(a[0]), k)
-    @inbounds c[k] = (a[k] - c[k]/k) / constant_term(r)
+    @inbounds sqr!(r, a, zero(a_coeffs[1]), k)
+    @inbounds c_coeffs[kk] = (a_coeffs[kk] - acc/k) / r_coeffs[1]
     return nothing
 end
 function atan!(c::Taylor1{T}, a::Taylor1{T}, r::Taylor1{T}, k::Int) where
@@ -748,21 +771,25 @@ end
 
 function sinhcosh!(s::Taylor1{T}, c::Taylor1{T}, a::Taylor1{T}, k::Int) where
         {T<:NumberNotSeries}
+    s_coeffs = s.coeffs
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
     if k == 0
-        @inbounds s[0] = sinh( constant_term(a) )
-        @inbounds c[0] = cosh( constant_term(a) )
+        @inbounds a0 = a_coeffs[1]
+        @inbounds s_coeffs[1] = sinh(a0)
+        @inbounds c_coeffs[1] = cosh(a0)
         return nothing
     end
-    x = a[1]
-    zero!(s, k)
-    zero!(c, k)
+    kk = k+1
+    @inbounds s_acc = zero(s_coeffs[kk])
+    @inbounds c_acc = zero(c_coeffs[kk])
     @inbounds for i = 1:k
-        x = i * a[i]
-        s[k] += x * c[k-i]
-        c[k] += x * s[k-i]
+        x = i * a_coeffs[i+1]
+        s_acc += x * c_coeffs[k-i+1]
+        c_acc += x * s_coeffs[k-i+1]
     end
-    @inbounds div!(s, s, k, k)
-    @inbounds div!(c, c, k, k)
+    @inbounds s_coeffs[kk] = s_acc / k
+    @inbounds c_coeffs[kk] = c_acc / k
     return nothing
 end
 function sinhcosh!(s::Taylor1{T}, c::Taylor1{T}, a::Taylor1{T}, k::Int) where
@@ -793,18 +820,22 @@ end
 
 function tanh!(c::Taylor1{T}, a::Taylor1{T}, c2::Taylor1{T}, k::Int) where
         {T<:NumberNotSeries}
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    c2_coeffs = c2.coeffs
     if k == 0
-        @inbounds aux = tanh( constant_term(a) )
-        @inbounds c[0] = aux
-        @inbounds c2[0] = aux^2
+        @inbounds aux = tanh(a_coeffs[1])
+        @inbounds c_coeffs[1] = aux
+        @inbounds c2_coeffs[1] = aux^2
         return nothing
     end
-    zero!(c, k)
-    @inbounds for i = 0:k-1
-        c[k] += (k-i) * a[k-i] * c2[i]
+    kk = k+1
+    @inbounds acc = zero(c_coeffs[kk])
+    @inbounds for i = 1:k
+        acc += (k-i+1) * a_coeffs[k-i+2] * c2_coeffs[i]
     end
-    @inbounds c[k] = a[k] - c[k]/k
-    sqr!(c2, c, zero(c[0]), k)
+    @inbounds c_coeffs[kk] = a_coeffs[kk] - acc/k
+    @inbounds sqr!(c2, c, zero(c_coeffs[1]), k)
     return nothing
 end
 function tanh!(c::Taylor1{T}, a::Taylor1{T}, c2::Taylor1{T}, k::Int) where
@@ -933,18 +964,22 @@ end
 
 function atanh!(c::Taylor1{T}, a::Taylor1{T}, r::Taylor1{T}, k::Int) where
         {T<:NumberNotSeries}
-    a0 = constant_term(a)
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    r_coeffs = r.coeffs
+    @inbounds a0 = a_coeffs[1]
     if k == 0
-        @inbounds c[0] = atanh( a0 )
-        @inbounds r[0] = 1 - a0^2
+        @inbounds c_coeffs[1] = atanh(a0)
+        @inbounds r_coeffs[1] = 1 - a0^2
         return nothing
     end
-    zero!(c, k)
+    kk = k+1
+    @inbounds acc = zero(c_coeffs[kk])
     @inbounds for i in 1:k-1
-        c[k] += (k-i) * r[i] * c[k-i]
+        acc += (k-i) * r_coeffs[i+1] * c_coeffs[k-i+1]
     end
-    sqr!(r, a, zero(a0), k)
-    @inbounds c[k] = (a[k] + c[k]/k) / constant_term(r)
+    @inbounds sqr!(r, a, zero(a0), k)
+    @inbounds c_coeffs[kk] = (a_coeffs[kk] + acc/k) / r_coeffs[1]
     return nothing
 end
 function atanh!(c::Taylor1{T}, a::Taylor1{T}, r::Taylor1{T}, k::Int) where
diff --git a/src/identity.jl b/src/identity.jl
index 7197ce20..568e8a9a 100644
--- a/src/identity.jl
+++ b/src/identity.jl
@@ -17,16 +17,21 @@ for T in (:Taylor1, :TaylorN)
         identity!(c[k], a[k])
 end
 
+function identity!(c::Taylor1{T}, a::Taylor1{T}) where {T<:NumberNotSeries}
+    copyto!(c.coeffs, a.coeffs)
+    return nothing
+end
+
 function identity!(c::Taylor1{T}, a::Taylor1{T}, k::Int) where
         {T<:NumberNotSeries}
-    c[k] = a[k]
+    @inbounds c.coeffs[k+1] = a.coeffs[k+1]
     return nothing
 end
 
 @inline function identity!(c::Taylor1{T}, a::T, k::Int) where {T<:NumberNotSeries}
-    zero!(c[k])
+    @inbounds c.coeffs[k+1] = zero(c.coeffs[k+1])
     k != 0 && return nothing
-    c[0] = a
+    @inbounds c.coeffs[1] = a
     return nothing
 end
 function identity!(c::Taylor1{T}, a::T, k::Int) where {T<:Number}
diff --git a/src/power.jl b/src/power.jl
index 28a30364..376ee6f5 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -867,8 +867,10 @@ coefficient, which must be even.
 function sqrt!(c::Taylor1{T}, a::Taylor1{T}, ::Taylor1{T}, k::Int, k0::Int=0) where
         {T<:NumberNotSeries}
     k < k0 && return nothing
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
     if k == k0
-        @inbounds c[k] = sqrt(a[2*k0])
+        @inbounds c_coeffs[k+1] = sqrt(a_coeffs[2*k0+1])
         return nothing
     end
     # Recursion formula
@@ -877,17 +879,19 @@ function sqrt!(c::Taylor1{T}, a::Taylor1{T}, ::Taylor1{T}, k::Int, k0::Int=0) wh
     kend = (k - k0 - 2 + kodd) >> 1
     imax = min(k0+kend, order(a))
     imin = max(k0+1, k+k0-order(a))
+    kk = k+1
+    @inbounds acc = zero(c_coeffs[kk])
     if k+k0 ≤ order(a)
-        @inbounds c[k] = a[k+k0]
+        @inbounds acc = a_coeffs[k+k0+1]
     end
     if kodd == 0
-        @inbounds c[k] -= (c[kend+k0+1])^2
+        @inbounds acc -= (c_coeffs[kend+k0+2])^2
     end
-    imin ≤ imax && ( @inbounds c[k] -= 2 * c[imin] * c[k+k0-imin] )
+    imin ≤ imax && ( @inbounds acc -= 2 * c_coeffs[imin+1] * c_coeffs[k+k0-imin+1] )
     @inbounds for i = imin+1:imax
-        c[k] -= 2 * c[i] * c[k+k0-i]
+        acc -= 2 * c_coeffs[i+1] * c_coeffs[k+k0-i+1]
     end
-    @inbounds c[k] = c[k] / (2*c[k0])
+    @inbounds c_coeffs[kk] = acc / (2*c_coeffs[k0+1])
     return nothing
 end
 
diff --git a/src/zero_one.jl b/src/zero_one.jl
index ca33a1ab..f7b7ed83 100644
--- a/src/zero_one.jl
+++ b/src/zero_one.jl
@@ -55,7 +55,7 @@ ones(::Type{HomogeneousPolynomial{T}}, order::Int) where {T<:Number} =
 
 # Recursive functions (homogeneous coefficients)
 function zero!(a::Taylor1{T}, k::Int) where {T<:NumberNotSeries}
-    a[k] = zero(a[k])
+    @inbounds a.coeffs[k+1] = zero(a.coeffs[k+1])
     return nothing
 end
 
@@ -66,6 +66,14 @@ function zero!(a::Taylor1{T}, k::Int) where {T<:Number}
     return nothing
 end
 
+function zero!(a::Taylor1{T}) where {T<:NumberNotSeries}
+    a_coeffs = a.coeffs
+    @inbounds for i in eachindex(a_coeffs)
+        a_coeffs[i] = zero(a_coeffs[i])
+    end
+    return nothing
+end
+
 function zero!(a::Taylor1{T}) where {T<:Number}
     for k in eachindex(a)
         zero!(a, k)
@@ -116,7 +124,8 @@ end
 
 #
 function one!(a::Taylor1{T}, k::Int) where {T<:NumberNotSeries}
-    a[k] = k == 0 ? one(a[0]) : zero(a[k])
+    a_coeffs = a.coeffs
+    @inbounds a_coeffs[k+1] = k == 0 ? one(a_coeffs[1]) : zero(a_coeffs[k+1])
     return nothing
 end
 
@@ -129,6 +138,15 @@ function one!(a::Taylor1{T}, k::Int) where {T<:Number}
     return nothing
 end
 
+function one!(a::Taylor1{T}) where {T<:NumberNotSeries}
+    a_coeffs = a.coeffs
+    @inbounds for i in eachindex(a_coeffs)
+        a_coeffs[i] = zero(a_coeffs[i])
+    end
+    @inbounds a_coeffs[1] = one(a_coeffs[1])
+    return nothing
+end
+
 function one!(a::Taylor1{T}) where {T<:Number}
     for k in eachindex(a)
         one!(a, k)
@@ -147,8 +165,10 @@ end
 
 # Taylor1 (including nested Taylor1s)
 function one!(c::Taylor1{T}, a::Taylor1{T}, k::Int) where {T<:NumberNotSeries}
-    zero!(c, k)
-    (k == 0) && (@inbounds c[0] = one(constant_term(a)))
+    c_coeffs = c.coeffs
+    a_coeffs = a.coeffs
+    @inbounds c_coeffs[k+1] = zero(c_coeffs[k+1])
+    (k == 0) && (@inbounds c_coeffs[1] = one(a_coeffs[1]))
     return nothing
 end
 function one!(c::Taylor1{T}, a::Taylor1{T}, k::Int) where {T<:Number}

From d515ac5dba4b86a21fe48cde3c2f9cc4ab498820 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Sat, 6 Jun 2026 16:43:38 -0400
Subject: [PATCH 05/20] Fix #411

---
 src/TaylorSeries.jl    |  2 +-
 src/arithmetic.jl      | 15 ++++++++-------
 src/auxiliary.jl       | 23 +++++++++++++++++++++++
 test/mutatingfuncts.jl | 32 ++++++++++++++++++++++++++++++++
 4 files changed, 64 insertions(+), 8 deletions(-)

diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl
index 0155292b..80ec52af 100644
--- a/src/TaylorSeries.jl
+++ b/src/TaylorSeries.jl
@@ -56,7 +56,7 @@ export getcoeff, derivative, integrate, differentiate,
     get_variable_names, get_variable_symbols,
     # jacobian, hessian, jacobian!, hessian!,
     ∇, taylor_expand, update!,
-    constant_term, linear_polynomial, nonlinear_polynomial,
+    constant_term, constant_term!, linear_polynomial, nonlinear_polynomial,
     normalize_taylor, norm
 
 const TS = TaylorSeries
diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index d08f3cf8..b36753c6 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -1824,13 +1824,14 @@ function div!(c::TaylorN, a::NumberNotSeries, b::TaylorN)
     return nothing
 end
 
-# # c[k] <- a[k]/b, where b is a scalar
-# function div!(c::TaylorN, a::TaylorN, b::NumberNotSeries)
-#     @inbounds for k in eachindex(c)
-#         div!(c[k], a[k], b)
-#     end
-#     return nothing
-# end
+# c[k] <- a[k]/b, where b is a scalar
+function div!(c::TaylorN{T}, a::TaylorN{T}, b::NumberNotSeries) where {T<:Number}
+    _check_same_space(c, a)
+    @inbounds for k in eachindex(c)
+        div!(c[k], a[k], b)
+    end
+    return nothing
+end
 
 # NOTE: Here `div!` *accumulates* the result of a / b in res[k] (k > 0)
 @inline function div!(c::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
diff --git a/src/auxiliary.jl b/src/auxiliary.jl
index 7c311de1..12167923 100644
--- a/src/auxiliary.jl
+++ b/src/auxiliary.jl
@@ -514,6 +514,29 @@ constant_term(a::Vector{T}) where {T<:Number} = constant_term.(a)
 
 constant_term(a::Number) = a
 
+"""
+    constant_term!(a, c)
+
+Update the constant term (zero order coefficient) of `a` to `c`,
+leaving higher order coefficients unchanged.
+"""
+function constant_term!(a::Taylor1{T}, c::T) where {T<:Number}
+    a[0] = c
+    return a
+end
+
+function constant_term!(a::TaylorN{T}, c::T) where {T<:Number}
+    a[0][1] = c
+    return a
+end
+
+function constant_term!(a::HomogeneousPolynomial{T}, c::T) where {T<:Number}
+    iszero(order(a)) ||
+        throw(ArgumentError("only zero-order HomogeneousPolynomial has a constant term"))
+    a[1] = c
+    return a
+end
+
 """
     linear_polynomial(a)
 
diff --git a/test/mutatingfuncts.jl b/test/mutatingfuncts.jl
index 6c6eeb56..a770594e 100644
--- a/test/mutatingfuncts.jl
+++ b/test/mutatingfuncts.jl
@@ -125,4 +125,36 @@ using Test
     TS.mul!(ninplace, nt2)
     @test ninplace ≈ nt1 * nt2
 
+    ct1 = 1 + t1 + t1^2
+    ct1_tail = collect(ct1[1:order(ct1)])
+    @test constant_term!(ct1, 3.5) === ct1
+    @test constant_term(ct1) == 3.5
+    @test collect(ct1[1:order(ct1)]) == ct1_tail
+
+    nct1 = deepcopy(nt1)
+    nct1_tail = [nct1[k] for k in 1:order(nct1)]
+    new_inner_constant = Taylor1([2.0, -1.0, 0.25], 2)
+    @test constant_term!(nct1, new_inner_constant) === nct1
+    @test constant_term(nct1) == new_inner_constant
+    @test all(nct1[k] == nct1_tail[k] for k in 1:order(nct1))
+
+    space = JetSpace(order=3, variables=[:x, :y])
+    x, y = variables(space)
+    tn = 1.0 + x + 2y + 3x*y
+    tn_tail = [tn[k] for k in 1:order(tn)]
+    @test constant_term!(tn, 7.5) === tn
+    @test constant_term(tn) == 7.5
+    @test all(tn[k] == tn_tail[k] for k in 1:order(tn))
+
+    hp0 = HomogeneousPolynomial(space, 1.0, 0)
+    @test constant_term!(hp0, 4.0) === hp0
+    @test hp0[1] == 4.0
+    hp1 = HomogeneousPolynomial(space, [1.0, 0.0], 1)
+    @test_throws ArgumentError constant_term!(hp1, 2.0)
+
+    tn_div = 2.0 + x + 4y + 6x*y
+    tn_res = zero(tn_div)
+    @test isnothing(TS.div!(tn_res, tn_div, 2.0))
+    @test tn_res == tn_div / 2.0
+
 end

From 7e546d2353e312477a6a88966e43191548f1213a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Sat, 6 Jun 2026 16:55:00 -0400
Subject: [PATCH 06/20] Fix docs build

---
 docs/src/api.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/src/api.md b/docs/src/api.md
index 85f4968b..67d2f03b 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -106,6 +106,7 @@ asinh!
 acosh!
 atanh!
 differentiate!
+constant_term!
 _isthinzero
 _internalmutfunc_call
 _dict_unary_ops

From 63f90ce5d4c5caa7fe5126ccc946411625d905c7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Sun, 7 Jun 2026 18:17:41 -0400
Subject: [PATCH 07/20] Try a similar approach for Taylor1{TaylorN{Float64}}

---
 src/arithmetic.jl      | 127 +++++++++++++++++++++++++++++++++--------
 src/calculus.jl        |   8 ++-
 src/functions.jl       | 107 +++++++++++++++++++++++-----------
 test/mutatingfuncts.jl |  46 +++++++++++++++
 4 files changed, 230 insertions(+), 58 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index b36753c6..9a7e682e 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -361,6 +361,50 @@ function subst!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
     return nothing
 end
 
+function add!(v::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
+        b::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
+    v_coeffs = v.coeffs
+    a_coeffs = a.coeffs
+    b_coeffs = b.coeffs
+    _check_same_space(v_coeffs[1], a_coeffs[1], b_coeffs[1])
+    @inbounds for i in eachindex(v_coeffs)
+        v_hps = v_coeffs[i].coeffs
+        a_hps = a_coeffs[i].coeffs
+        b_hps = b_coeffs[i].coeffs
+        for j in eachindex(v_hps)
+            v_hp = v_hps[j].coeffs
+            a_hp = a_hps[j].coeffs
+            b_hp = b_hps[j].coeffs
+            for k in eachindex(v_hp)
+                v_hp[k] = a_hp[k] + b_hp[k]
+            end
+        end
+    end
+    return nothing
+end
+
+function subst!(v::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
+        b::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
+    v_coeffs = v.coeffs
+    a_coeffs = a.coeffs
+    b_coeffs = b.coeffs
+    _check_same_space(v_coeffs[1], a_coeffs[1], b_coeffs[1])
+    @inbounds for i in eachindex(v_coeffs)
+        v_hps = v_coeffs[i].coeffs
+        a_hps = a_coeffs[i].coeffs
+        b_hps = b_coeffs[i].coeffs
+        for j in eachindex(v_hps)
+            v_hp = v_hps[j].coeffs
+            a_hp = a_hps[j].coeffs
+            b_hp = b_hps[j].coeffs
+            for k in eachindex(v_hp)
+                v_hp[k] = a_hp[k] - b_hp[k]
+            end
+        end
+    end
+    return nothing
+end
+
 function +(a::Taylor1{T}, b::Taylor1{T}) where {T<:NumberNotSeries}
     if order(a) != order(b)
         a, b = fixorder(a, b)
@@ -409,45 +453,70 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
                 a, b = fixorder(a, b)
             end
             c = zero(a)
-            for k in eachindex(a)
-                ($fc)(c, a, b, k)
-            end
+            ($fc)(c, a, b)
             return c
         end
         function ($fc)(v::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
                 b::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries}
-            @inbounds for i in eachindex(v[k])
-                for j in eachindex(v[k][i])
-                    v[k][i][j] = ($f)(a[k][i][j], b[k][i][j])
+            v_k = v.coeffs[k+1]
+            a_k = a.coeffs[k+1]
+            b_k = b.coeffs[k+1]
+            v_hps = v_k.coeffs
+            a_hps = a_k.coeffs
+            b_hps = b_k.coeffs
+            @inbounds for i in eachindex(v_hps)
+                v_hp = v_hps[i].coeffs
+                a_hp = a_hps[i].coeffs
+                b_hp = b_hps[i].coeffs
+                for j in eachindex(v_hp)
+                    v_hp[j] = ($f)(a_hp[j], b_hp[j])
                 end
             end
             return nothing
         end
         function ($fc)(v::Taylor1{TaylorN{T}}, a::NumberNotSeries,
                 b::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries}
-            @inbounds for i in eachindex(v[k])
-                aaa = ifelse(k == 0 && i == 0, a, zero(a))
-                for j in eachindex(v[k][i])
-                    v[k][i][j] = ($f)(aaa, b[k][i][j])
+            v_k = v.coeffs[k+1]
+            b_k = b.coeffs[k+1]
+            v_hps = v_k.coeffs
+            b_hps = b_k.coeffs
+            @inbounds for i in eachindex(v_hps)
+                aaa = ifelse(k == 0 && i == 1, a, zero(a))
+                v_hp = v_hps[i].coeffs
+                b_hp = b_hps[i].coeffs
+                for j in eachindex(v_hp)
+                    v_hp[j] = ($f)(aaa, b_hp[j])
                 end
             end
             return nothing
         end
         function ($fc)(v::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}},
                 a::NumberNotSeries, k::Int) where {T<:NumberNotSeries}
-            @inbounds for i in eachindex(v[k])
-                aaa = ifelse(k == 0 && i == 0, a, zero(a))
-                for j in eachindex(v[k][i])
-                    v[k][i][j] = ($f)(b[k][i][j], aaa)
+            v_k = v.coeffs[k+1]
+            b_k = b.coeffs[k+1]
+            v_hps = v_k.coeffs
+            b_hps = b_k.coeffs
+            @inbounds for i in eachindex(v_hps)
+                aaa = ifelse(k == 0 && i == 1, a, zero(a))
+                v_hp = v_hps[i].coeffs
+                b_hp = b_hps[i].coeffs
+                for j in eachindex(v_hp)
+                    v_hp[j] = ($f)(b_hp[j], aaa)
                 end
             end
             return nothing
         end
         function ($fc)(v::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
                 k::Int) where {T<:NumberNotSeries}
-            @inbounds for l in eachindex(v[k])
-                for m in eachindex(v[k][l])
-                    v[k][l][m] = ($f)(a[k][l][m])
+            v_k = v.coeffs[k+1]
+            a_k = a.coeffs[k+1]
+            v_hps = v_k.coeffs
+            a_hps = a_k.coeffs
+            @inbounds for l in eachindex(v_hps)
+                v_hp = v_hps[l].coeffs
+                a_hp = a_hps[l].coeffs
+                for m in eachindex(v_hp)
+                    v_hp[m] = ($f)(a_hp[m])
                 end
             end
             return nothing
@@ -1010,9 +1079,15 @@ end
 
 function mul!(res::Taylor1{TaylorN{T}}, a::NumberNotSeries,
         b::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries}
-    for l in eachindex(b[k])
-        for m in eachindex(b[k][l])
-            res[k][l][m] = a*b[k][l][m]
+    res_k = res.coeffs[k+1]
+    b_k = b.coeffs[k+1]
+    res_hps = res_k.coeffs
+    b_hps = b_k.coeffs
+    @inbounds for l in eachindex(res_hps)
+        res_hp = res_hps[l].coeffs
+        b_hp = b_hps[l].coeffs
+        for m in eachindex(res_hp)
+            res_hp[m] = a*b_hp[m]
         end
     end
     return nothing
@@ -1878,9 +1953,15 @@ end
 
 @inline function div!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
         b::NumberNotSeries, k::Int) where {T<:NumberNotSeries}
-    for l in eachindex(a[k])
-        for m in eachindex(a[k][l])
-            res[k][l][m] = a[k][l][m]/b
+    res_k = res.coeffs[k+1]
+    a_k = a.coeffs[k+1]
+    res_hps = res_k.coeffs
+    a_hps = a_k.coeffs
+    @inbounds for l in eachindex(res_hps)
+        res_hp = res_hps[l].coeffs
+        a_hp = a_hps[l].coeffs
+        for m in eachindex(res_hp)
+            res_hp[m] = a_hp[m]/b
         end
     end
     return nothing
diff --git a/src/calculus.jl b/src/calculus.jl
index 728bab68..a639fa41 100644
--- a/src/calculus.jl
+++ b/src/calculus.jl
@@ -73,9 +73,11 @@ end
 function differentiate!(p::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, k::Int) where
         {T<:NumberNotSeries}
     k >= order(a) && return nothing
-    @inbounds for ord in eachindex(p[k])
-        zero!(p[k], ord)
-        mul!(p[k], a[k+1], k+1, ord)
+    p_k = p.coeffs[k+1]
+    a_next = a.coeffs[k+2]
+    @inbounds for ord in eachindex(p_k)
+        zero!(p_k, ord)
+        mul!(p_k, a_next, k+1, ord)
     end
     return nothing
 end
diff --git a/src/functions.jl b/src/functions.jl
index 8f6ee70b..089e433f 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -232,6 +232,18 @@ for T in (:Taylor1, :TaylorN)
     end
 end
 
+function log1p(a::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
+    aux = log1p(constant_term(a))
+    aa = convert(Taylor1{typeof(aux)}, a)
+    c = _constant_series_like(a, aux, order(a))
+    denom = _constant_series_like(aux, zero(aux[0][1]), order(aux))
+    _log1p_denominator!(denom, aa.coeffs[1])
+    for k in eachindex(a)
+        log1p!(c, aa, k, denom)
+    end
+    return c
+end
+
 # Taylor1 mutating functions
 function exp!(c::Taylor1{T}, a::Taylor1{T}, k::Int) where
         {T<:NumberNotSeries}
@@ -1366,72 +1378,103 @@ end
 
 function log!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
         k::Int) where {T<:NumberNotSeries}
+    res_coeffs = res.coeffs
+    a_coeffs = a.coeffs
+    res_k = res_coeffs[k+1]
+    a0 = a_coeffs[1]
     if k == 0
-        @inbounds for ordQ in eachindex(a[0])
+        @inbounds for ordQ in eachindex(a0)
             # zero!(res[k], a[0], ordQ)
-            log!(res[k], a[0], ordQ)
+            log!(res_k, a0, ordQ)
         end
         return nothing
     elseif k == 1
-        @inbounds for ordQ in eachindex(a[0])
-            zero!(res[k][ordQ])
-            div!(res[k], a[1], a[0], ordQ)
+        a1 = a_coeffs[2]
+        @inbounds for ordQ in eachindex(a0)
+            zero!(res_k[ordQ])
+            div!(res_k, a1, a0, ordQ)
         end
         return nothing
     end
     # The recursion formula
-    tmp = _constant_series_like(a[k], zero(a[k][0][1]), order(a[0]))
-    zero!(res[k])
+    zero!(res_k)
     for i = 1:k-1
-        @inbounds for ordQ in eachindex(a[0])
-            tmp[ordQ] = (k-i) * res[k-i][ordQ]
-            mul!(res[k], tmp, a[i], ordQ)
+        res_ki = res_coeffs[k-i+1]
+        a_i = a_coeffs[i+1]
+        @inbounds for ordQ in eachindex(a0)
+            mul_scalar!(res_k, k-i, res_ki, a_i, ordQ)
         end
     end
     div!(res, res, k, k)
-    @inbounds for ordQ in eachindex(a[0])
-        subst!(tmp, a[k], res[k], ordQ)
-        zero!(res[k][ordQ])
-        div!(res[k], tmp, a[0], ordQ)
+    a_k = a_coeffs[k+1]
+    @inbounds for ordQ in eachindex(a0)
+        subst!(res_k, a_k, res_k, ordQ)
+        div!(res_k, a0, ordQ)
     end
     return nothing
 end
 
 function log1p!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
         k::Int) where {T<:NumberNotSeries}
+    res_coeffs = res.coeffs
+    a_coeffs = a.coeffs
+    res_k = res_coeffs[k+1]
+    a0 = a_coeffs[1]
     if k == 0
-        @inbounds for ordQ in eachindex(a[0])
+        @inbounds for ordQ in eachindex(a0)
             # zero!(res[k], a[0], ordQ)
-            log1p!(res[k], a[0], ordQ)
+            log1p!(res_k, a0, ordQ)
         end
         return nothing
     end
-    tmp1 = _constant_series_like(a[k], zero(a[k][0][1]), order(a[0]))
-    zero!(res[k])
-    @inbounds for ordQ in eachindex(a[0])
-        # zero!(res[k], a[0], ordQ)
-        one!(tmp1, a[0], ordQ)
-        add!(tmp1, tmp1, a[0], ordQ)
+    tmp1 = _constant_series_like(a_coeffs[k+1], zero(a_coeffs[k+1][0][1]), order(a0))
+    _log1p_denominator!(tmp1, a0)
+    log1p!(res, a, k, tmp1)
+    return nothing
+end
+
+function _log1p_denominator!(tmp::TaylorN{T}, a0::TaylorN{T}) where {T<:Number}
+    zero!(tmp)
+    @inbounds for ordQ in eachindex(a0)
+        one!(tmp, a0, ordQ)
+        add!(tmp, tmp, a0, ordQ)
     end
+    return tmp
+end
+
+function log1p!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
+        k::Int, tmp1::TaylorN{T}) where {T<:NumberNotSeries}
+    res_coeffs = res.coeffs
+    a_coeffs = a.coeffs
+    res_k = res_coeffs[k+1]
+    a0 = a_coeffs[1]
+    if k == 0
+        @inbounds for ordQ in eachindex(a0)
+            log1p!(res_k, a0, ordQ)
+        end
+        return nothing
+    end
+    zero!(res_k)
     if k == 1
-        @inbounds for ordQ in eachindex(a[0])
-            div!(res[k], a[1], tmp1, ordQ)
+        a1 = a_coeffs[2]
+        @inbounds for ordQ in eachindex(a0)
+            div!(res_k, a1, tmp1, ordQ)
         end
         return nothing
     end
     # The recursion formula
-    tmp = _constant_series_like(a[k], zero(a[k][0][1]), order(a[0]))
     for i = 1:k-1
-        @inbounds for ordQ in eachindex(a[0])
-            tmp[ordQ] = (k-i) * res[k-i][ordQ]
-            mul!(res[k], tmp, a[i], ordQ)
+        res_ki = res_coeffs[k-i+1]
+        a_i = a_coeffs[i+1]
+        @inbounds for ordQ in eachindex(a0)
+            mul_scalar!(res_k, k-i, res_ki, a_i, ordQ)
         end
     end
     div!(res, res, k, k)
-    @inbounds for ordQ in eachindex(a[0])
-        subst!(tmp, a[k], res[k], ordQ)
-        zero!(res[k][ordQ])
-        div!(res[k], tmp, tmp1, ordQ)
+    a_k = a_coeffs[k+1]
+    @inbounds for ordQ in eachindex(a0)
+        subst!(res_k, a_k, res_k, ordQ)
+        div!(res_k, tmp1, ordQ)
     end
     return nothing
 end
diff --git a/test/mutatingfuncts.jl b/test/mutatingfuncts.jl
index a770594e..8c991d69 100644
--- a/test/mutatingfuncts.jl
+++ b/test/mutatingfuncts.jl
@@ -157,4 +157,50 @@ using Test
     @test isnothing(TS.div!(tn_res, tn_div, 2.0))
     @test tn_res == tn_div / 2.0
 
+    t1tn1 = Taylor1([
+        (2.0 + 0.1i) + (0.2 + 0.01i)*x + (0.1 - 0.02i)*y +
+            (0.03 + 0.01i)*x*y for i in 0:4], 4)
+    t1tn2 = Taylor1([
+        (1.0 - 0.05i) + (0.3 - 0.02i)*x - (0.04 + 0.01i)*y +
+            (0.02 - 0.005i)*x^2 for i in 0:4], 4)
+    t1tn_res = zero(t1tn1)
+    TS.add!(t1tn_res, t1tn1, t1tn2)
+    @test t1tn_res == t1tn1 + t1tn2
+    TS.subst!(t1tn_res, t1tn1, t1tn2)
+    @test t1tn_res == t1tn1 - t1tn2
+    for k in eachindex(t1tn_res)
+        TS.mul!(t1tn_res, 2.5, t1tn1, k)
+    end
+    @test t1tn_res == 2.5 * t1tn1
+    for k in eachindex(t1tn_res)
+        TS.div!(t1tn_res, t1tn1, 2.5, k)
+    end
+    t1tn_div_expected = zero(t1tn1)
+    for k in eachindex(t1tn_div_expected)
+        TS.mul!(t1tn_div_expected, 0.4, t1tn1, k)
+    end
+    @test t1tn_res ≈ t1tn_div_expected
+    d_t1tn = Taylor1(zero(t1tn1[0]), order(t1tn1)-1)
+    TS.differentiate!(d_t1tn, t1tn1)
+    @test d_t1tn == differentiate(t1tn1)
+    log_t1tn = zero(t1tn1)
+    for k in eachindex(log_t1tn)
+        TS.log!(log_t1tn, t1tn1, k)
+    end
+    @test exp(log_t1tn) ≈ t1tn1
+    small_t1tn = zero(t1tn1)
+    for k in eachindex(small_t1tn)
+        TS.div!(small_t1tn, t1tn1, 20.0, k)
+    end
+    log1p_t1tn = zero(small_t1tn)
+    for k in eachindex(log1p_t1tn)
+        TS.log1p!(log1p_t1tn, small_t1tn, k)
+    end
+    one_plus_small = zero(small_t1tn)
+    for k in eachindex(one_plus_small)
+        TS.add!(one_plus_small, 1.0, small_t1tn, k)
+    end
+    @test exp(log1p_t1tn) ≈ one_plus_small
+    @test log1p(small_t1tn) ≈ log1p_t1tn
+
 end

From 21bc23460ced68bd93c8c3b77e86301da8a5b9d9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20P=C3=A9rez?= <PerezHz@users.noreply.github.com>
Date: Tue, 9 Jun 2026 22:08:39 -0400
Subject: [PATCH 08/20] Apply suggestion from @lbenet

Co-authored-by: Luis Benet <lbenet@users.noreply.github.com>
---
 src/arithmetic.jl | 21 +++++++++++++++++++++
 1 file changed, 21 insertions(+)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 9a7e682e..9acf7bdb 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -445,6 +445,27 @@ end
 
 for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
     @eval begin
+        function ($fc)(v::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
+                b::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
+            v_coeffs = v.coeffs
+            a_coeffs = a.coeffs
+            b_coeffs = b.coeffs
+            _check_same_space(v_coeffs[1], a_coeffs[1], b_coeffs[1])
+            @inbounds for i in eachindex(v_coeffs)
+                v_hps = v_coeffs[i].coeffs
+                a_hps = a_coeffs[i].coeffs
+                b_hps = b_coeffs[i].coeffs
+                for j in eachindex(v_hps)
+                    v_hp = v_hps[j].coeffs
+                    a_hp = a_hps[j].coeffs
+                    b_hp = b_hps[j].coeffs
+                    for k in eachindex(v_hp)
+                        v_hp[k] = ($f)(a_hp[k], b_hp[k])
+                    end
+                end
+            end
+            return nothing
+        end
         function ($f)(a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}}) where
                 {T<:NumberNotSeries}
             _check_same_space(a[0], b[0])

From bc7c4ea971927a87e23de1fb47d4df627f1f36f0 Mon Sep 17 00:00:00 2001
From: Jorge Perez-Hernandez <jorge.perez_hernandez@cfa.harvard.edu>
Date: Tue, 9 Jun 2026 22:10:05 -0400
Subject: [PATCH 09/20] Remove methods (follow-up to last commit)

---
 src/arithmetic.jl | 44 --------------------------------------------
 1 file changed, 44 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 9acf7bdb..6f3422e0 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -361,50 +361,6 @@ function subst!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
     return nothing
 end
 
-function add!(v::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
-        b::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
-    v_coeffs = v.coeffs
-    a_coeffs = a.coeffs
-    b_coeffs = b.coeffs
-    _check_same_space(v_coeffs[1], a_coeffs[1], b_coeffs[1])
-    @inbounds for i in eachindex(v_coeffs)
-        v_hps = v_coeffs[i].coeffs
-        a_hps = a_coeffs[i].coeffs
-        b_hps = b_coeffs[i].coeffs
-        for j in eachindex(v_hps)
-            v_hp = v_hps[j].coeffs
-            a_hp = a_hps[j].coeffs
-            b_hp = b_hps[j].coeffs
-            for k in eachindex(v_hp)
-                v_hp[k] = a_hp[k] + b_hp[k]
-            end
-        end
-    end
-    return nothing
-end
-
-function subst!(v::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
-        b::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
-    v_coeffs = v.coeffs
-    a_coeffs = a.coeffs
-    b_coeffs = b.coeffs
-    _check_same_space(v_coeffs[1], a_coeffs[1], b_coeffs[1])
-    @inbounds for i in eachindex(v_coeffs)
-        v_hps = v_coeffs[i].coeffs
-        a_hps = a_coeffs[i].coeffs
-        b_hps = b_coeffs[i].coeffs
-        for j in eachindex(v_hps)
-            v_hp = v_hps[j].coeffs
-            a_hp = a_hps[j].coeffs
-            b_hp = b_hps[j].coeffs
-            for k in eachindex(v_hp)
-                v_hp[k] = a_hp[k] - b_hp[k]
-            end
-        end
-    end
-    return nothing
-end
-
 function +(a::Taylor1{T}, b::Taylor1{T}) where {T<:NumberNotSeries}
     if order(a) != order(b)
         a, b = fixorder(a, b)

From 1593d3f2c020cf7b835f9c0753f12a341d865da0 Mon Sep 17 00:00:00 2001
From: Jorge Perez-Hernandez <jorge.perez_hernandez@cfa.harvard.edu>
Date: Wed, 10 Jun 2026 22:24:37 -0400
Subject: [PATCH 10/20] Move add! subst! methods

---
 src/arithmetic.jl | 24 ++++++++++++------------
 1 file changed, 12 insertions(+), 12 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 6f3422e0..c756aaea 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -361,6 +361,18 @@ function subst!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
     return nothing
 end
 
+function add!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
+        b::Taylor1{Taylor1{T}}, k::Int) where {T<:NumberNotSeries}
+    @inbounds add!(v.coeffs[k+1], a.coeffs[k+1], b.coeffs[k+1])
+    return nothing
+end
+
+function subst!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
+        b::Taylor1{Taylor1{T}}, k::Int) where {T<:NumberNotSeries}
+    @inbounds subst!(v.coeffs[k+1], a.coeffs[k+1], b.coeffs[k+1])
+    return nothing
+end
+
 function +(a::Taylor1{T}, b::Taylor1{T}) where {T<:NumberNotSeries}
     if order(a) != order(b)
         a, b = fixorder(a, b)
@@ -793,18 +805,6 @@ function mul_scalar!(c::TaylorN{T}, scalar::NumberNotSeries, a::TaylorN{T},
 end
 
 # Nested Taylor1s
-function add!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
-        b::Taylor1{Taylor1{T}}, k::Int) where {T<:NumberNotSeries}
-    @inbounds add!(v.coeffs[k+1], a.coeffs[k+1], b.coeffs[k+1])
-    return nothing
-end
-
-function subst!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
-        b::Taylor1{Taylor1{T}}, k::Int) where {T<:NumberNotSeries}
-    @inbounds subst!(v.coeffs[k+1], a.coeffs[k+1], b.coeffs[k+1])
-    return nothing
-end
-
 function mul!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, b::Taylor1{Taylor1{T}},
         k::Int) where {T<:NumberNotSeries}
     c_coeffs = c.coeffs

From 15abfd83ac6cc967c49fee10632ca973997818cc Mon Sep 17 00:00:00 2001
From: Jorge Perez-Hernandez <jorge.perez_hernandez@cfa.harvard.edu>
Date: Thu, 11 Jun 2026 00:03:13 -0400
Subject: [PATCH 11/20] Add missing method (fix #411)

---
 src/arithmetic.jl      | 13 +++++++++++--
 test/mutatingfuncts.jl |  5 +++++
 2 files changed, 16 insertions(+), 2 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index c756aaea..fb9fd90a 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -1836,6 +1836,16 @@ end
     return nothing
 end
 
+# c[k] <- a[k]/b, where b is a scalar
+@inline function div!(c::TaylorN{T}, a::TaylorN{T}, b::NumberNotSeries,
+        k::Int) where {T<:Number}
+    _check_same_space(c, a)
+    @inbounds for i in eachindex(c[k])
+        c[k][i] = a[k][i] / b
+    end
+    return nothing
+end
+
 # in-place division c <- c/a (assumes equal order among TaylorNs)
 function div!(c::TaylorN, a::TaylorN)
     @inbounds for k in eachindex(c)
@@ -1878,9 +1888,8 @@ end
 
 # c[k] <- a[k]/b, where b is a scalar
 function div!(c::TaylorN{T}, a::TaylorN{T}, b::NumberNotSeries) where {T<:Number}
-    _check_same_space(c, a)
     @inbounds for k in eachindex(c)
-        div!(c[k], a[k], b)
+        div!(c, a, b, k)
     end
     return nothing
 end
diff --git a/test/mutatingfuncts.jl b/test/mutatingfuncts.jl
index 8c991d69..81ac2bc2 100644
--- a/test/mutatingfuncts.jl
+++ b/test/mutatingfuncts.jl
@@ -156,6 +156,11 @@ using Test
     tn_res = zero(tn_div)
     @test isnothing(TS.div!(tn_res, tn_div, 2.0))
     @test tn_res == tn_div / 2.0
+    tn_res_k = zero(tn_div)
+    for k in eachindex(tn_res_k)
+        @test isnothing(TS.div!(tn_res_k, tn_div, 2.0, k))
+    end
+    @test tn_res_k == tn_div / 2.0
 
     t1tn1 = Taylor1([
         (2.0 + 0.1i) + (0.2 + 0.01i)*x + (0.1 - 0.02i)*y +

From 2fe3eb0a16565f034edc58718ffb8e8fa7197711 Mon Sep 17 00:00:00 2001
From: Jorge Perez-Hernandez <jorge.perez_hernandez@cfa.harvard.edu>
Date: Thu, 11 Jun 2026 12:01:54 -0400
Subject: [PATCH 12/20] Add test to catch warning from variables! with
 nowarn=false

---
 test/manyvariables.jl | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

diff --git a/test/manyvariables.jl b/test/manyvariables.jl
index b9d16085..8dcd2d17 100644
--- a/test/manyvariables.jl
+++ b/test/manyvariables.jl
@@ -52,6 +52,23 @@ using Test
     @test sp.index_table[2][1] == 7
     @test TS.in_base(order(), [2,1]) == 15
     @test sp.pos_table[4][15] == 2
+
+    default_space_warning = r"Updating TaylorSeries.default_space\[\]"
+    @test_logs (:warn, default_space_warning) variables!("w", order=2,
+        numvars=2)
+    @test_logs (:warn, default_space_warning) variables!(:w, order=2,
+        numvars=2)
+    @test_logs (:warn, default_space_warning) variables!(["w", "z"],
+        order=2)
+    @test_logs (:warn, default_space_warning) variables!([:w, :z], order=2)
+    @test_logs (:warn, default_space_warning) variables!(Int, "w", order=2,
+        numvars=2)
+    @test_logs (:warn, default_space_warning) variables!(Int, :w, order=2,
+        numvars=2)
+    @test_logs (:warn, default_space_warning) variables!(BigInt, ["w", "z"],
+        order=2)
+    @test_logs (:warn, default_space_warning) variables!(BigInt, [:w, :z],
+        order=2)
 end
 
 @testset "Explicit JetSpaces" begin

From 791c6602fde95ded8f11d14c0b98e1f2d91553db Mon Sep 17 00:00:00 2001
From: Jorge Perez-Hernandez <jorge.perez_hernandez@cfa.harvard.edu>
Date: Thu, 11 Jun 2026 12:09:30 -0400
Subject: [PATCH 13/20] Rename _pow_taylor1_cached! -> _pow_cached!

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

diff --git a/src/power.jl b/src/power.jl
index 376ee6f5..f3c4960e 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -93,7 +93,7 @@ function _pow(a::Taylor1{T}, r::S) where {T<:NumberNotSeries, S<:Real}
     order_a = order(a)
     lastnz = isinteger(r) && r > 0 ? findlast(a) : order_a
     for k in eachindex(c)
-        _pow_taylor1_cached!(c, a, r, k, l0, lnull, lastnz, order_a)
+        _pow_cached!(c, a, r, k, l0, lnull, lastnz, order_a)
     end
     return c
 end
@@ -231,7 +231,7 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`.
 
 """ pow!
 
-@inline function _pow_taylor1_cached!(c::Taylor1{T}, a::Taylor1{T},
+@inline function _pow_cached!(c::Taylor1{T}, a::Taylor1{T},
         r::S, k::Int, l0::Int, lnull::Int, lastnz::Int,
         order_a::Int) where {T<:NumberNotSeries, S<:Real}
     zero!(c, k)
@@ -279,7 +279,7 @@ function pow!(c::Taylor1{T}, a::Taylor1{T}, aux::Taylor1{T},
     lnull = trunc(Int, r*l0)
     order_a = order(a)
     lastnz = isinteger(r) && r > 0 ? findlast(a) : order_a
-    return _pow_taylor1_cached!(c, a, r, k, l0, lnull, lastnz, order_a)
+    return _pow_cached!(c, a, r, k, l0, lnull, lastnz, order_a)
 end
 
 function pow!(c::TaylorN{T}, a::TaylorN{T}, aux::TaylorN{T},

From d26e622e01d65d669f737bff906db420bbc350e1 Mon Sep 17 00:00:00 2001
From: Jorge Perez-Hernandez <jorge.perez_hernandez@cfa.harvard.edu>
Date: Thu, 11 Jun 2026 12:32:19 -0400
Subject: [PATCH 14/20] Inline identity! for scalar coeffs case, and switch
 from copyto! to identity! in call site

---
 src/identity.jl | 2 +-
 src/power.jl    | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/identity.jl b/src/identity.jl
index 568e8a9a..1c1cceb3 100644
--- a/src/identity.jl
+++ b/src/identity.jl
@@ -17,7 +17,7 @@ for T in (:Taylor1, :TaylorN)
         identity!(c[k], a[k])
 end
 
-function identity!(c::Taylor1{T}, a::Taylor1{T}) where {T<:NumberNotSeries}
+@inline function identity!(c::Taylor1{T}, a::Taylor1{T}) where {T<:NumberNotSeries}
     copyto!(c.coeffs, a.coeffs)
     return nothing
 end
diff --git a/src/power.jl b/src/power.jl
index f3c4960e..7162ed75 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -407,7 +407,7 @@ function pow!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
         add!(c_k, c_k, aux_k)
     end
     # c[k] = c[k] / (kprime * a[l0])
-    @inbounds copyto!(aux_k.coeffs, c_k.coeffs)
+    identity!(aux_k, c_k)
     @inbounds for j in eachindex(a_l0)
         div!(c_k, aux_k, a_l0, j)
     end

From 9b6d3271db8253fb49099c8b59cf0b739461aada Mon Sep 17 00:00:00 2001
From: Jorge Perez-Hernandez <jorge.perez_hernandez@cfa.harvard.edu>
Date: Thu, 11 Jun 2026 13:44:13 -0400
Subject: [PATCH 15/20] Apply suggestion by @lbenet

---
 src/evaluate.jl | 9 +--------
 1 file changed, 1 insertion(+), 8 deletions(-)

diff --git a/src/evaluate.jl b/src/evaluate.jl
index 44676e80..dd2f60fd 100644
--- a/src/evaluate.jl
+++ b/src/evaluate.jl
@@ -65,14 +65,7 @@ function evaluate(a::Taylor1{T}, x::Taylor1{T}) where {T<:NumberNotSeries}
     end
     suma = zero(x)
     aux = zero(x)
-    suma_coeffs = suma.coeffs
-    aux_coeffs = aux.coeffs
-    a_coeffs = a.coeffs
-    @inbounds for k in reverse(eachindex(a_coeffs))
-        mul!(aux, suma, x)
-        copyto!(suma_coeffs, aux_coeffs)
-        suma_coeffs[1] += a_coeffs[k]
-    end
+    _horner!(suma, a, x, aux)
     return suma
 end
 

From f2794f861dbd1cadd5839b974b7b540b3014b196 Mon Sep 17 00:00:00 2001
From: Jorge Perez-Hernandez <jorge.perez_hernandez@cfa.harvard.edu>
Date: Thu, 11 Jun 2026 15:45:59 -0400
Subject: [PATCH 16/20] Apply suggeston by @lbenet

---
 src/arithmetic.jl | 32 +++++++++++++++++++++++---------
 1 file changed, 23 insertions(+), 9 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index fb9fd90a..bc4be4a3 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -317,24 +317,38 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
     end
 end
 
-function add!(v::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}) where
-        {T<:NumberNotSeries}
+@inline function add!(v::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T},
+        k::Int) where {T<:NumberNotSeries}
     v_coeffs = v.coeffs
     a_coeffs = a.coeffs
     b_coeffs = b.coeffs
-    @inbounds for i in eachindex(v_coeffs)
-        v_coeffs[i] = a_coeffs[i] + b_coeffs[i]
-    end
+    kk = k + 1
+    @inbounds v_coeffs[kk] = a_coeffs[kk] + b_coeffs[kk]
     return nothing
 end
 
-function subst!(v::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}) where
-        {T<:NumberNotSeries}
+@inline function subst!(v::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T},
+        k::Int) where {T<:NumberNotSeries}
     v_coeffs = v.coeffs
     a_coeffs = a.coeffs
     b_coeffs = b.coeffs
-    @inbounds for i in eachindex(v_coeffs)
-        v_coeffs[i] = a_coeffs[i] - b_coeffs[i]
+    kk = k + 1
+    @inbounds v_coeffs[kk] = a_coeffs[kk] - b_coeffs[kk]
+    return nothing
+end
+
+function add!(v::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}) where
+        {T<:NumberNotSeries}
+    for k in eachindex(v)
+        add!(v, a, b, k)
+    end
+    return nothing
+end
+
+function subst!(v::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}) where
+        {T<:NumberNotSeries}
+    for k in eachindex(v)
+        subst!(v, a, b, k)
     end
     return nothing
 end

From 1b747aa80d67fe6f64527c175397510bfe473f72 Mon Sep 17 00:00:00 2001
From: Jorge Perez-Hernandez <jorge.perez_hernandez@cfa.harvard.edu>
Date: Thu, 11 Jun 2026 16:01:12 -0400
Subject: [PATCH 17/20] Remove redundant methods; edit existing methods
 appropriately

---
 src/arithmetic.jl | 28 ++++++----------------------
 1 file changed, 6 insertions(+), 22 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index bc4be4a3..c93e3e8e 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -69,8 +69,12 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
                     return nothing
                 end
 
-                function ($fc)(v::$T, a::$T, b::$T, k::Int)
-                    @inbounds v[k] = ($f)(a[k], b[k])
+                @inline function ($fc)(v::$T, a::$T, b::$T, k::Int)
+                    v_coeffs = v.coeffs
+                    a_coeffs = a.coeffs
+                    b_coeffs = b.coeffs
+                    kk = k + 1
+                    @inbounds v_coeffs[kk] = ($f)(a_coeffs[kk], b_coeffs[kk])
                     return nothing
                 end
 
@@ -317,26 +321,6 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
     end
 end
 
-@inline function add!(v::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T},
-        k::Int) where {T<:NumberNotSeries}
-    v_coeffs = v.coeffs
-    a_coeffs = a.coeffs
-    b_coeffs = b.coeffs
-    kk = k + 1
-    @inbounds v_coeffs[kk] = a_coeffs[kk] + b_coeffs[kk]
-    return nothing
-end
-
-@inline function subst!(v::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T},
-        k::Int) where {T<:NumberNotSeries}
-    v_coeffs = v.coeffs
-    a_coeffs = a.coeffs
-    b_coeffs = b.coeffs
-    kk = k + 1
-    @inbounds v_coeffs[kk] = a_coeffs[kk] - b_coeffs[kk]
-    return nothing
-end
-
 function add!(v::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}) where
         {T<:NumberNotSeries}
     for k in eachindex(v)

From b3a000a759d8c41ba52dcbab619a7092a87f929d Mon Sep 17 00:00:00 2001
From: Jorge Perez-Hernandez <jorge.perez_hernandez@cfa.harvard.edu>
Date: Thu, 11 Jun 2026 17:27:31 -0400
Subject: [PATCH 18/20] Move +,add!;-,subst! dispatches into Taylor1
 metaprogramming block

---
 src/arithmetic.jl | 132 ++++++++++++++++------------------------------
 1 file changed, 44 insertions(+), 88 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index c93e3e8e..0cd19e80 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -147,6 +147,50 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
                     return nothing
                 end
 
+                function ($fc)(v::$T{T}, a::$T{T}, b::$T{T}) where
+                        {T<:NumberNotSeries}
+                    for k in eachindex(v)
+                        ($fc)(v, a, b, k)
+                    end
+                    return nothing
+                end
+
+                function ($fc)(v::$T{$T{T}}, a::$T{$T{T}},
+                        b::$T{$T{T}}) where {T<:NumberNotSeries}
+                    v_coeffs = v.coeffs
+                    a_coeffs = a.coeffs
+                    b_coeffs = b.coeffs
+                    @inbounds for i in eachindex(v_coeffs)
+                        ($fc)(v_coeffs[i], a_coeffs[i], b_coeffs[i])
+                    end
+                    return nothing
+                end
+
+                function ($fc)(v::$T{$T{T}}, a::$T{$T{T}},
+                        b::$T{$T{T}}, k::Int) where {T<:NumberNotSeries}
+                    @inbounds ($fc)(v.coeffs[k+1], a.coeffs[k+1], b.coeffs[k+1])
+                    return nothing
+                end
+
+                function ($f)(a::$T{T}, b::$T{T}) where {T<:NumberNotSeries}
+                    if order(a) != order(b)
+                        a, b = fixorder(a, b)
+                    end
+                    c = zero(a)
+                    ($fc)(c, a, b)
+                    return c
+                end
+
+                function ($f)(a::$T{$T{T}}, b::$T{$T{T}}) where
+                        {T<:NumberNotSeries}
+                    if order(a) != order(b)
+                        a, b = fixorder(a, b)
+                    end
+                    c = zero(a)
+                    ($fc)(c, a, b)
+                    return c
+                end
+
             end
         else
             @eval begin
@@ -321,94 +365,6 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
     end
 end
 
-function add!(v::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}) where
-        {T<:NumberNotSeries}
-    for k in eachindex(v)
-        add!(v, a, b, k)
-    end
-    return nothing
-end
-
-function subst!(v::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}) where
-        {T<:NumberNotSeries}
-    for k in eachindex(v)
-        subst!(v, a, b, k)
-    end
-    return nothing
-end
-
-function add!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
-        b::Taylor1{Taylor1{T}}) where {T<:NumberNotSeries}
-    v_coeffs = v.coeffs
-    a_coeffs = a.coeffs
-    b_coeffs = b.coeffs
-    @inbounds for i in eachindex(v_coeffs)
-        add!(v_coeffs[i], a_coeffs[i], b_coeffs[i])
-    end
-    return nothing
-end
-
-function subst!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
-        b::Taylor1{Taylor1{T}}) where {T<:NumberNotSeries}
-    v_coeffs = v.coeffs
-    a_coeffs = a.coeffs
-    b_coeffs = b.coeffs
-    @inbounds for i in eachindex(v_coeffs)
-        subst!(v_coeffs[i], a_coeffs[i], b_coeffs[i])
-    end
-    return nothing
-end
-
-function add!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
-        b::Taylor1{Taylor1{T}}, k::Int) where {T<:NumberNotSeries}
-    @inbounds add!(v.coeffs[k+1], a.coeffs[k+1], b.coeffs[k+1])
-    return nothing
-end
-
-function subst!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
-        b::Taylor1{Taylor1{T}}, k::Int) where {T<:NumberNotSeries}
-    @inbounds subst!(v.coeffs[k+1], a.coeffs[k+1], b.coeffs[k+1])
-    return nothing
-end
-
-function +(a::Taylor1{T}, b::Taylor1{T}) where {T<:NumberNotSeries}
-    if order(a) != order(b)
-        a, b = fixorder(a, b)
-    end
-    c = zero(a)
-    add!(c, a, b)
-    return c
-end
-
-function -(a::Taylor1{T}, b::Taylor1{T}) where {T<:NumberNotSeries}
-    if order(a) != order(b)
-        a, b = fixorder(a, b)
-    end
-    c = zero(a)
-    subst!(c, a, b)
-    return c
-end
-
-function +(a::Taylor1{Taylor1{T}}, b::Taylor1{Taylor1{T}}) where
-        {T<:NumberNotSeries}
-    if order(a) != order(b)
-        a, b = fixorder(a, b)
-    end
-    c = zero(a)
-    add!(c, a, b)
-    return c
-end
-
-function -(a::Taylor1{Taylor1{T}}, b::Taylor1{Taylor1{T}}) where
-        {T<:NumberNotSeries}
-    if order(a) != order(b)
-        a, b = fixorder(a, b)
-    end
-    c = zero(a)
-    subst!(c, a, b)
-    return c
-end
-
 for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
     @eval begin
         function ($fc)(v::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},

From f97c60113f397669c25bb118987a01f2099e6c7a Mon Sep 17 00:00:00 2001
From: Jorge Perez-Hernandez <jorge.perez_hernandez@cfa.harvard.edu>
Date: Thu, 11 Jun 2026 21:16:01 -0400
Subject: [PATCH 19/20] Use constant_term! internally

---
 src/arithmetic.jl      | 28 ++++++++++++++++------------
 src/calculus.jl        |  4 ++--
 src/functions.jl       |  2 +-
 src/other_functions.jl | 36 ++++++++++++++++++------------------
 4 files changed, 37 insertions(+), 33 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 0cd19e80..7dc0885a 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -37,13 +37,13 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
 
             function ($f)(a::$T{T}, b::T) where {T<:Number}
                 c = $T(a.coeffs[:])
-                c[0] = $f(a[0], b)
+                constant_term!(c, $f(constant_term(a), b))
                 return c
             end
 
             function ($f)(b::T, a::$T{T}) where {T<:Number}
                 c = $T($f(a.coeffs[:]))
-                c[0] = $f(b, a[0])
+                constant_term!(c, $f(b, constant_term(a)))
                 return c
             end
 
@@ -199,7 +199,7 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
                     aa = convert(TaylorN{R}, a)
                     bb = convert(R, b)
                     c = TaylorN(aa.coeffs[:], order(aa))
-                    c[0][1] = ($f)(aa[0][1], bb)
+                    constant_term!(c, ($f)(constant_term(aa), bb))
                     return c
                 end
 
@@ -208,7 +208,7 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
                     aa = convert(TaylorN{R}, a)
                     bb = convert(R, b)
                     c = TaylorN(($f)(aa.coeffs[:]), order(aa))
-                    c[0][1] = ($f)(bb, aa[0][1])
+                    constant_term!(c, ($f)(bb, constant_term(aa)))
                     return c
                 end
 
@@ -285,8 +285,9 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
             coeffs = FixedSizeVectorDefault{HomogeneousPolynomial{Taylor1{R}}}(
                     undef, order(a)+1)
             coeffs .= a.coeffs
-            @inbounds coeffs[1] = aux
-            return TaylorN(coeffs, order(a))
+            c = TaylorN(coeffs, order(a))
+            constant_term!(c, aux)
+            return c
         end
 
         function ($f)(b::S, a::TaylorN{Taylor1{T}}) where
@@ -296,8 +297,9 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
             coeffs = FixedSizeVectorDefault{HomogeneousPolynomial{Taylor1{R}}}(
                     undef, order(a)+1)
             coeffs .= $f.(a.coeffs)
-            @inbounds coeffs[1] = aux
-            return TaylorN(coeffs, order(a))
+            c = TaylorN(coeffs, order(a))
+            constant_term!(c, aux)
+            return c
         end
 
         function ($f)(a::TaylorN{Taylor1{T}}, b::Taylor1{S}) where
@@ -306,8 +308,9 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
             R = TS.numtype(aux)
             coeffs = FixedSizeVectorDefault{HomogeneousPolynomial{Taylor1{R}}}(undef, order(a)+1)
             coeffs .= a.coeffs
-            @inbounds coeffs[1] = aux
-            return TaylorN(coeffs, order(a))
+            c = TaylorN(coeffs, order(a))
+            constant_term!(c, aux)
+            return c
         end
 
         function ($f)(b::Taylor1{S}, a::TaylorN{Taylor1{T}}) where
@@ -316,8 +319,9 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
             R = TS.numtype(aux)
             coeffs = FixedSizeVectorDefault{HomogeneousPolynomial{Taylor1{R}}}(undef, order(a)+1)
             coeffs .= $f.(a.coeffs)
-            @inbounds coeffs[1] = aux
-            return TaylorN(coeffs, order(a))
+            c = TaylorN(coeffs, order(a))
+            constant_term!(c, aux)
+            return c
         end
 
         function ($f)(a::Taylor1{TaylorN{T}}, b::S) where
diff --git a/src/calculus.jl b/src/calculus.jl
index a639fa41..77f75c39 100644
--- a/src/calculus.jl
+++ b/src/calculus.jl
@@ -149,7 +149,7 @@ function integrate!(res::Taylor1{T}, a::Taylor1{T}, x::T) where {T<:Number}
     @inbounds for i = 1:order
         res_coeffs[i+1] = a_coeffs[i] / i
     end
-    @inbounds res_coeffs[1] = x
+    constant_term!(res, x)
     return nothing
 end
 integrate!(res::Taylor1{T}, a::Taylor1{T}) where {T<:Number} =
@@ -163,7 +163,7 @@ function integrate!(res::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
     @inbounds for i = 1:order
         div!(res_coeffs[i+1], a_coeffs[i], i)
     end
-    @inbounds res_coeffs[1] = x
+    constant_term!(res, x)
     return nothing
 end
 
diff --git a/src/functions.jl b/src/functions.jl
index 089e433f..27d82aa1 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -148,7 +148,7 @@ for T in (:Taylor1, :TaylorN)
 
         function atan(a::$T, b::$T)
             c = atan(a/b)
-            c[0] = atan(constant_term(a), constant_term(b))
+            constant_term!(c, atan(constant_term(a), constant_term(b)))
             return c
         end
 
diff --git a/src/other_functions.jl b/src/other_functions.jl
index bb81ae32..4b80d7fe 100644
--- a/src/other_functions.jl
+++ b/src/other_functions.jl
@@ -29,9 +29,9 @@ for op in (:mod, :rem)
     for T in (:Taylor1, :TaylorN)
         @eval begin
             function ($op)(a::$T{T}, x::T) where {T<:Real}
-                coeffs = copy(a.coeffs)
-                @inbounds coeffs[1] = ($op)(constant_term(a), x)
-                return $T(coeffs, order(a))
+                c = $T(copy(a.coeffs), order(a))
+                constant_term!(c, ($op)(constant_term(a), x))
+                return c
             end
 
             function ($op)(a::$T{T}, x::S) where {T<:Real,S<:Real}
@@ -44,9 +44,9 @@ for op in (:mod, :rem)
 
     @eval begin
         function ($op)(a::TaylorN{Taylor1{T}}, x::T) where {T<:Real}
-            coeffs = copy(a.coeffs)
-            @inbounds coeffs[1] = ($op)(constant_term(a), x)
-            return TaylorN( coeffs, order(a) )
+            c = TaylorN(copy(a.coeffs), order(a))
+            constant_term!(c, ($op)(constant_term(a), x))
+            return c
         end
 
         function ($op)(a::TaylorN{Taylor1{T}}, x::S) where {T<:Real,S<:Real}
@@ -56,9 +56,9 @@ for op in (:mod, :rem)
         end
 
         function ($op)(a::Taylor1{TaylorN{T}}, x::T) where {T<:Real}
-            coeffs = copy(a.coeffs)
-            @inbounds coeffs[1] = ($op)(constant_term(a), x)
-            return Taylor1( coeffs, order(a) )
+            c = Taylor1(copy(a.coeffs), order(a))
+            constant_term!(c, ($op)(constant_term(a), x))
+            return c
         end
 
         @inbounds function ($op)(a::Taylor1{TaylorN{T}}, x::S) where {T<:Real,S<:Real}
@@ -74,9 +74,9 @@ end
 for T in (:Taylor1, :TaylorN)
     @eval begin
         function mod2pi(a::$T{T}) where {T<:Real}
-            coeffs = copy(a.coeffs)
-            @inbounds coeffs[1] = mod2pi( constant_term(a) )
-            return $T(coeffs, order(a))
+            c = $T(copy(a.coeffs), order(a))
+            constant_term!(c, mod2pi(constant_term(a)))
+            return c
         end
 
         function abs(a::$T{T}) where {T<:Real}
@@ -98,15 +98,15 @@ for T in (:Taylor1, :TaylorN)
 end
 
 function mod2pi(a::TaylorN{Taylor1{T}}) where {T<:Real}
-    coeffs = copy(a.coeffs)
-    @inbounds coeffs[1] = mod2pi( constant_term(a) )
-    return TaylorN( coeffs, order(a) )
+    c = TaylorN(copy(a.coeffs), order(a))
+    constant_term!(c, mod2pi(constant_term(a)))
+    return c
 end
 
 function mod2pi(a::Taylor1{TaylorN{T}}) where {T<:Real}
-    coeffs = copy(a.coeffs)
-    @inbounds coeffs[1] = mod2pi( constant_term(a) )
-    return Taylor1( coeffs, order(a) )
+    c = Taylor1(copy(a.coeffs), order(a))
+    constant_term!(c, mod2pi(constant_term(a)))
+    return c
 end
 
 function abs(a::TaylorN{Taylor1{T}}) where {T<:Real}

From 48fcef199304e1a39327a7612c88e3d6a812a747 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?=
 <perez.hz@gmail.com>
Date: Sun, 7 Jun 2026 23:31:42 -0400
Subject: [PATCH 20/20] Refactor evaluate methods

---
 src/evaluate.jl     | 164 ++++++++++++++++++++++++++++++++++++++++++--
 test/mixtures.jl    |  18 ++++-
 test/onevariable.jl |   2 +
 3 files changed, 176 insertions(+), 8 deletions(-)

diff --git a/src/evaluate.jl b/src/evaluate.jl
index dd2f60fd..87636a09 100644
--- a/src/evaluate.jl
+++ b/src/evaluate.jl
@@ -44,6 +44,33 @@ representing the dependent variables of an ODE, at *time* δt. Note that the
 syntax `x(δt)` is equivalent to `evaluate(x, δt)`, and `x()`
 is equivalent to `evaluate(x)`.
 """
+function evaluate(x::AbstractArray{Taylor1{T}}, δt::S) where
+        {T<:NumberNotSeries, S<:NumberNotSeries}
+    dest = similar(x, promote_type(T, S))
+    evaluate!(x, δt, dest)
+    return dest
+end
+
+function evaluate(x::AbstractArray{Taylor1{Taylor1{T}}}, δt::S) where
+        {T<:NumberNotSeries, S<:NumberNotSeries}
+    dest = similar(x, Taylor1{T})
+    @inbounds for i in eachindex(x, dest)
+        dest[i] = zero(x[i][0])
+    end
+    evaluate!(x, δt, dest)
+    return dest
+end
+
+function evaluate(x::AbstractArray{Taylor1{TaylorN{T}}}, δt::S) where
+        {T<:Number, S<:NumberNotSeries}
+    dest = similar(x, TaylorN{T})
+    @inbounds for i in eachindex(x, dest)
+        dest[i] = zero(x[i][0])
+    end
+    evaluate!(x, δt, dest)
+    return dest
+end
+
 evaluate(x::AbstractArray{Taylor1{T}}, δt::S) where
     {T<:Number, S<:Number} = evaluate.(x, δt)
 
@@ -511,17 +538,140 @@ of an ODE at *time* `δt`. It updates the vector `dest` with the
 computed values.
 """
 function evaluate!(x::AbstractArray{Taylor1{T}}, δt::S,
-        dest::AbstractArray{T}) where {T<:NumberNotSeriesN, S<:Number}
-    dest .= evaluate.( x, δt )
+        dest::AbstractArray{R}) where
+        {T<:NumberNotSeries, S<:NumberNotSeries, R<:NumberNotSeries}
+    @inbounds for i in eachindex(x, dest)
+        dest[i] = evaluate(x[i], δt)
+    end
+    return nothing
+end
+
+function evaluate!(a::Taylor1{Taylor1{T}}, δt::S,
+        dest::Taylor1{T}) where {T<:NumberNotSeries, S<:NumberNotSeries}
+    zero!(dest)
+    @inbounds for k in reverse(eachindex(a))
+        mul!(dest, dest, δt)
+        add!(dest, dest, a[k])
+    end
+    return nothing
+end
+
+function _muladd_taylor1_evaluate!(dest::Taylor1{T}, δt::Taylor1{T},
+        coeff::Taylor1{T}) where {T<:NumberNotSeries}
+    @inbounds for ord in reverse(eachindex(dest))
+        acc = zero(dest[ord])
+        for i in 0:ord
+            acc += dest[i] * δt[ord-i]
+        end
+        dest[ord] = acc + coeff[ord]
+    end
+    return nothing
+end
+
+function evaluate!(a::Taylor1{Taylor1{T}}, δt::Taylor1{T},
+        dest::Taylor1{T}) where {T<:NumberNotSeries}
+    zero!(dest)
+    @inbounds for k in reverse(eachindex(a))
+        _muladd_taylor1_evaluate!(dest, δt, a[k])
+    end
+    return nothing
+end
+
+function evaluate!(x::AbstractArray{Taylor1{Taylor1{T}}}, δt::S,
+        dest::AbstractArray{Taylor1{T}}) where {T<:NumberNotSeries, S<:NumberNotSeries}
+    @inbounds for i in eachindex(x, dest)
+        evaluate!(x[i], δt, dest[i])
+    end
+    return nothing
+end
+
+function evaluate!(x::AbstractArray{Taylor1{Taylor1{T}}}, δt::Taylor1{T},
+        dest::AbstractArray{Taylor1{T}}) where {T<:NumberNotSeries}
+    @inbounds for i in eachindex(x, dest)
+        evaluate!(x[i], δt, dest[i])
+    end
+    return nothing
+end
+
+function evaluate!(a::Taylor1{TaylorN{T}}, δt::S,
+        dest::TaylorN{T}) where {T<:Number, S<:NumberNotSeries}
+    _check_same_space(dest, a[0])
+    zero!(dest)
+    @inbounds for k in reverse(eachindex(a))
+        for ordQ in eachindex(dest)
+            dest_hp = dest[ordQ].coeffs
+            a_hp = a[k][ordQ].coeffs
+            for j in eachindex(dest_hp)
+                dest_hp[j] = dest_hp[j] * δt + a_hp[j]
+            end
+        end
+    end
+    return nothing
+end
+
+@inline function _taylorN_product_coeff(dest::TaylorN{T}, δt::TaylorN{T},
+        ordQ::Int, pos::Int) where {T<:Number}
+    acc = zero(dest[ordQ][pos])
+    @inbounds for i in 0:ordQ
+        dest_hp = dest[i]
+        δt_hp = δt[ordQ-i]
+        if i == 0
+            acc += dest_hp[1] * δt_hp[pos]
+        elseif i == ordQ
+            acc += δt_hp[1] * dest_hp[pos]
+        else
+            table = _init_output_major_product_table!(dest_hp.space, i, ordQ-i)
+            offsets = table.output_offsets
+            output_pairs = table.output_pairs
+            num_right = table.num_right
+            dest_coeffs = dest_hp.coeffs
+            δt_coeffs = δt_hp.coeffs
+            for csr_pos in offsets[pos]:(offsets[pos+1]-1)
+                pair = Int(output_pairs[csr_pos]) - 1
+                na = pair ÷ num_right + 1
+                nb = pair - (na-1) * num_right + 1
+                acc += dest_coeffs[na] * δt_coeffs[nb]
+            end
+        end
+    end
+    return acc
+end
+
+function _muladd_taylorN_evaluate!(dest::TaylorN{T}, δt::TaylorN{T},
+        coeff::TaylorN{T}) where {T<:Number}
+    _check_same_space(dest, δt, coeff)
+    @inbounds for ordQ in reverse(eachindex(dest))
+        dest_hp = dest[ordQ]
+        coeff_hp = coeff[ordQ]
+        for pos in eachindex(dest_hp)
+            dest_hp[pos] = _taylorN_product_coeff(dest, δt, ordQ, pos) + coeff_hp[pos]
+        end
+    end
+    return nothing
+end
+
+function evaluate!(a::Taylor1{TaylorN{T}}, δt::TaylorN{T},
+        dest::TaylorN{T}) where {T<:Number}
+    _check_same_space(dest, δt, a[0])
+    zero!(dest)
+    @inbounds for k in reverse(eachindex(a))
+        _muladd_taylorN_evaluate!(dest, δt, a[k])
+    end
     return nothing
 end
 
 function evaluate!(x::AbstractArray{Taylor1{TaylorN{T}}}, δt::S,
-        dest::AbstractArray{TaylorN{T}}) where {T<:Number, S<:Number}
-    @inbounds for i in eachindex(x)
-        zero!(dest[i])
-        aux = zero(dest[i])
-        _horner!(dest[i], x[i], δt, aux)
+        dest::AbstractArray{TaylorN{T}}) where {T<:Number, S<:NumberNotSeries}
+    @inbounds for i in eachindex(x, dest)
+        evaluate!(x[i], δt, dest[i])
+    end
+    return nothing
+end
+
+function evaluate!(x::AbstractArray{Taylor1{TaylorN{T}}}, δt::TaylorN{T},
+        dest::AbstractArray{TaylorN{T}}) where {T<:Number}
+    @inbounds for i in eachindex(x, dest)
+        evaluate!(x[i], δt, dest[i])
     end
     return nothing
 end
diff --git a/test/mixtures.jl b/test/mixtures.jl
index 89245a1a..d59d4540 100644
--- a/test/mixtures.jl
+++ b/test/mixtures.jl
@@ -160,8 +160,15 @@ using Test
     @test string(evaluate(t1N^2, 1.0)) == " 1.0 + 2.0 x₁ + 1.0 x₁² + 2.0 x₂² + 𝒪(‖x‖³)"
     @test string((t1N^(2//1))(1.0)) == " 1.0 + 2.0 x₁ + 1.0 x₁² + 2.0 x₂² + 𝒪(‖x‖³)"
     v = TaylorN.(zeros(Float64, 2), 2)
-    @test isnothing(evaluate!([t1N, t1N^2], 0.0, v))
+    vt1N = [t1N, t1N^2]
+    @test isnothing(evaluate!(vt1N, 0.0, v))
     @test v == [TaylorN(1), TaylorN(1)^2]
+    evaluate!(vt1N, 0.5, v)
+    @test (@allocated evaluate!(vt1N, 0.5, v)) == 0
+    δtN = 0.5 + t1N[0]
+    evaluate!(vt1N, δtN, v)
+    @test isapprox(v, evaluate(vt1N, δtN))
+    @test (@allocated evaluate!(vt1N, δtN, v)) == 0
     @test tN1() == t
     @test evaluate(tN1, :x₁ => 1.0) == TaylorN([HomogeneousPolynomial([1.0+t]), zero(xHt), yHt^2])
     @test evaluate(tN1, 1, 1.0) == TaylorN([HomogeneousPolynomial([1.0+t]), zero(xHt), yHt^2])
@@ -553,6 +560,15 @@ end
     titii = ti * tii
     @test string(titii) == " ( 1.0 t + 𝒪(t⁴)) s + 𝒪(s¹⁰)"
     @test titii == Taylor1([zero(ti), ti], 9)
+    vtii = [tii, titii]
+    vtii_dest = [zero(ti), zero(ti)]
+    evaluate!(vtii, 0.5, vtii_dest)
+    @test vtii_dest == evaluate(vtii, 0.5)
+    @test (@allocated evaluate!(vtii, 0.5, vtii_dest)) == 0
+    δt1 = 0.5 + ti
+    evaluate!(vtii, δt1, vtii_dest)
+    @test isapprox(vtii_dest, evaluate(vtii, δt1))
+    @test (@allocated evaluate!(vtii, δt1, vtii_dest)) == 0
     @test titii / tii == ti
     @test order(titii/tii) == order(tii)-1
     @test titii / ti == tii
diff --git a/test/onevariable.jl b/test/onevariable.jl
index e6051d0e..d9512c13 100644
--- a/test/onevariable.jl
+++ b/test/onevariable.jl
@@ -471,6 +471,8 @@ Base.iszero(::SymbNumber) = false
     @test evaluate!(v, 0.0, view(vv, 1:2)) == nothing
     @test vv == [0.0,1.0]
     @test evaluate(v) == vv
+    evaluate!(v, 0.2, vv)
+    @test (@allocated evaluate!(v, 0.2, vv)) == 0
     @test isapprox(evaluate(v, complex(0.0,0.2)),
         [complex(0.0,sinh(0.2)),complex(cos(0.2),sin(-0.2))], atol=eps(), rtol=0.0)
     m = [sin(t) exp(-t); cos(t) exp(t)]