From d8989a7ead6580d1f0c67e4c7b9b3cc3a6f4a51b 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: Wed, 31 Jan 2024 23:16:05 +0100
Subject: [PATCH 01/23] WIP: proof of concept: avoid allocations in inplace
 methods

---
 src/arithmetic.jl      | 28 +++++++++++++++++++++++++++-
 test/mutatingfuncts.jl |  4 ++++
 2 files changed, 31 insertions(+), 1 deletion(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 04b0e2f6..a5142909 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -203,6 +203,30 @@ end
 ones(::Type{HomogeneousPolynomial{T}}, order::Int) where {T<:Number} =
     ones( HomogeneousPolynomial([one(T)], 0), order)
 
+function zero!(a::Taylor1{T}) where {T<:NumberNotSeries}
+    a.coeffs .= zero.(a.coeffs)
+    return nothing
+end
+
+function zero!(a::HomogeneousPolynomial{T}) where {T<:NumberNotSeries}
+    a.coeffs .= zero.(a.coeffs)
+    return nothing
+end
+
+function zero!(a::TaylorN{T}) where {T<:NumberNotSeries}
+    for i in 0:a.order
+        zero!(a[i])
+    end
+    return nothing
+end
+
+function zero!(a::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
+    for i in 0:a.order
+        zero!(a[i])
+    end
+    return nothing
+end
+
 
 
 ## Addition and subtraction ##
@@ -479,10 +503,12 @@ for T in (:Taylor1, :TaylorN)
     end
 end
 
+########
+
 function mul!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}},
         ordT::Int) where {T<:NumberNotSeries}
     # Sanity
-    zero!(res, a, ordT)
+    zero!(res[ordT])
     for k in 0:ordT
         @inbounds for ordQ in eachindex(a[ordT])
             mul!(res[ordT], a[k], b[ordT-k], ordQ)
diff --git a/test/mutatingfuncts.jl b/test/mutatingfuncts.jl
index 98de7745..6599391f 100644
--- a/test/mutatingfuncts.jl
+++ b/test/mutatingfuncts.jl
@@ -88,4 +88,8 @@ using Test
     TS.abs2!(res, t1, 2)
     @test res[2] == 1.0
 
+    t2 = Taylor1(Int,15)
+    TaylorSeries.zero!(t2)
+    @test TaylorSeries.iszero(t2)
+
 end

From e1f39b3d2617db6e8f8ddea3d526a50ba4466ff5 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: Thu, 1 Feb 2024 22:54:05 +0100
Subject: [PATCH 02/23] Move zero! methods from arithmetic.jl to functions.jl

---
 src/arithmetic.jl | 26 --------------------------
 src/functions.jl  | 24 ++++++++++++++++++++++++
 2 files changed, 24 insertions(+), 26 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index a5142909..09790541 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -203,30 +203,6 @@ end
 ones(::Type{HomogeneousPolynomial{T}}, order::Int) where {T<:Number} =
     ones( HomogeneousPolynomial([one(T)], 0), order)
 
-function zero!(a::Taylor1{T}) where {T<:NumberNotSeries}
-    a.coeffs .= zero.(a.coeffs)
-    return nothing
-end
-
-function zero!(a::HomogeneousPolynomial{T}) where {T<:NumberNotSeries}
-    a.coeffs .= zero.(a.coeffs)
-    return nothing
-end
-
-function zero!(a::TaylorN{T}) where {T<:NumberNotSeries}
-    for i in 0:a.order
-        zero!(a[i])
-    end
-    return nothing
-end
-
-function zero!(a::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
-    for i in 0:a.order
-        zero!(a[i])
-    end
-    return nothing
-end
-
 
 
 ## Addition and subtraction ##
@@ -503,8 +479,6 @@ for T in (:Taylor1, :TaylorN)
     end
 end
 
-########
-
 function mul!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}},
         ordT::Int) where {T<:NumberNotSeries}
     # Sanity
diff --git a/src/functions.jl b/src/functions.jl
index 4af728bc..de92e7aa 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -655,6 +655,30 @@ end
 
 
 # Mutating functions for Taylor1{TaylorN{T}}
+@inline function zero!(a::Taylor1{T}) where {T<:NumberNotSeries}
+    a.coeffs .= zero.(a.coeffs)
+    return nothing
+end
+
+@inline function zero!(a::HomogeneousPolynomial{T}) where {T<:NumberNotSeries}
+    a.coeffs .= zero.(a.coeffs)
+    return nothing
+end
+
+@inline function zero!(a::TaylorN{T}) where {T<:NumberNotSeries}
+    for i in 0:a.order
+        zero!(a[i])
+    end
+    return nothing
+end
+
+@inline function zero!(a::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
+    for i in 0:a.order
+        zero!(a[i])
+    end
+    return nothing
+end
+
 @inline function exp!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
         k::Int) where {T<:NumberNotSeries}
     if k == 0

From fb339f697f63026ab53f1fef48eccf8a872c0c3d 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: Thu, 1 Feb 2024 23:10:34 +0100
Subject: [PATCH 03/23] Substitute zero! methods where appropriate

---
 src/arithmetic.jl |  2 +-
 src/functions.jl  | 56 +++++++++++++++++++++++------------------------
 src/power.jl      |  6 ++---
 3 files changed, 32 insertions(+), 32 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 09790541..194ab219 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -758,7 +758,7 @@ end
         @inbounds res[0] = cdivfact
         return nothing
     end
-    zero!(res, a, ordT)
+    zero!(res[ordT])
     imin = max(0, ordT+ordfact-b.order)
     aux = TaylorN(zero(a[ordT][0][1]), a[ordT].order)
     for k in imin:ordT-1
diff --git a/src/functions.jl b/src/functions.jl
index de92e7aa..5bc7c765 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -690,7 +690,7 @@ end
     end
     # The recursion formula
     tmp = TaylorN( zero(a[k][0][1]), a[0].order)
-    zero!(res, a, k)
+    zero!(res[k])
     for i = 0:k-1
         @inbounds for ordQ in eachindex(a[0])
             tmp[ordQ] = (k-i) * res[i][ordQ]
@@ -712,7 +712,7 @@ end
     end
     # The recursion formula
     tmp = TaylorN( zero(a[k][0][1]), a[0].order)
-    zero!(res, a, k)
+    zero!(res[k])
     # i=0 term of sum
     @inbounds for ordQ in eachindex(a[0])
         one!(tmp, a[0], ordQ)
@@ -741,14 +741,14 @@ end
         return nothing
     elseif k == 1
         @inbounds for ordQ in eachindex(a[0])
-            zero!(res[k], a[0], ordQ)
+            zero!(res[k][ordQ])
             div!(res[k], a[1], a[0], ordQ)
         end
         return nothing
     end
     # The recursion formula
     tmp = TaylorN( zero(a[k][0][1]), a[0].order)
-    zero!(res, a, k)
+    zero!(res[k])
     for i = 1:k-1
         @inbounds for ordQ in eachindex(a[0])
             tmp[ordQ] = (k-i) * res[k-i][ordQ]
@@ -758,7 +758,7 @@ end
     div!(res, res, k, k)
     @inbounds for ordQ in eachindex(a[0])
         subst!(tmp, a[k], res[k], ordQ)
-        zero!(res[k], a[0], ordQ)
+        zero!(res[k][ordQ])
         div!(res[k], tmp, a[0], ordQ)
     end
     return nothing
@@ -774,7 +774,7 @@ end
         return nothing
     end
     tmp1 = TaylorN( zero(a[k][0][1]), a[0].order)
-    zero!(res, a, k)
+    zero!(res[k])
     @inbounds for ordQ in eachindex(a[0])
         # zero!(res[k], a[0], ordQ)
         one!(tmp1, a[0], ordQ)
@@ -797,7 +797,7 @@ end
     div!(res, res, k, k)
     @inbounds for ordQ in eachindex(a[0])
         subst!(tmp, a[k], res[k], ordQ)
-        zero!(res[k], a[0], ordQ)
+        zero!(res[k][ordQ])
         div!(res[k], tmp, tmp1, ordQ)
     end
     return nothing
@@ -813,8 +813,8 @@ end
     end
     # The recursion formula
     x = TaylorN( a[1][0][1], a[0].order )
-    zero!(s, a, k)
-    zero!(c, a, k)
+    zero!(s[k])
+    zero!(c[k])
     @inbounds for i = 1:k
         for ordQ in eachindex(a[0])
             x[ordQ] = i * a[i][ordQ]
@@ -855,7 +855,7 @@ end
     end
     # The recursion formula
     tmp = TaylorN( zero(a[0][0][1]), a[0].order)
-    zero!(res, a, k)
+    zero!(res[k])
     for i = 0:k-1
         @inbounds for ordQ in eachindex(a[0])
             tmp[ordQ] = (k-i) * res2[i][ordQ]
@@ -889,7 +889,7 @@ end
     end
     # The recursion formula
     tmp = TaylorN( zero(a[0][0][1]), a[0].order)
-    zero!(res, a, k)
+    zero!(res[k])
     for i in 1:k-1
         @inbounds for ordQ in eachindex(a[0])
             tmp[ordQ] = (k-i) * res[k-i][ordQ]
@@ -899,7 +899,7 @@ end
     div!(res, res, k, k)
     @inbounds for ordQ in eachindex(a[0])
         subst!(tmp, a[k], res[k], ordQ)
-        zero!(res[k], a[0], ordQ)
+        zero!(res[k][ordQ])
         div!(res[k], tmp, r[0], ordQ)
     end
     # Compute auxiliary term s=1-a^2
@@ -934,7 +934,7 @@ end
     end
     # The recursion formula
     tmp = TaylorN( zero(a[0][0][1]), a[0].order)
-    zero!(res, a, k)
+    zero!(res[k])
     for i in 1:k-1
         @inbounds for ordQ in eachindex(a[0])
             tmp[ordQ] = (k-i) * res[k-i][ordQ]
@@ -945,7 +945,7 @@ end
     @inbounds for ordQ in eachindex(a[0])
         add!(tmp, a[k], res[k], ordQ)
         subst!(tmp, tmp, ordQ)
-        zero!(res[k], a[0], ordQ)
+        zero!(res[k][ordQ])
         div!(res[k], tmp, r[0], ordQ)
     end
     # Compute auxiliary term s=1-a^2
@@ -974,7 +974,7 @@ end
     end
     # The recursion formula
     tmp = TaylorN( zero(a[0][0][1]), a[0].order )
-    zero!(res, a, k)
+    zero!(res[k])
     for i in 1:k-1
         @inbounds for ordQ in eachindex(a[0])
             tmp[ordQ] = (k-i) * res[k-i][ordQ]
@@ -985,10 +985,10 @@ end
         # zero!(tmp, res[k], ordQ)
         tmp[ordQ] = - res[k][ordQ] / k
         add!(tmp, a[k], tmp, ordQ)
-        zero!(res[k], a[0], ordQ)
+        zero!(res[k][ordQ])
         div!(res[k], tmp, r[0], ordQ)
     end
-    zero!(r, a, k)
+    zero!(r[k])
     sqr!(r, a, k)
     return nothing
 end
@@ -1003,8 +1003,8 @@ end
     end
     # The recursion formula
     x = TaylorN( a[k][0][1], a[0].order )
-    zero!(s, a, k)
-    zero!(c, a, k)
+    zero!(s[k])
+    zero!(c[k])
     @inbounds for i = 1:k
         for ordQ in eachindex(a[0])
             x[ordQ] = i * a[i][ordQ]
@@ -1027,7 +1027,7 @@ end
     end
     # The recursion formula
     tmp = TaylorN( zero(a[0][0][1]), a[0].order)
-    zero!(res, a, k)
+    zero!(res[k])
     for i = 0:k-1
         @inbounds for ordQ in eachindex(a[0])
             tmp[ordQ] = (k-i) * res2[i][ordQ]
@@ -1039,7 +1039,7 @@ end
         tmp[ordQ] = res[k][ordQ] / k
         subst!(res[k], a[k], tmp, ordQ)
     end
-    zero!(res2, res, k)
+    zero!(res2[k])
     sqr!(res2, res, k)
     return nothing
 end
@@ -1061,7 +1061,7 @@ end
     end
     # The recursion formula
     tmp = TaylorN( zero(a[0][0][1]), a[0].order)
-    zero!(res, a, k)
+    zero!(res[k])
     for i in 1:k-1
         @inbounds for ordQ in eachindex(a[0])
             tmp[ordQ] = (k-i) * res[k-i][ordQ]
@@ -1071,7 +1071,7 @@ end
     div!(res, res, k, k)
     @inbounds for ordQ in eachindex(a[0])
         subst!(tmp, a[k], res[k], ordQ)
-        zero!(res[k], a[0], ordQ)
+        zero!(res[k][ordQ])
         div!(res[k], tmp, r[0], ordQ)
     end
     # Compute auxiliary term s=1+a^2
@@ -1105,7 +1105,7 @@ end
     end
     # The recursion formula
     tmp = TaylorN( zero(a[0][0][1]), a[0].order)
-    zero!(res, a, k)
+    zero!(res[k])
     for i in 1:k-1
         @inbounds for ordQ in eachindex(a[0])
             tmp[ordQ] = (k-i) * res[k-i][ordQ]
@@ -1115,7 +1115,7 @@ end
     div!(res, res, k, k)
     @inbounds for ordQ in eachindex(a[0])
         subst!(tmp, a[k], res[k], ordQ)
-        zero!(res[k], a[0], ordQ)
+        zero!(res[k][ordQ])
         div!(res[k], tmp, r[0], ordQ)
     end
     # Compute auxiliary term s=a^2-1
@@ -1143,7 +1143,7 @@ end
     end
     # The recursion formula
     tmp = TaylorN( zero(a[0][0][1]), a[0].order )
-    zero!(res, a, k)
+    zero!(res[k])
     for i in 1:k-1
         @inbounds for ordQ in eachindex(a[0])
             tmp[ordQ] = (k-i) * res[k-i][ordQ]
@@ -1154,10 +1154,10 @@ end
         # zero!(tmp, res[k], ordQ)
         tmp[ordQ] = res[k][ordQ] / k
         add!(tmp, a[k], tmp, ordQ)
-        zero!(res[k], a[0], ordQ)
+        zero!(res[k][ordQ])
         div!(res[k], tmp, r[0], ordQ)
     end
-    zero!(r, a, k)
+    zero!(r[k])
     sqr!(r, a, k)
     return nothing
 end
diff --git a/src/power.jl b/src/power.jl
index 8026eb44..1889e7d9 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -276,7 +276,7 @@ end
     end
 
     # Sanity
-    zero!(res, a, ordT)
+    zero!(res[ordT])
 
     # First non-zero coefficient
     l0 = findfirst(a)
@@ -411,7 +411,7 @@ end
 @inline function sqr!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
         ordT::Int) where {T<:NumberNotSeries}
     # Sanity
-    zero!(res, a, ordT)
+    zero!(res[ordT])
     if ordT == 0
         @inbounds for ordQ in eachindex(a[0])
             @inbounds sqr!(res[0], a[0], ordQ)
@@ -628,7 +628,7 @@ end
 @inline function sqrt!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, ordT::Int,
         ordT0::Int=0) where {T<:NumberNotSeries}
     # Sanity
-    zero!(res, a, ordT)
+    zero!(res[ordT])
     ordT < ordT0 && return nothing
 
     if ordT == ordT0

From bca5bb0303cd6d5c0fef046721c8b4e062d3a356 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: Thu, 1 Feb 2024 23:15:50 +0100
Subject: [PATCH 04/23] Substitute 0:a.order -> eachindex(a)

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

diff --git a/src/functions.jl b/src/functions.jl
index 5bc7c765..0d8f2cf0 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -666,14 +666,14 @@ end
 end
 
 @inline function zero!(a::TaylorN{T}) where {T<:NumberNotSeries}
-    for i in 0:a.order
+    for i in eachindex(a)
         zero!(a[i])
     end
     return nothing
 end
 
 @inline function zero!(a::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
-    for i in 0:a.order
+    for i in eachindex(a)
         zero!(a[i])
     end
     return nothing

From cdbb75eabea6ca0d8e1e2fcc3343b70738b47d2d 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: Fri, 2 Feb 2024 18:44:35 +0100
Subject: [PATCH 05/23] Add one-order zero! methods

---
 src/functions.jl | 36 ++++++++++++++++++++++++++++++------
 1 file changed, 30 insertions(+), 6 deletions(-)

diff --git a/src/functions.jl b/src/functions.jl
index 0d8f2cf0..51f8359e 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -655,26 +655,50 @@ end
 
 
 # Mutating functions for Taylor1{TaylorN{T}}
+@inline function zero!(a::Taylor1{T}, k::Int) where {T<:NumberNotSeries}
+    a[k] = zero(a[k])
+    return nothing
+end
+
 @inline function zero!(a::Taylor1{T}) where {T<:NumberNotSeries}
-    a.coeffs .= zero.(a.coeffs)
+    for k in eachindex(a)
+        zero!(a, k)
+    end
+    return nothing
+end
+
+@inline function zero!(a::HomogeneousPolynomial{T}, k::Int) where {T<:NumberNotSeries}
+    a[k] = zero(a[k])
     return nothing
 end
 
 @inline function zero!(a::HomogeneousPolynomial{T}) where {T<:NumberNotSeries}
-    a.coeffs .= zero.(a.coeffs)
+    for k in eachindex(a)
+        zero!(a, k)
+    end
+    return nothing
+end
+
+@inline function zero!(a::TaylorN{T}, k::Int) where {T<:NumberNotSeries}
+    zero!(a[k])
     return nothing
 end
 
 @inline function zero!(a::TaylorN{T}) where {T<:NumberNotSeries}
-    for i in eachindex(a)
-        zero!(a[i])
+    for k in eachindex(a)
+        zero!(a, k)
     end
     return nothing
 end
 
+@inline function zero!(a::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries}
+    zero!(a[k])
+    return nothing
+end
+
 @inline function zero!(a::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
-    for i in eachindex(a)
-        zero!(a[i])
+    for k in eachindex(a)
+        zero!(a, k)
     end
     return nothing
 end

From 7b3bf251b59c401348fb0d6fcd864295fdef888b 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: Fri, 2 Feb 2024 22:32:25 +0100
Subject: [PATCH 06/23] Try allocation-free zero! for non-mixtures

---
 src/functions.jl | 50 +++++++++++++++++++++++++++++++-----------------
 src/power.jl     | 10 +++++-----
 2 files changed, 37 insertions(+), 23 deletions(-)

diff --git a/src/functions.jl b/src/functions.jl
index 51f8359e..eb2306d2 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -247,7 +247,7 @@ for T in (:Taylor1, :TaylorN)
         end
 
         @inline function one!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
-            zero!(c, a, k)
+            zero!(c, k)
             if k == 0
                 @inbounds c[0] = one(a[0])
             end
@@ -275,7 +275,7 @@ for T in (:Taylor1, :TaylorN)
                 @inbounds c[0] = exp(constant_term(a))
                 return nothing
             end
-            zero!(c, a, k)
+            zero!(c, k)
             if $T == Taylor1
                 @inbounds c[k] = k * a[k] * c[0]
             else
@@ -297,7 +297,7 @@ for T in (:Taylor1, :TaylorN)
                 @inbounds c[0] = expm1(constant_term(a))
                 return nothing
             end
-            zero!(c, a, k)
+            zero!(c, k)
             c0 = c[0]+one(c[0])
             if $T == Taylor1
                 @inbounds c[k] = k * a[k] * c0
@@ -323,7 +323,7 @@ for T in (:Taylor1, :TaylorN)
                 @inbounds c[1] = a[1] / constant_term(a)
                 return nothing
             end
-            zero!(c, a, k)
+            zero!(c, k)
             if $T == Taylor1
                 @inbounds c[k] = (k-1) * a[1] * c[k-1]
             else
@@ -353,7 +353,7 @@ for T in (:Taylor1, :TaylorN)
             end
             a0 = constant_term(a)
             a0p1 = a0+one(a0)
-            zero!(c, a, k)
+            zero!(c, k)
             if $T == Taylor1
                 @inbounds c[k] = (k-1) * a[1] * c[k-1]
             else
@@ -377,8 +377,8 @@ for T in (:Taylor1, :TaylorN)
                 return nothing
             end
             x = a[1]
-            zero!(s, a, k)
-            zero!(c, a, k)
+            zero!(s, k)
+            zero!(c, k)
             if $T == Taylor1
                 @inbounds s[k] = x * c[k-1]
                 @inbounds c[k] = -x * s[k-1]
@@ -423,7 +423,7 @@ for T in (:Taylor1, :TaylorN)
                 @inbounds c2[0] = aux^2
                 return nothing
             end
-            zero!(c, a, k)
+            zero!(c, k)
             if $T == Taylor1
                 @inbounds c[k] = k * a[k] * c2[0]
             else
@@ -449,7 +449,7 @@ for T in (:Taylor1, :TaylorN)
                 @inbounds r[0] = sqrt( 1 - a0^2 )
                 return nothing
             end
-            zero!(c, a, k)
+            zero!(c, k)
             if $T == Taylor1
                 @inbounds c[k] = (k-1) * r[1] * c[k-1]
             else
@@ -474,7 +474,7 @@ for T in (:Taylor1, :TaylorN)
                 @inbounds r[0] = sqrt( 1 - a0^2 )
                 return nothing
             end
-            zero!(c, a, k)
+            zero!(c, k)
             if $T == Taylor1
                 @inbounds c[k] = (k-1) * r[1] * c[k-1]
             else
@@ -499,7 +499,7 @@ for T in (:Taylor1, :TaylorN)
                 @inbounds r[0] = 1 + a0^2
                 return nothing
             end
-            zero!(c, a, k)
+            zero!(c, k)
             if $T == Taylor1
                 @inbounds c[k] = (k-1) * r[1] * c[k-1]
             else
@@ -524,8 +524,8 @@ for T in (:Taylor1, :TaylorN)
                 return nothing
             end
             x = a[1]
-            zero!(s, a, k)
-            zero!(c, a, k)
+            zero!(s, k)
+            zero!(c, k)
             if $T == Taylor1
                 @inbounds s[k] = x * c[k-1]
                 @inbounds c[k] = x * s[k-1]
@@ -555,7 +555,7 @@ for T in (:Taylor1, :TaylorN)
                 @inbounds c2[0] = aux^2
                 return nothing
             end
-            zero!(c, a, k)
+            zero!(c, k)
             if $T == Taylor1
                 @inbounds c[k] = k * a[k] * c2[0]
             else
@@ -581,7 +581,7 @@ for T in (:Taylor1, :TaylorN)
                 @inbounds r[0] = sqrt( a0^2 + 1 )
                 return nothing
             end
-            zero!(c, a, k)
+            zero!(c, k)
             if $T == Taylor1
                 @inbounds c[k] = (k-1) * r[1] * c[k-1]
             else
@@ -606,7 +606,7 @@ for T in (:Taylor1, :TaylorN)
                 @inbounds r[0] = sqrt( a0^2 - 1 )
                 return nothing
             end
-            zero!(c, a, k)
+            zero!(c, k)
             if $T == Taylor1
                 @inbounds c[k] = (k-1) * r[1] * c[k-1]
             else
@@ -631,7 +631,7 @@ for T in (:Taylor1, :TaylorN)
                 @inbounds r[0] = 1 - a0^2
                 return nothing
             end
-            zero!(c, a, k)
+            zero!(c, k)
             if $T == Taylor1
                 @inbounds c[k] = (k-1) * r[1] * c[k-1]
             else
@@ -660,6 +660,13 @@ end
     return nothing
 end
 
+@inline function zero!(a::Taylor1{Taylor1{T}}, k::Int) where {T<:NumberNotSeries}
+    for l in eachindex(a[k])
+        zero!(a[k], l)
+    end
+    return nothing
+end
+
 @inline function zero!(a::Taylor1{T}) where {T<:NumberNotSeries}
     for k in eachindex(a)
         zero!(a, k)
@@ -703,6 +710,13 @@ end
     return nothing
 end
 
+@inline function zero!(a::TaylorN{Taylor1{T}}, k::Int) where {T<:NumberNotSeries}
+    for l in eachindex(a[k])
+        zero!(a[k][l])
+    end
+    return nothing
+end
+
 @inline function exp!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
         k::Int) where {T<:NumberNotSeries}
     if k == 0
@@ -841,7 +855,7 @@ end
     zero!(c[k])
     @inbounds for i = 1:k
         for ordQ in eachindex(a[0])
-            x[ordQ] = i * a[i][ordQ]
+            x[ordQ].coeffs .= i .* a[i][ordQ].coeffs
             mul!(s[k], x, c[k-i], ordQ)
             mul!(c[k], x, s[k-i], ordQ)
         end
diff --git a/src/power.jl b/src/power.jl
index 1889e7d9..838255bf 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -196,7 +196,7 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`.
     end
 
     # Sanity
-    zero!(c, a, k)
+    zero!(c, k)
 
     # First non-zero coefficient
     l0 = findfirst(a)
@@ -250,7 +250,7 @@ end
     end
 
     # Sanity
-    zero!(c, a, k)
+    zero!(c, k)
 
     # The recursion formula
     for i = 0:k-1
@@ -382,7 +382,7 @@ for T = (:Taylor1, :TaylorN)
             end
 
             # Sanity
-            zero!(c, a, k)
+            zero!(c, k)
 
             # Recursion formula
             kodd = k%2
@@ -569,7 +569,7 @@ coefficient, which must be even.
 
 @inline function sqrt!(c::Taylor1{T}, a::Taylor1{T}, k::Int, k0::Int=0) where {T<:Number}
     # Sanity
-    zero!(c, a, k)
+    zero!(c, k)
 
     k < k0 && return nothing
 
@@ -603,7 +603,7 @@ end
 
 @inline function sqrt!(c::TaylorN{T}, a::TaylorN{T}, k::Int) where {T<:NumberNotSeriesN}
     # Sanity
-    zero!(c, a, k)
+    zero!(c, k)
 
     if k == 0
         @inbounds c[0] = sqrt( constant_term(a) )

From f60296e43f0904a2ee4b755608d7608e3a577c6c 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: Thu, 8 Feb 2024 20:01:34 +0100
Subject: [PATCH 07/23] WIP: debugging zero stuff

---
 Project.toml        | 13 +++++++------
 src/TaylorSeries.jl |  2 +-
 src/arithmetic.jl   |  8 ++++++--
 src/power.jl        |  9 ++++++---
 4 files changed, 20 insertions(+), 12 deletions(-)

diff --git a/Project.toml b/Project.toml
index b781e891..550508ef 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,11 +4,18 @@ repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
 version = "0.16.0"
 
 [deps]
+InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
+[weakdeps]
+IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+
+[extensions]
+TaylorSeriesIAExt = "IntervalArithmetic"
+
 [compat]
 Aqua = "0.8"
 IntervalArithmetic = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20"
@@ -19,9 +26,6 @@ SparseArrays = "<0.0.1, 1"
 Test = "<0.0.1, 1"
 julia = "1"
 
-[extensions]
-TaylorSeriesIAExt = "IntervalArithmetic"
-
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
@@ -31,6 +35,3 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["IntervalArithmetic", "LinearAlgebra", "SparseArrays", "Test", "Aqua"]
-
-[weakdeps]
-IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl
index dca66013..4c3988a8 100644
--- a/src/TaylorSeries.jl
+++ b/src/TaylorSeries.jl
@@ -16,7 +16,7 @@ see also [`HomogeneousPolynomial`](@ref).
 """
 module TaylorSeries
 
-
+using InteractiveUtils
 using SparseArrays: SparseMatrixCSC
 using Markdown
 
diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 194ab219..99b70d22 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -482,7 +482,7 @@ end
 function mul!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}},
         ordT::Int) where {T<:NumberNotSeries}
     # Sanity
-    zero!(res[ordT])
+    zero!(res, ordT)
     for k in 0:ordT
         @inbounds for ordQ in eachindex(a[ordT])
             mul!(res[ordT], a[k], b[ordT-k], ordQ)
@@ -736,15 +736,19 @@ div!(v::Taylor1, b::NumberNotSeries, a::Taylor1, k::Int) =
 
 # NOTE: Here `div!` *accumulates* the result of a / b in c[k] (k > 0)
 @inline function div!(c::TaylorN, a::TaylorN, b::TaylorN, k::Int)
+    @show c
     if k==0
         @inbounds c[0] = a[0] / constant_term(b)
         return nothing
     end
+    @show c
 
     @inbounds for i = 0:k-1
         mul!(c[k], c[i], b[k-i])
     end
+    @show c
     @inbounds c[k] = (a[k] - c[k]) / constant_term(b)
+    @show c
     return nothing
 end
 
@@ -758,7 +762,7 @@ end
         @inbounds res[0] = cdivfact
         return nothing
     end
-    zero!(res[ordT])
+    zero!(res, ordT)
     imin = max(0, ordT+ordfact-b.order)
     aux = TaylorN(zero(a[ordT][0][1]), a[ordT].order)
     for k in imin:ordT-1
diff --git a/src/power.jl b/src/power.jl
index 838255bf..abcb6b85 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -276,7 +276,7 @@ end
     end
 
     # Sanity
-    zero!(res[ordT])
+    zero!(res, ordT)
 
     # First non-zero coefficient
     l0 = findfirst(a)
@@ -311,6 +311,7 @@ end
     end
     @inbounds for ordQ in eachindex(tmp)
         tmp1[ordQ] = tmp[ordQ]/kprime
+        @show @which div!(res[ordT], tmp1, a[l0], ordQ)
         div!(res[ordT], tmp1, a[l0], ordQ)
     end
 
@@ -411,7 +412,7 @@ end
 @inline function sqr!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
         ordT::Int) where {T<:NumberNotSeries}
     # Sanity
-    zero!(res[ordT])
+    zero!(res, ordT)
     if ordT == 0
         @inbounds for ordQ in eachindex(a[0])
             @inbounds sqr!(res[0], a[0], ordQ)
@@ -628,7 +629,7 @@ end
 @inline function sqrt!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, ordT::Int,
         ordT0::Int=0) where {T<:NumberNotSeries}
     # Sanity
-    zero!(res[ordT])
+    zero!(res, ordT)
     ordT < ordT0 && return nothing
 
     if ordT == ordT0
@@ -656,6 +657,7 @@ end
         end
     end
 
+    @show res
     @inbounds for ordQ in eachindex(a[0])
         tmp[ordQ] = -2 * tmp[ordQ]
         if ordT+ordT0 ≤ a.order
@@ -668,6 +670,7 @@ end
         tmp1[ordQ] = 2*res[ordT0][ordQ]
         div!(res[ordT], tmp, tmp1, ordQ)
     end
+    @show res
 
     return nothing
 end

From 75628c7796755fadbfcdb49573133acccad3df03 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: Fri, 9 Feb 2024 20:03:37 +0100
Subject: [PATCH 08/23] WIP: workaround TaylorIntegration bug

---
 Project.toml        |  1 -
 src/TaylorSeries.jl |  1 -
 src/arithmetic.jl   | 25 +++++++++++++++++++++----
 src/power.jl        |  3 ---
 4 files changed, 21 insertions(+), 9 deletions(-)

diff --git a/Project.toml b/Project.toml
index 550508ef..134ce35a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,7 +4,6 @@ repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
 version = "0.16.0"
 
 [deps]
-InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl
index 4c3988a8..440fad9b 100644
--- a/src/TaylorSeries.jl
+++ b/src/TaylorSeries.jl
@@ -16,7 +16,6 @@ see also [`HomogeneousPolynomial`](@ref).
 """
 module TaylorSeries
 
-using InteractiveUtils
 using SparseArrays: SparseMatrixCSC
 using Markdown
 
diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 99b70d22..8a1d43b5 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -736,19 +736,15 @@ div!(v::Taylor1, b::NumberNotSeries, a::Taylor1, k::Int) =
 
 # NOTE: Here `div!` *accumulates* the result of a / b in c[k] (k > 0)
 @inline function div!(c::TaylorN, a::TaylorN, b::TaylorN, k::Int)
-    @show c
     if k==0
         @inbounds c[0] = a[0] / constant_term(b)
         return nothing
     end
-    @show c
 
     @inbounds for i = 0:k-1
         mul!(c[k], c[i], b[k-i])
     end
-    @show c
     @inbounds c[k] = (a[k] - c[k]) / constant_term(b)
-    @show c
     return nothing
 end
 
@@ -787,7 +783,28 @@ end
     return nothing
 end
 
+@inline 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]
+        end
+    end
+    return nothing
+end
+
+mul!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
+    b::NumberNotSeries, k::Int) where {T<:NumberNotSeries} = mul!(res, b, a, k)
 
+@inline function div!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
+        b::NumberNotSeries, k::Int) where {T<:NumberNotSeries}
+    for l in eachindex(b[k])
+        for m in eachindex(b[k][l])
+            res[k][l][m] = a[k][l][m]/b
+        end
+    end
+    return nothing
+end
 
 
 """
diff --git a/src/power.jl b/src/power.jl
index abcb6b85..cb3da8db 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -311,7 +311,6 @@ end
     end
     @inbounds for ordQ in eachindex(tmp)
         tmp1[ordQ] = tmp[ordQ]/kprime
-        @show @which div!(res[ordT], tmp1, a[l0], ordQ)
         div!(res[ordT], tmp1, a[l0], ordQ)
     end
 
@@ -657,7 +656,6 @@ end
         end
     end
 
-    @show res
     @inbounds for ordQ in eachindex(a[0])
         tmp[ordQ] = -2 * tmp[ordQ]
         if ordT+ordT0 ≤ a.order
@@ -670,7 +668,6 @@ end
         tmp1[ordQ] = 2*res[ordT0][ordQ]
         div!(res[ordT], tmp, tmp1, ordQ)
     end
-    @show res
 
     return nothing
 end

From 20771fc1ada3b8e21774b41f11a00a8cb6604752 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: Fri, 9 Feb 2024 20:48:14 +0100
Subject: [PATCH 09/23] Add tests

---
 test/mixtures.jl | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/test/mixtures.jl b/test/mixtures.jl
index 02792a25..242d946d 100644
--- a/test/mixtures.jl
+++ b/test/mixtures.jl
@@ -480,6 +480,16 @@ using Test
     c = Taylor1(constant_term(b),0)
     c[0][0][1] = 0.0
     b == bcopy
+
+    #347
+    a = Taylor1([1.0+X,-X, Y, X-Y,X])
+    b = deepcopy(a)
+    b[0] = zero(a[0])
+    b.coeffs[2:end] .= zero(b.coeffs[1])
+    @test iszero(b)
+    @test b.coeffs[2] === b.coeffs[3]
+    b.coeffs[2:end] .= zero.(b.coeffs[1])
+    @test !(b.coeffs[2] === b.coeffs[3])
 end
 
 @testset "Tests with nested Taylor1s" begin

From 0187b2d5ef6f0e83f472296efc98dbf8351f7967 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: Fri, 9 Feb 2024 20:54:00 +0100
Subject: [PATCH 10/23] Fix typo

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

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 8a1d43b5..15e09a4c 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -798,8 +798,8 @@ mul!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
 
 @inline function div!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
         b::NumberNotSeries, k::Int) where {T<:NumberNotSeries}
-    for l in eachindex(b[k])
-        for m in eachindex(b[k][l])
+    for l in eachindex(a[k])
+        for m in eachindex(a[k][l])
             res[k][l][m] = a[k][l][m]/b
         end
     end

From 90c50967a089712bc856c376de24073490efa848 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, 10 Feb 2024 21:55:07 +0100
Subject: [PATCH 11/23] Small fix in div!

---
 src/arithmetic.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 15e09a4c..ed2fafab 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -732,7 +732,7 @@ end
 end
 
 div!(v::Taylor1, b::NumberNotSeries, a::Taylor1, k::Int) =
-    div!(v::Taylor1, Taylor1(b, a.order), a, k)
+    div!(v::Taylor1, b*one(a), a, k)
 
 # NOTE: Here `div!` *accumulates* the result of a / b in c[k] (k > 0)
 @inline function div!(c::TaylorN, a::TaylorN, b::TaylorN, k::Int)

From 31ce07fe3c738f7f4c865c303cc441cf4a59d777 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, 11 Feb 2024 12:29:29 +0100
Subject: [PATCH 12/23] Allocate less in
 div!(::Taylor1,::NumberNotSeries,::Taylor1,Int)

---
 src/arithmetic.jl | 41 +++++++++++++++++++++++++++++++++++++++--
 1 file changed, 39 insertions(+), 2 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index ed2fafab..688124f0 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -731,8 +731,24 @@ end
     return nothing
 end
 
-div!(v::Taylor1, b::NumberNotSeries, a::Taylor1, k::Int) =
-    div!(v::Taylor1, b*one(a), a, k)
+@inline function div!(c::Taylor1, a::NumberNotSeries,
+        b::Taylor1, k::Int)
+    # 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]
+        return nothing
+    end
+
+    imin = max(0, k-b.order)
+    @inbounds c[k] = c[imin] * b[k-imin]
+    @inbounds for i = imin+1:k-1
+        c[k] += c[i] * b[k-i]
+    end
+    @inbounds c[k] = -c[k]/b[0]
+    return nothing
+end
 
 # NOTE: Here `div!` *accumulates* the result of a / b in c[k] (k > 0)
 @inline function div!(c::TaylorN, a::TaylorN, b::TaylorN, k::Int)
@@ -806,6 +822,27 @@ mul!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
     return nothing
 end
 
+# ### TODO: use non-allocating div!(::TaylorN,::TaylorN,::TaylorN)
+# @inline function div!(c::Taylor1{TaylorN{T}}, a::NumberNotSeries,
+#         b::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries}
+#     # 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]
+#         return nothing
+#     end
+
+#     imin = max(0, k-b.order)
+#     @inbounds c[k] = c[imin] * b[k-imin]
+#     @inbounds for i = imin+1:k-1
+#         c[k] += c[i] * b[k-i]
+#     end
+#     @inbounds c[k] = -c[k]/b[0]
+#     return nothing
+# end
+
+
 
 """
     mul!(Y, A, B)

From 8443ff4cf28963552a2e254e60c80e53ef1890ea 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, 11 Feb 2024 14:27:40 +0100
Subject: [PATCH 13/23] Add missing div! method

---
 src/functions.jl | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/src/functions.jl b/src/functions.jl
index eb2306d2..50672f57 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -667,6 +667,13 @@ end
     return nothing
 end
 
+@inline function zero!(a::Taylor1{Taylor1{T}}) where {T<:NumberNotSeries}
+    for k in eachindex(a)
+        zero!(a, k)
+    end
+    return nothing
+end
+
 @inline function zero!(a::Taylor1{T}) where {T<:NumberNotSeries}
     for k in eachindex(a)
         zero!(a, k)

From 83c581c108dea0107164e7cb0cb51f0bb7e91a93 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, 11 Feb 2024 15:23:52 +0100
Subject: [PATCH 14/23] Add
 /(::NumberNotSeries,::Taylor1{TaylorN{<:NumberNotSeries}}

---
 src/arithmetic.jl | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 688124f0..9264ca55 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -665,6 +665,19 @@ function /(a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}}) where {T<:NumberNotSe
     return res
 end
 
+function /(a::NumberNotSeries, b::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
+    iszero(a) && !iszero(b) && return zero(a)
+
+    # 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]
+
+    res = Taylor1(a/b[0], b.order)
+    for ordT in eachindex(res)
+        div!(res, a, b, ordT)
+    end
+    return res
+end
 
 @inline function divfactorization(a1::Taylor1, b1::Taylor1)
     # order of first factorized term; a1 and b1 assumed to be of the same order
@@ -733,6 +746,7 @@ end
 
 @inline function div!(c::Taylor1, a::NumberNotSeries,
         b::Taylor1, k::Int)
+    iszero(a) && !iszero(b) && zero!(c, k)
     # 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]

From e6d9167b3410ec8b376d4bd001fe3f260166888e 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: Mon, 12 Feb 2024 23:19:52 +0100
Subject: [PATCH 15/23] Update comments

---
 src/arithmetic.jl | 21 +--------------------
 1 file changed, 1 insertion(+), 20 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 9264ca55..cd6ca7fe 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -746,6 +746,7 @@ end
 
 @inline function div!(c::Taylor1, a::NumberNotSeries,
         b::Taylor1, k::Int)
+    ### TODO: use non-allocating div!(::TaylorN,::TaylorN,::TaylorN)
     iszero(a) && !iszero(b) && zero!(c, k)
     # order and coefficient of first factorized term
     # In this case, since a[k]=0 for k>0, we can simplify to:
@@ -836,26 +837,6 @@ mul!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
     return nothing
 end
 
-# ### TODO: use non-allocating div!(::TaylorN,::TaylorN,::TaylorN)
-# @inline function div!(c::Taylor1{TaylorN{T}}, a::NumberNotSeries,
-#         b::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries}
-#     # 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]
-#         return nothing
-#     end
-
-#     imin = max(0, k-b.order)
-#     @inbounds c[k] = c[imin] * b[k-imin]
-#     @inbounds for i = imin+1:k-1
-#         c[k] += c[i] * b[k-i]
-#     end
-#     @inbounds c[k] = -c[k]/b[0]
-#     return nothing
-# end
-
 
 
 """

From c78eb49e19ecd778441f24b4274c778c37c04263 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: Wed, 14 Feb 2024 07:02:36 +0100
Subject: [PATCH 16/23] Address review by @lbenet

---
 src/arithmetic.jl | 53 ++++++++++++++++++++++++++++++++---------------
 src/functions.jl  |  3 ++-
 2 files changed, 38 insertions(+), 18 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index cd6ca7fe..fed9e1fb 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -665,8 +665,9 @@ function /(a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}}) where {T<:NumberNotSe
     return res
 end
 
-function /(a::NumberNotSeries, b::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
-    iszero(a) && !iszero(b) && return zero(a)
+function /(a::S, b::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries, S<:Union{T,TaylorN{T}}}
+    R = promote_type(S, TaylorN{T})
+    iszero(a) && !iszero(b) && return convert(Taylor1{R}, zero(b))
 
     # order and coefficient of first factorized term
     # In this case, since a[k]=0 for k>0, we can simplify to:
@@ -744,8 +745,27 @@ end
     return nothing
 end
 
-@inline function div!(c::Taylor1, a::NumberNotSeries,
-        b::Taylor1, k::Int)
+@inline function div!(c::Taylor1{T}, a::NumberNotSeries,
+        b::Taylor1{T}, k::Int) where {T<:Number}
+    iszero(a) && !iszero(b) && zero!(c, k)
+    # 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]
+        return nothing
+    end
+
+    @inbounds c[k] = c[0] * b[k]
+    @inbounds for i = 1:k-1
+        c[k] += c[i] * b[k-i]
+    end
+    @inbounds c[k] = -c[k]/b[0]
+    return nothing
+end
+
+@inline function div!(c::Taylor1{TaylorN{T}}, a::NumberNotSeries,
+        b::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries}
     ### TODO: use non-allocating div!(::TaylorN,::TaylorN,::TaylorN)
     iszero(a) && !iszero(b) && zero!(c, k)
     # order and coefficient of first factorized term
@@ -756,9 +776,8 @@ end
         return nothing
     end
 
-    imin = max(0, k-b.order)
-    @inbounds c[k] = c[imin] * b[k-imin]
-    @inbounds for i = imin+1:k-1
+    @inbounds c[k] = c[0] * b[k]
+    @inbounds for i = 1:k-1
         c[k] += c[i] * b[k-i]
     end
     @inbounds c[k] = -c[k]/b[0]
@@ -814,6 +833,16 @@ end
     return nothing
 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
+        end
+    end
+    return nothing
+end
+
 @inline function mul!(res::Taylor1{TaylorN{T}}, a::NumberNotSeries,
         b::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries}
     for l in eachindex(b[k])
@@ -827,16 +856,6 @@ end
 mul!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
     b::NumberNotSeries, k::Int) where {T<:NumberNotSeries} = mul!(res, b, a, k)
 
-@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
-        end
-    end
-    return nothing
-end
-
 
 
 """
diff --git a/src/functions.jl b/src/functions.jl
index 50672f57..db1e6540 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -242,7 +242,8 @@ for T in (:Taylor1, :TaylorN)
         end
 
         @inline function zero!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
-            @inbounds c[k] = zero(a[k])
+            # @inbounds c[k] = zero(a[k])
+            @inbounds zero!(c, k)
             return nothing
         end
 

From b9f9b134ad09849b6e3675e36cadbdc3ce01334f 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: Thu, 15 Feb 2024 12:47:20 +0100
Subject: [PATCH 17/23] WIP: try less allocations in div!

---
 Project.toml        |  2 ++
 src/TaylorSeries.jl |  1 +
 src/arithmetic.jl   | 72 ++++++++++++++++++++++++++++++++++++++++++---
 3 files changed, 71 insertions(+), 4 deletions(-)

diff --git a/Project.toml b/Project.toml
index 134ce35a..2bbdbfc0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,6 +4,7 @@ repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
 version = "0.16.0"
 
 [deps]
+InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
@@ -17,6 +18,7 @@ TaylorSeriesIAExt = "IntervalArithmetic"
 
 [compat]
 Aqua = "0.8"
+InteractiveUtils = "1"
 IntervalArithmetic = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20"
 LinearAlgebra = "<0.0.1, 1"
 Markdown = "<0.0.1, 1"
diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl
index 440fad9b..4c3988a8 100644
--- a/src/TaylorSeries.jl
+++ b/src/TaylorSeries.jl
@@ -16,6 +16,7 @@ see also [`HomogeneousPolynomial`](@ref).
 """
 module TaylorSeries
 
+using InteractiveUtils
 using SparseArrays: SparseMatrixCSC
 using Markdown
 
diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index fed9e1fb..c6b2eda3 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -772,32 +772,96 @@ end
     # 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 div!(c[0], a, b[0])
         return nothing
     end
 
-    @inbounds c[k] = c[0] * b[k]
+    @inbounds mul!(c[k], c[0], b[k])
     @inbounds for i = 1:k-1
         c[k] += c[i] * b[k-i]
     end
     @inbounds c[k] = -c[k]/b[0]
+    # @inbounds divsubst!(c[k], b[0])
     return nothing
 end
 
 # NOTE: Here `div!` *accumulates* the result of a / b in c[k] (k > 0)
 @inline function div!(c::TaylorN, a::TaylorN, b::TaylorN, k::Int)
     if k==0
-        @inbounds c[0] = a[0] / constant_term(b)
+        @inbounds c[0][1] = constant_term(a) / constant_term(b)
         return nothing
     end
 
     @inbounds for i = 0:k-1
         mul!(c[k], c[i], b[k-i])
     end
-    @inbounds c[k] = (a[k] - c[k]) / constant_term(b)
+    @inbounds for i in eachindex(c[k])
+        c[k][i] = (a[k][i] - c[k][i]) / constant_term(b)
+    end
     return nothing
 end
 
+# NOTE: Here `div!` *accumulates* the result of a / b in c[k] (k > 0)
+@inline function div!(c::TaylorN, a::NumberNotSeries, b::TaylorN, k::Int)
+    if k==0
+        @inbounds c[0][1] = a / constant_term(b)
+        return nothing
+    end
+
+    @inbounds for i = 0:k-1
+        mul!(c[k], c[i], b[k-i])
+    end
+    @inbounds for i in eachindex(c[k])
+        c[k][i] = ( -c[k][i] ) / constant_term(b)
+    end
+    return nothing
+end
+
+# in-place division (assumes equal order among TaylorNs)
+function div!(c::TaylorN, a::TaylorN, b::TaylorN)
+    for k in eachindex(c)
+        div!(c, a, b, k)
+    end
+end
+
+function div!(c::TaylorN, a::NumberNotSeries, b::TaylorN)
+    for k in eachindex(c)
+        div!(c, a, b, k)
+    end
+end
+
+# inplace multiplication TaylorNs
+function mul!(c::TaylorN, a::TaylorN, b::TaylorN)
+    for k in eachindex(c)
+        mul!(c, a, b, k)
+    end
+end
+
+# # inplace method for computing c = -c/b
+# # NOTE: Here `divsubst!` *accumulates* the result inplace
+# @inline function divsubst!(c::TaylorN, b::TaylorN, k::Int)
+#     if k==0
+#         @inbounds for i in eachindex(c[0])
+#             c[0][i] = -c[0][i] / constant_term(b)
+#         end
+#         return nothing
+#     end
+
+#     @inbounds for i = 0:k-1
+#         mul!(c[k], c[i], b[k-i])
+#     end
+#     @inbounds for i in eachindex(c[k])
+#         c[k][i] = (-2c[k][i]) / constant_term(b)
+#     end
+#     return nothing
+# end
+
+# function divsubst!(c::TaylorN, b::TaylorN)
+#     for k in eachindex(c)
+#         divsubst!(c, b, k)
+#     end
+# end
+
 # NOTE: Here `div!` *accumulates* the result of a / b in res[k] (k > 0)
 @inline function div!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
         b::Taylor1{TaylorN{T}}, ordT::Int) where {T<:NumberNotSeriesN}

From 6f05d471bd6023aebbb88403bd950ef843c74bfb 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: Fri, 16 Feb 2024 14:27:31 +0100
Subject: [PATCH 18/23] Fixes in division methods

---
 Project.toml        |  2 --
 src/TaylorSeries.jl |  1 -
 src/arithmetic.jl   | 23 +++++++++++++++++++----
 test/mixtures.jl    |  4 ++++
 4 files changed, 23 insertions(+), 7 deletions(-)

diff --git a/Project.toml b/Project.toml
index 2bbdbfc0..134ce35a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,7 +4,6 @@ repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
 version = "0.16.0"
 
 [deps]
-InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
@@ -18,7 +17,6 @@ TaylorSeriesIAExt = "IntervalArithmetic"
 
 [compat]
 Aqua = "0.8"
-InteractiveUtils = "1"
 IntervalArithmetic = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20"
 LinearAlgebra = "<0.0.1, 1"
 Markdown = "<0.0.1, 1"
diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl
index 4c3988a8..440fad9b 100644
--- a/src/TaylorSeries.jl
+++ b/src/TaylorSeries.jl
@@ -16,7 +16,6 @@ see also [`HomogeneousPolynomial`](@ref).
 """
 module TaylorSeries
 
-using InteractiveUtils
 using SparseArrays: SparseMatrixCSC
 using Markdown
 
diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index c6b2eda3..70585972 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -665,21 +665,36 @@ function /(a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}}) where {T<:NumberNotSe
     return res
 end
 
-function /(a::S, b::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries, S<:Union{T,TaylorN{T}}}
-    R = promote_type(S, TaylorN{T})
-    iszero(a) && !iszero(b) && return convert(Taylor1{R}, zero(b))
+function /(a::S, b::Taylor1{TaylorN{T}}) where {S<:NumberNotSeries, T<:NumberNotSeries}
+    R = promote_type(TaylorN{S}, TaylorN{T})
+    res = convert(Taylor1{R}, zero(b))
+    iszero(a) && !iszero(b) && return res
 
     # 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]
 
-    res = Taylor1(a/b[0], b.order)
     for ordT in eachindex(res)
         div!(res, a, b, ordT)
     end
     return res
 end
 
+function /(a::TaylorN{T}, b::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
+    res = zero(b)
+    iszero(a) && !iszero(b) && return res
+
+    # 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]
+
+    aa = Taylor1(a, b.order)
+    for ordT in eachindex(res)
+        div!(res, aa, b, ordT)
+    end
+    return res
+end
+
 @inline function divfactorization(a1::Taylor1, b1::Taylor1)
     # order of first factorized term; a1 and b1 assumed to be of the same order
     a1nz = findfirst(a1)
diff --git a/test/mixtures.jl b/test/mixtures.jl
index 242d946d..a55ff8ec 100644
--- a/test/mixtures.jl
+++ b/test/mixtures.jl
@@ -490,6 +490,10 @@ using Test
     @test b.coeffs[2] === b.coeffs[3]
     b.coeffs[2:end] .= zero.(b.coeffs[1])
     @test !(b.coeffs[2] === b.coeffs[3])
+    x = Taylor1([1.0+X,-X, Y, X-Y,X])
+    z = zero(x)
+    two = 2one(x[0])
+    @test two/x == 2/x == 2.0/x
 end
 
 @testset "Tests with nested Taylor1s" begin

From f07594f8121c4630ce434ed974283eeaaf64ee3e 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: Fri, 16 Feb 2024 23:22:09 +0100
Subject: [PATCH 19/23] Cleanup and add test

---
 src/arithmetic.jl | 26 --------------------------
 test/mixtures.jl  |  1 +
 2 files changed, 1 insertion(+), 26 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 70585972..ffc4364e 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -796,7 +796,6 @@ end
         c[k] += c[i] * b[k-i]
     end
     @inbounds c[k] = -c[k]/b[0]
-    # @inbounds divsubst!(c[k], b[0])
     return nothing
 end
 
@@ -852,31 +851,6 @@ function mul!(c::TaylorN, a::TaylorN, b::TaylorN)
     end
 end
 
-# # inplace method for computing c = -c/b
-# # NOTE: Here `divsubst!` *accumulates* the result inplace
-# @inline function divsubst!(c::TaylorN, b::TaylorN, k::Int)
-#     if k==0
-#         @inbounds for i in eachindex(c[0])
-#             c[0][i] = -c[0][i] / constant_term(b)
-#         end
-#         return nothing
-#     end
-
-#     @inbounds for i = 0:k-1
-#         mul!(c[k], c[i], b[k-i])
-#     end
-#     @inbounds for i in eachindex(c[k])
-#         c[k][i] = (-2c[k][i]) / constant_term(b)
-#     end
-#     return nothing
-# end
-
-# function divsubst!(c::TaylorN, b::TaylorN)
-#     for k in eachindex(c)
-#         divsubst!(c, b, k)
-#     end
-# end
-
 # NOTE: Here `div!` *accumulates* the result of a / b in res[k] (k > 0)
 @inline function div!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
         b::Taylor1{TaylorN{T}}, ordT::Int) where {T<:NumberNotSeriesN}
diff --git a/test/mixtures.jl b/test/mixtures.jl
index a55ff8ec..2cb0078b 100644
--- a/test/mixtures.jl
+++ b/test/mixtures.jl
@@ -494,6 +494,7 @@ using Test
     z = zero(x)
     two = 2one(x[0])
     @test two/x == 2/x == 2.0/x
+    @test (2one(x))/x == 2/x
 end
 
 @testset "Tests with nested Taylor1s" begin

From 20c20d02faa6edf516147b7c40d9cb13881202b7 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: Tue, 20 Feb 2024 17:35:34 +0100
Subject: [PATCH 20/23] Cleanup

---
 src/arithmetic.jl |  33 ++++++-----
 src/functions.jl  | 139 +++++++++++++++++++++++-----------------------
 2 files changed, 88 insertions(+), 84 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index ffc4364e..e3d16a3a 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -491,6 +491,20 @@ function mul!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, b::Taylor1{Taylo
     return nothing
 end
 
+@inline 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]
+    end
+end
+return nothing
+end
+
+mul!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
+b::NumberNotSeries, k::Int) where {T<:NumberNotSeries} = mul!(res, b, a, k)
+
+
 
 @doc doc"""
     mul!(c, a, b, k::Int) --> nothing
@@ -779,9 +793,13 @@ end
     return nothing
 end
 
+# TODO: get rid of remaining allocations here.
+#       This can either be achieved via pre-allocating auxiliary variables
+#       with can be passed here as arguments, or adding allocation-free in-place
+#       methods for operations such as `a += a*b` and `a = -a/b`, where a and b are
+#       TaylorN variables. See #347 for further discussion.
 @inline function div!(c::Taylor1{TaylorN{T}}, a::NumberNotSeries,
         b::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries}
-    ### TODO: use non-allocating div!(::TaylorN,::TaylorN,::TaylorN)
     iszero(a) && !iszero(b) && zero!(c, k)
     # order and coefficient of first factorized term
     # In this case, since a[k]=0 for k>0, we can simplify to:
@@ -896,19 +914,6 @@ end
     return nothing
 end
 
-@inline 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]
-        end
-    end
-    return nothing
-end
-
-mul!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
-    b::NumberNotSeries, k::Int) where {T<:NumberNotSeries} = mul!(res, b, a, k)
-
 
 
 """
diff --git a/src/functions.jl b/src/functions.jl
index db1e6540..87f5ea9c 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -234,6 +234,75 @@ end
 
 
 # Recursive functions (homogeneous coefficients)
+@inline function zero!(a::Taylor1{T}, k::Int) where {T<:NumberNotSeries}
+    a[k] = zero(a[k])
+    return nothing
+end
+
+@inline function zero!(a::Taylor1{T}) where {T<:NumberNotSeries}
+    for k in eachindex(a)
+        zero!(a, k)
+    end
+    return nothing
+end
+
+@inline function zero!(a::HomogeneousPolynomial{T}, k::Int) where {T<:NumberNotSeries}
+    a[k] = zero(a[k])
+    return nothing
+end
+
+@inline function zero!(a::HomogeneousPolynomial{T}) where {T<:NumberNotSeries}
+    for k in eachindex(a)
+        zero!(a, k)
+    end
+    return nothing
+end
+
+@inline function zero!(a::TaylorN{T}, k::Int) where {T<:NumberNotSeries}
+    zero!(a[k])
+    return nothing
+end
+
+@inline function zero!(a::TaylorN{T}) where {T<:NumberNotSeries}
+    for k in eachindex(a)
+        zero!(a, k)
+    end
+    return nothing
+end
+
+@inline function zero!(a::Taylor1{Taylor1{T}}, k::Int) where {T<:NumberNotSeries}
+    for l in eachindex(a[k])
+        zero!(a[k], l)
+    end
+    return nothing
+end
+
+@inline function zero!(a::Taylor1{Taylor1{T}}) where {T<:NumberNotSeries}
+    for k in eachindex(a)
+        zero!(a, k)
+    end
+    return nothing
+end
+
+@inline function zero!(a::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries}
+    zero!(a[k])
+    return nothing
+end
+
+@inline function zero!(a::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
+    for k in eachindex(a)
+        zero!(a, k)
+    end
+    return nothing
+end
+
+@inline function zero!(a::TaylorN{Taylor1{T}}, k::Int) where {T<:NumberNotSeries}
+    for l in eachindex(a[k])
+        zero!(a[k][l])
+    end
+    return nothing
+end
+
 for T in (:Taylor1, :TaylorN)
     @eval begin
         @inline function identity!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
@@ -242,7 +311,6 @@ for T in (:Taylor1, :TaylorN)
         end
 
         @inline function zero!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
-            # @inbounds c[k] = zero(a[k])
             @inbounds zero!(c, k)
             return nothing
         end
@@ -656,75 +724,6 @@ end
 
 
 # Mutating functions for Taylor1{TaylorN{T}}
-@inline function zero!(a::Taylor1{T}, k::Int) where {T<:NumberNotSeries}
-    a[k] = zero(a[k])
-    return nothing
-end
-
-@inline function zero!(a::Taylor1{Taylor1{T}}, k::Int) where {T<:NumberNotSeries}
-    for l in eachindex(a[k])
-        zero!(a[k], l)
-    end
-    return nothing
-end
-
-@inline function zero!(a::Taylor1{Taylor1{T}}) where {T<:NumberNotSeries}
-    for k in eachindex(a)
-        zero!(a, k)
-    end
-    return nothing
-end
-
-@inline function zero!(a::Taylor1{T}) where {T<:NumberNotSeries}
-    for k in eachindex(a)
-        zero!(a, k)
-    end
-    return nothing
-end
-
-@inline function zero!(a::HomogeneousPolynomial{T}, k::Int) where {T<:NumberNotSeries}
-    a[k] = zero(a[k])
-    return nothing
-end
-
-@inline function zero!(a::HomogeneousPolynomial{T}) where {T<:NumberNotSeries}
-    for k in eachindex(a)
-        zero!(a, k)
-    end
-    return nothing
-end
-
-@inline function zero!(a::TaylorN{T}, k::Int) where {T<:NumberNotSeries}
-    zero!(a[k])
-    return nothing
-end
-
-@inline function zero!(a::TaylorN{T}) where {T<:NumberNotSeries}
-    for k in eachindex(a)
-        zero!(a, k)
-    end
-    return nothing
-end
-
-@inline function zero!(a::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries}
-    zero!(a[k])
-    return nothing
-end
-
-@inline function zero!(a::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
-    for k in eachindex(a)
-        zero!(a, k)
-    end
-    return nothing
-end
-
-@inline function zero!(a::TaylorN{Taylor1{T}}, k::Int) where {T<:NumberNotSeries}
-    for l in eachindex(a[k])
-        zero!(a[k][l])
-    end
-    return nothing
-end
-
 @inline function exp!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
         k::Int) where {T<:NumberNotSeries}
     if k == 0

From a5fceb8dddfacead50d3b8fd5badd57d4a1f5677 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: Tue, 20 Feb 2024 17:41:46 +0100
Subject: [PATCH 21/23] More cleanup

---
 src/TaylorSeries.jl |  1 +
 src/arithmetic.jl   | 14 +++++++-------
 2 files changed, 8 insertions(+), 7 deletions(-)

diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl
index 440fad9b..dca66013 100644
--- a/src/TaylorSeries.jl
+++ b/src/TaylorSeries.jl
@@ -16,6 +16,7 @@ see also [`HomogeneousPolynomial`](@ref).
 """
 module TaylorSeries
 
+
 using SparseArrays: SparseMatrixCSC
 using Markdown
 
diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index e3d16a3a..3969592b 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -504,6 +504,13 @@ end
 mul!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
 b::NumberNotSeries, k::Int) where {T<:NumberNotSeries} = mul!(res, b, a, k)
 
+# in-place division (assumes equal order among TaylorNs)
+function mul!(c::TaylorN, a::TaylorN, b::TaylorN)
+    for k in eachindex(c)
+        mul!(c, a, b, k)
+    end
+end
+
 
 
 @doc doc"""
@@ -862,13 +869,6 @@ function div!(c::TaylorN, a::NumberNotSeries, b::TaylorN)
     end
 end
 
-# inplace multiplication TaylorNs
-function mul!(c::TaylorN, a::TaylorN, b::TaylorN)
-    for k in eachindex(c)
-        mul!(c, a, b, k)
-    end
-end
-
 # NOTE: Here `div!` *accumulates* the result of a / b in res[k] (k > 0)
 @inline function div!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
         b::Taylor1{TaylorN{T}}, ordT::Int) where {T<:NumberNotSeriesN}

From 21b99b7599dcfd42d4e22040d4253ae196aff7f5 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: Tue, 20 Feb 2024 17:50:28 +0100
Subject: [PATCH 22/23] Remove comments

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

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 3969592b..bb22a986 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -691,10 +691,6 @@ function /(a::S, b::Taylor1{TaylorN{T}}) where {S<:NumberNotSeries, T<:NumberNot
     res = convert(Taylor1{R}, zero(b))
     iszero(a) && !iszero(b) && return res
 
-    # 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]
-
     for ordT in eachindex(res)
         div!(res, a, b, ordT)
     end
@@ -705,10 +701,6 @@ function /(a::TaylorN{T}, b::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
     res = zero(b)
     iszero(a) && !iszero(b) && return res
 
-    # 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]
-
     aa = Taylor1(a, b.order)
     for ordT in eachindex(res)
         div!(res, aa, b, ordT)

From 40b5267b831197e75ba1e456986b7929868d03e7 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: Fri, 23 Feb 2024 23:19:48 +0100
Subject: [PATCH 23/23] Update Project.toml: bump minor version and use hyphen
 instead of commas in IntervalArithmetic compat entry

---
 Project.toml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index 134ce35a..f2d4fd2c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "TaylorSeries"
 uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
-version = "0.16.0"
+version = "0.17.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -17,7 +17,7 @@ TaylorSeriesIAExt = "IntervalArithmetic"
 
 [compat]
 Aqua = "0.8"
-IntervalArithmetic = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20"
+IntervalArithmetic = "0.15 - 0.20"
 LinearAlgebra = "<0.0.1, 1"
 Markdown = "<0.0.1, 1"
 Requires = "0.5.2, 1"