Improve nested Taylors (#381)

* Better display of Taylor1{Taylor1{...}}, and minor fixes in pretty_print

* Fixes in promotion for nested Taylor1s, and total ordering

* Add promotion methods for nested Taylor1, and refine tests

* Add various methods involving nested Taylor1s and untangle some metaprogramming stuff

* Methods for nested Taylor1s functions

* Minor fixes in tests

* Add tests for functions on nested Taylor1s

* Minor fixes

* Add specialized methods for pow!, sqr! and sqrt!

* Comment seemingly unused (arithmetic) methods; minor fixes in abs

* Minor improvements in power

* Fix a bug in pow!

* More minor fixes reducing use of memory

* Improvements in evaluate and differentiate

Author: lbenet

src/arithmetic.jl

@@ -291,6 +291,17 @@ for T in (:HomogeneousPolynomial, :TaylorN)
 end
 
 
+## minlength
+@inline function minlength(a::Taylor1, b::Taylor1)
+    length(eachindex(a)) < length(eachindex(b)) && return a
+    return b
+end
+@inline function minlength(a::Taylor1, b::Taylor1, c::Taylor1)
+    length(minlength(a, c)) < length(eachindex(b)) && return minlength(a, c)
+    return b
+end
+
+
 ## _isthinzero
 """
     _isthinzero(x)

src/auxiliary.jl

@@ -57,15 +57,28 @@ p_k = (k+1) a_{k+1}.
 ```
 
 """
-function differentiate!(p::Taylor1, a::Taylor1, k::Int)
+function differentiate!(p::Taylor1{T}, a::Taylor1{T}, k::Int) where
+        {T<:NumberNotSeries}
     k >= a.order && return nothing
     @inbounds p[k] = (k+1)*a[k+1]
     return nothing
 end
 function differentiate!(p::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, k::Int) where
         {T<:NumberNotSeries}
     k >= a.order && return nothing
-    @inbounds p[k] = (k+1)*a[k+1]
+    @inbounds for ord in eachindex(p[k])
+        zero!(p[k], ord)
+        mul!(p[k], a[k+1], k+1, ord)
+    end
+    return nothing
+end
+function differentiate!(p::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}, k::Int) where
+        {T<:NumberNotSeriesN}
+    k >= a.order && return nothing
+    # @inbounds p[k] = (k+1)*a[k+1]
+    @inbounds for ord in eachindex(p[k])
+        mul!(p[k], a[k+1], k+1, ord)
+    end
     return nothing
 end
 

src/calculus.jl

@@ -214,3 +214,15 @@ const NumberNotSeriesN = Union{Real,Complex,Taylor1}
 ## Additional Taylor1 and TaylorN outer constructor ##
 Taylor1{T}(x::S) where {T<:Number,S<:NumberNotSeries} = Taylor1([convert(T,x)], 0)
 TaylorN{T}(x::S) where {T<:Number,S<:NumberNotSeries} = TaylorN(convert(T, x), TaylorSeries.get_order())
+
+
+# """
+#     get_numvars
+#
+# Return the number of variables of a `Taylor1`, `HomogeneousPolynomial`
+# or `TaylorN` object.
+# """
+get_numvars(t::Number) = 0
+get_numvars(t::Taylor1) = 1
+get_numvars(t::Taylor1{Taylor1{T}}) where {T<:Number} = get_numvars(t[0])+1
+get_numvars(::T) where {T<:Union{HomogeneousPolynomial, TaylorN}} = _params_TaylorN_.num_vars

src/constructors.jl

@@ -187,20 +187,41 @@ promote_rule(::Type{Taylor1{TaylorN{T}}}, ::Type{TaylorN{Taylor1{S}}}) where
 promote_rule(::Type{TaylorN{Taylor1{T}}}, ::Type{Taylor1{TaylorN{S}}}) where
     {T<:NumberNotSeries, S<:NumberNotSeries} = Taylor1{TaylorN{promote_type(T,S)}}
 
-# Nested Taylor1's
-function promote(a::Taylor1{Taylor1{T}}, b::Taylor1{T}) where {T<:NumberNotSeriesN}
-    order_a = get_order(a)
-    order_b = get_order(b)
-    zb = zero(b)
-    new_bcoeffs = similar(a.coeffs)
-    new_bcoeffs[1] = b
-    @inbounds for ind in 2:order_a+1
-        new_bcoeffs[ind] = zb
+# Different nested Taylor1's promotion methods; consider to 4 levels of nesting
+let
+    strucTT = :(Taylor1{T})
+    strucTR = :(Taylor1{R})
+    strucT = :(T)
+    strucS = :(S)
+    strucR = :(R)
+    expr_elem = :(bnew[0])
+    for j = 0:5
+        for i = 1+j:6
+            @eval function Base._promote(a::$strucTT, b::$strucS) where {T<:NumberNotSeries,
+                    S<:NumberNotSeries}
+                R = promote_type(T, S)
+                anew = convert($strucTR, a)
+                bnew = zero(anew)
+                $(expr_elem) = convert($strucR, b)
+                return anew, bnew
+            end
+            @eval Base._promote(b::$strucS, a::$strucTT) where {T<:NumberNotSeries,
+                S<:NumberNotSeries} = reverse(Base._promote(a, b))
+            strucTT = :(Taylor1{$strucTT})
+            strucTR = :(Taylor1{$strucTR})
+            expr_elem = :(($(expr_elem))[0])
+        end
+        strucS = :(Taylor1{$strucS})
+        strucR = :(Taylor1{$strucR})
+        strucT = :(Taylor1{$strucT})
+        strucTT = deepcopy(strucT)
+        strucTT = :(Taylor1{$strucT})
+        strucTR = deepcopy(strucR)
+        strucTR = :(Taylor1{$strucTR})
+        expr_elem = :(bnew[0])
     end
-    return a, Taylor1(b, order_a)
 end
-promote(b::Taylor1{T}, a::Taylor1{Taylor1{T}}) where {T<:NumberNotSeriesN} =
-    reverse(promote(a, b))
+
 
 # float
 float(::Type{Taylor1{T}}) where T<:Number = Taylor1{float(T)}