New method of div! for the division of a NumberNotSeries by a Taylor1{T} (#387)

* Allocations reduction in div!(::Taylor1{Taylor1{T}}, ::NumberNotSeries, ::Taylor1{Taylor1{T}}, ::Int)

* Allocations reduction in subst!(::Taylor1{Taylor1{T}}, ::Taylor1{Taylor1{T}}, ::Int)

* Fix to div_scalar!(::Taylor1{T}, ::NumberNotSeries, ::Taylor1{T}, ::Int) suggested by @lbenet

* Add tests for mutating functions of nested Taylor1s

* Add Random as extra dependency and fix nested Taylor1s tests

* Bump patch version

Author: LuEdRaMo

Project.toml

@@ -1,7 +1,7 @@
 name = "TaylorSeries"
 uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
-version = "0.20.4"
+version = "0.20.5"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -26,6 +26,7 @@ IntervalArithmetic = "0.23"
 JLD2 = "0.5"
 LinearAlgebra = "<0.0.1, 1"
 Markdown = "<0.0.1, 1"
+Random = "1"
 RecursiveArrayTools = "2, 3"
 SparseArrays = "<0.0.1, 1"
 StaticArrays = "1"
@@ -37,10 +38,11 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["IntervalArithmetic", "JLD2", "LinearAlgebra", "RecursiveArrayTools", "SparseArrays", "StaticArrays", "Test", "Aqua"]
+test = ["IntervalArithmetic", "JLD2", "LinearAlgebra", "RecursiveArrayTools", "SparseArrays", "StaticArrays", "Test", "Aqua", "Random"]

src/arithmetic.jl

@@ -473,6 +473,13 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
     end
 end
 
+function subst!(v::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, k::Int) where {T <: TS.NumberNotSeries}
+    @inbounds for i in eachindex(v[k])
+        v[k][i] = -a[k][i]
+    end
+    return nothing
+end
+
 for T in (:Taylor1, :TaylorN)
     @eval begin
         function sum!(v::$T{S}, a::AbstractArray{$T{S}}) where {S <: Number}
@@ -640,6 +647,15 @@ function muladd!(c::Taylor1{T}, a::Taylor1{T}, b::Taylor1{T}, k::Int) where
     end
     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)
+    end
+    return nothing
+end
+
 # function muladd!(v::Taylor1{T}, a::Taylor1{T}, b::NumberNotSeries, k::Int) where
 #         {T<:NumberNotSeries}
 #     @inbounds v[k] += a[k] * b
@@ -1194,7 +1210,7 @@ end
 
 # @inline
 function div!(c::Taylor1{T}, a::NumberNotSeries, b::Taylor1{T}, k::Int) where
-        {T<:Number}
+        {T<:NumberNotSeries}
     zero!(c, k)
     iszero(a) && !iszero(b) && return nothing
     # order and coefficient of first factorized term
@@ -1212,6 +1228,34 @@ function div!(c::Taylor1{T}, a::NumberNotSeries, b::Taylor1{T}, k::Int) where
     return nothing
 end
 
+function div!(c::Taylor1, a::NumberNotSeries, b::Taylor1)
+    @inbounds for k in eachindex(c)
+        div!(c, a, b, k)
+    end
+    return nothing
+end
+
+@inline function div!(c::Taylor1{Taylor1{T}}, a::NumberNotSeries,
+        b::Taylor1{Taylor1{T}}, k::Int) where {T<:NumberNotSeries}
+    zero!(c, k)
+    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 div!(c[0], a, b[0])
+        return nothing
+    end
+    @inbounds mul!(c[k], c[0], b[k])
+    @inbounds for i = 1:k-1
+        # c[k] += c[i] * b[k-i]
+        muladd!(c[k], c[i], b[k-i])
+    end
+    # @inbounds c[k] = -c[k]/b[0]
+    @inbounds div_scalar!(c[k], -1, b[0])
+    return nothing
+end
+
 #
 # @inline
 function div!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}},
@@ -1261,26 +1305,6 @@ end
 #     return nothing
 # end
 
-# @inline function div!(c::Taylor1{Taylor1{T}}, a::NumberNotSeries,
-#         b::Taylor1{Taylor1{T}}, k::Int) where {T<:Number}
-#     zero!(c, k)
-#     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]
-#         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}
     zero!(c, k)
@@ -1381,6 +1405,23 @@ end
     return nothing
 end
 
+@inline function div_scalar!(c::Taylor1{T}, scalar::NumberNotSeries,
+        a::Taylor1{T}, k::Int) where {T <: NumberNotSeries}
+    if k==0
+        @inbounds c[0] = scalar*constant_term(c) / constant_term(a)
+        return nothing
+    end
+
+    aux = c[k] = scalar * c[k]
+    @inbounds zero!(c, k)
+    @inbounds for i = 0:k-1
+        c[k] -= c[i] * a[k-i]
+    end
+    c[k] += aux
+    @inbounds c[k] = c[k] / constant_term(a)
+    return nothing
+end
+
 # NOTE: Here `div!` *accumulates* the result of a[k] / b[k] in c[k] (k > 0)
 @inline function div!(c::TaylorN, a::NumberNotSeries, b::TaylorN, k::Int)
     if k==0
@@ -1413,6 +1454,14 @@ function div_scalar!(c::TaylorN, scalar::NumberNotSeries, a::TaylorN)
     return nothing
 end
 
+# in-place division c <- scalar*c/a (assumes equal order among TaylorNs)
+function div_scalar!(c::Taylor1, scalar::NumberNotSeries, a::Taylor1)
+    @inbounds for k in eachindex(c)
+        div_scalar!(c, scalar, a, k)
+    end
+    return nothing
+end
+
 # c[k] <- (a/b)[k]
 function div!(c::TaylorN, a::TaylorN, b::TaylorN)
     @inbounds for k in eachindex(c)

test/mixtures.jl

@@ -636,6 +636,43 @@ end
 
     end
 
+    @testset "Tests for mutating functions of nested Taylor1s" begin
+        import Random
+        Random.seed!(1234)
+
+        for (f, fc) in ((:+, :(add!)), (:-, :(subst!)), (:*, :(mul!)), (:/, :(div!)))
+            @eval begin
+
+                local inorder = 6
+                local outorder = 25
+
+                c = rand()
+                x = Taylor1([Taylor1(rand(inorder+1), inorder) for _ in 1:outorder+1], outorder)
+                y = Taylor1([Taylor1(rand(inorder+1), inorder) for _ in 1:outorder+1], outorder)
+                z = zero(x)
+
+                w = $f(x, y)
+                for k in eachindex(z)
+                    TS.$fc(z, x, y, k)
+                end
+                @test norm(z - w, Inf) == 0.0
+
+                w = $f(c, y)
+                for k in eachindex(z)
+                    TS.$fc(z, c, y, k)
+                end
+                @test norm(z - w, Inf) == 0.0
+
+                w = $f(y, c)
+                for k in eachindex(z)
+                    TS.$fc(z, y, c, k)
+                end
+                @test norm(z - w, Inf) == 0.0
+
+            end
+        end
+    end
+
     # Back to default
     set_taylor1_varname(1, "t")
     @test TS._params_Taylor1_.var_name == ["t"]