More minor fixes reducing use of memory

Author: lbenet

src/arithmetic.jl

@@ -163,8 +163,8 @@ For many variables, the ordering includes a lexicographical convention in order
 total. We have opted for the simplest one: the *larger* variable appears *before*
 for the `TaylorN` variables are defined (e.g., through [`set_variables`](@ref)).
 
-For nested `Taylor1`, i.e. `Taylor1{Taylor1{...}}`s, the variables are compared
-wrt to the output of `get_numvars`
+For nested `Taylor1{Taylor1{...}}`s, the ordering is established by which one
+is a `constant_term` of the other.
 
 Refs:
 - M. Berz, AIP Conference Proceedings 177, 275 (1988); https://doi.org/10.1063/1.37800
@@ -749,7 +749,7 @@ 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)
-    # c[k] = scalar * c[k]
+    # c[k] <- scalar * c[k]
     for ord in eachindex(c[k])
         mul!(c[k], c[k], scalar, ord)
     end
@@ -781,7 +781,7 @@ function mul!(a::TaylorN{T}, b::TaylorN{T}) where {T<:Number}
     end
     return nothing
 end
-function mul!(a::Taylor1{T}, b::Taylor1{T}) where {T<:NumberNotSeriesN}
+function mul!(a::Taylor1{T}, b::Taylor1{T}) where {T<:NumberNotSeries}
     @inbounds for k in reverse(eachindex(a))
         # a[k] <- a[k]*b[0]
         mul!(a, a, b[0], k)
@@ -1218,9 +1218,9 @@ function div!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
     zero!(c, k)
     iszero(a) && !iszero(b) && return nothing
     # order and coefficient of first factorized term
-    ordfact, cdivfact = divfactorization(a, b)
+    ordfact, aux = divfactorization(a, b)
     if k == 0
-        identity!(c, cdivfact, k)
+        identity!(c, aux, k)
         return nothing
     end
     imin = max(0, k+ordfact-b.order)
@@ -1232,20 +1232,22 @@ function div!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
             muladd!(c[k], c[i], b[k+ordfact-i], ord)
         end
     end
-    aux = zero(c[k])
+    zero!(aux)
     if k+ordfact ≤ b.order
-        # @inbounds aux = a[k+ordfact] - c[k]
+        # @inbounds aux <- a[k+ordfact] - c[k]
         for ord in eachindex(minlength(aux, a[k+ordfact]))
             subst!(aux, a[k+ordfact], c[k], ord)
         end
     else
-        # @inbounds aux = - c[k]
+        # @inbounds aux <- - c[k]
         for ord in eachindex(minlength(aux, a[k+ordfact]))
             subst!(aux, c[k], ord)
         end
     end
-    # TODO: optimize next line
-    c[k] = aux / b[ordfact]
+    # c[k] <- aux / b[ordfact]
+    for ord in eachindex(c[k])
+        div!(c[k], aux, b[ordfact], ord)
+    end
     return nothing
 end
 

src/functions.jl

@@ -380,7 +380,13 @@ end
             zero!(aux, j)
         end
     end
-    div!(c, c, k, k)
+    # div!(c, c, k, k)
+    for i in eachindex(c[k])
+        identity!(aux, c[k], i)
+    end
+    for i in eachindex(c[k])
+        div!(c[k], aux, k, i)
+    end
     return nothing
 end
 

src/power.jl

@@ -360,14 +360,14 @@ end
         end
     end
     # @inbounds c[k] = c[k] / (kprime * a[l0])
-    aux2 = zero(c[k])
-    for j in eachindex(a[l0])
-        zero!(aux[k], j)
-        mul!(aux[k], kprime, a[l0], j)
-        identity!(aux2, c[k], j)
+    @inbounds for j in eachindex(c[k])
+        identity!(aux[k], c[k], j)
     end
-    for j in eachindex(a[l0])
-        div!(c[k], aux2, aux[k], j)
+    @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
@@ -432,8 +432,8 @@ end
     end
     return nothing
 end
-@inline function sqr_orderzero!(c::Taylor1{T}, a::Taylor1{T}) where
-        {T<:AbstractSeries}
+@inline function sqr_orderzero!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}) where
+        {T<:Number}
     @inbounds for ord in eachindex(c[0])
         sqr!(c[0], a[0], ord)
     end
@@ -553,8 +553,8 @@ end
     return nothing
 end
 
-@inline function sqr!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, k::Int) where
-        {T<:NumberNotSeriesN}
+@inline function sqr!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
+        k::Int) where {T<:NumberNotSeriesN}
     if k == 0
         sqr_orderzero!(c, a)
         return nothing
@@ -702,7 +702,8 @@ coefficient, which must be even.
 
 """ sqrt!
 
-@inline function sqrt!(c::Taylor1{T}, a::Taylor1{T}, k::Int, k0::Int=0) where {T<:Number}
+@inline function sqrt!(c::Taylor1{T}, a::Taylor1{T}, k::Int, k0::Int=0) where
+        {T<:NumberNotSeries}
     k < k0 && return nothing
     if k == k0
         @inbounds c[k] = sqrt(a[2*k0])
@@ -728,7 +729,8 @@ coefficient, which must be even.
     return nothing
 end
 
-@inline function sqrt!(c::TaylorN{T}, a::TaylorN{T}, k::Int) where {T<:NumberNotSeriesN}
+@inline function sqrt!(c::TaylorN{T}, a::TaylorN{T}, k::Int) where
+        {T<:NumberNotSeriesN}
 
     if k == 0
         @inbounds c[0][1] = sqrt( constant_term(a) )
@@ -800,7 +802,8 @@ end
     return nothing
 end
 
-@inline function sqrt!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, k::Int, k0::Int=0) where {T<:Number}
+@inline function sqrt!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, k::Int,
+        k0::Int=0) where {T<:Number}
     k < k0 && return nothing
     if k == k0
         @inbounds c[k] = sqrt(a[2*k0])
@@ -834,10 +837,14 @@ end
         end
     end
     # @inbounds c[k] = c[k] / (2*c[k0])
+    @inbounds for j in eachindex(c[k])
+        identity!(aux, c[k], j)
+    end
     @inbounds for j in eachindex(c[k0])
-        zero!(aux, j)
-        mul!(aux, 2, c[k0], j)
+        div!(c[k], aux, c[k0], j)
+    end
+    @inbounds for j in eachindex(c[k0])
+        div!(c[k], c[k], 2, j)
     end
-    c[k] = c[k] / aux
     return nothing
 end