Reimplement pow! and div! for Taylor1s (#236)

* Reimplement pow! for Taylor1

The first non-zero element is computed there, instead of being passed as an argument.

* Remove carets from Project.toml, and upgrade versions of some packages

* Readd IntervalArithmetics for Julia 1.0 compatibility

* Reimplement div! in a similar way as pow!

* Bump patch version

Project.toml

@@ -1,7 +1,7 @@
 name = "TaylorSeries"
 uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
-version = "0.10.1"
+version = "0.10.2"
 
 [deps]
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
@@ -11,8 +11,8 @@ Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-IntervalArithmetic = "^0.15.0"
-Requires = "^v0.5.2"
+IntervalArithmetic = "0.15, 0.16"
+Requires = "0.5.2, 1.0"
 julia = "1"
 
 [extras]

src/arithmetic.jl

@@ -422,7 +422,7 @@ function /(a::Taylor1{T}, b::Taylor1{T}) where {T<:Number}
 
     c = Taylor1(cdivfact, a.order)
     for ord in eachindex(c)
-        div!(c, a, b, ord, ordfact) # updates c[ord]
+        div!(c, a, b, ord) # updates c[ord]
     end
 
     return c
@@ -449,7 +449,7 @@ function /(a::TaylorN{T}, b::TaylorN{T}) where {T<:NumberNotSeriesN}
 end
 
 
-function divfactorization(a1::Taylor1, b1::Taylor1)
+@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)
     b1nz = findfirst(b1)
@@ -472,7 +472,7 @@ end
 
 # Homogeneous coefficient for the division
 @doc doc"""
-    div!(c, a, b, k::Int, ordfact::Int=0)
+    div!(c, a, b, k::Int)
 
 Compute the `k-th` expansion coefficient `c[k]` of `c = a / b`,
 where all `c`, `a` and `b` are either `Taylor1` or `TaylorN`.
@@ -483,13 +483,16 @@ The coefficients are given by
 c_k =  \frac{1}{b_0} \big(a_k - \sum_{j=0}^{k-1} c_j b_{k-j}\big).
 ```
 
-For `Taylor1` polynomials, `ordfact` is the order of the factorized
-term of the denominator.
+For `Taylor1` polynomials, a similar formula is implemented which
+exploits `k_0`, the order of the first non-zero coefficient of `a`.
 """ div!
 
-@inline function div!(c::Taylor1, a::Taylor1, b::Taylor1, k::Int, ordfact::Int=0)
+@inline function div!(c::Taylor1, a::Taylor1, b::Taylor1, k::Int)
+
+    # order and coefficient of first factorized term
+    ordfact, cdivfact = divfactorization(a, b)
     if k == 0
-        @inbounds c[0] = a[ordfact] / b[ordfact]
+        @inbounds c[0] = cdivfact
         return nothing
     end
 

src/power.jl

@@ -117,22 +117,9 @@ function ^(a::Taylor1, r::S) where {S<:Real}
     r == 2 && return square(aa)
     r == 1/2 && return sqrt(aa)
 
-    # First non-zero coefficient
-    l0nz = findfirst(a)
-    l0nz < 0 && return Taylor1(zero(aux), a.order)
-
-    # The first non-zero coefficient of the result; must be integer
-    !isinteger(r*l0nz) && throw(ArgumentError(
-        """The 0th order Taylor1 coefficient must be non-zero
-        to raise the Taylor1 polynomial to a non-integer exponent."""))
-    lnull = trunc(Int, r*l0nz )
-    lnull > a.order && return Taylor1(zero(aux), a.order)
-
-    k0 = lnull+l0nz
     c = Taylor1( zero(aux), a.order)
-    @inbounds c[lnull] = ( a[l0nz] )^r
-    for k = k0+1:a.order+l0nz
-        pow!(c, aa, r, k, l0nz)
+    for k = 0:a.order
+        pow!(c, aa, r, k)
     end
 
     return c
@@ -166,7 +153,7 @@ end
 
 # Homogeneous coefficients for real power
 @doc doc"""
-    pow!(c, a, r::Real, k::Int, k0::Int=0)
+    pow!(c, a, r::Real, k::Int)
 
 Update the `k`-th expansion coefficient `c[k]` of `c = a^r`, for
 both `c` and `a` either `Taylor1` or `TaylorN`.
@@ -177,12 +164,12 @@ The coefficients are given by
 c_k = \frac{1}{k a_0} \sum_{j=0}^{k-1} \big(r(k-j) -j\big)a_{k-j} c_j.
 ```
 
-For `Taylor1` polynomials, `k0` is the order of the first non-zero
-coefficient of `a`.
+For `Taylor1` polynomials, a similar formula is implemented which
+exploits `k_0`, the order of the first non-zero coefficient of `a`.
 
 """ pow!
 
-@inline function pow!(c::Taylor1{T}, a::Taylor1{T}, r::S, k::Int, l0::Int=0) where
+@inline function pow!(c::Taylor1{T}, a::Taylor1{T}, r::S, k::Int) where
         {T<:Number, S<:Real}
 
     if r == 0
@@ -195,27 +182,47 @@ coefficient of `a`.
         return sqrt!(c, a, k)
     end
 
-    if k == l0
-        @inbounds c[0] = ( a[l0] )^r
+    # First non-zero coefficient
+    l0 = findfirst(a)
+    if l0 < 0
+        c[k] = zero(a[0])
+        return nothing
+    end
+
+    # The first non-zero coefficient of the result; must be integer
+    !isinteger(r*l0) && throw(ArgumentError(
+        """The 0th 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
+    if (kprime < 0) || (lnull > a.order)
+        @inbounds c[k] = zero(a[0])
         return nothing
     end
 
     # Relevant for positive integer r, to avoid round-off errors
-    if isinteger(r) && (k-l0 > r*findlast(a) > 0)
-        @inbounds c[k-l0] = zero(c[0])
+    if isinteger(r) && (k > r*findlast(a))
+        @inbounds c[k] = zero(a[0])
         return nothing
     end
 
-    imin = max(0,k-a.order)
-    aux = r*(k-imin) - imin
-    @inbounds c[k-l0] = aux * a[k-imin] * c[imin]
-    for i = imin+1:k-l0-1
-        aux = r*(k-i) - i
-        @inbounds c[k-l0] += aux * a[k-i] * c[i]
+    if k == lnull
+        @inbounds c[k] = (a[l0])^r
+        return nothing
     end
-    aux = k - l0*(r+1)
-    @inbounds c[k-l0] = c[k-l0] / (aux * a[l0])
 
+    # The recursion formula
+    if l0+kprime ≤ a.order
+        @inbounds c[k] = r * kprime * c[lnull] * a[l0+kprime]
+    else
+        @inbounds c[k] = zero(a[0])
+    end
+    for i = 1:k-lnull-1
+        ((i +lnull) > a.order || (l0+kprime-i > a.order)) && continue
+        aux = r*(kprime-i) - i
+        @inbounds c[k] += aux * c[i+lnull] * a[l0+kprime-i]
+    end
+    @inbounds c[k] = c[k] / (kprime * a[l0])
     return nothing
 end
 

test/onevariable.jl

@@ -319,14 +319,18 @@ eeuler = Base.MathConstants.e
     iind, cind = TaylorSeries.divfactorization(ut, ut)
     @test iind == 1
     @test cind == 1.0
-    TaylorSeries.div!(tt, ut, ut, 0, iind)
+    TaylorSeries.div!(tt, ut, ut, 0)
     @test tt[0] == cind
     TaylorSeries.div!(tt, 1+ut, 1+ut, 0)
     @test tt[0] == 1.0
     TaylorSeries.div!(tt, 1, 1+ut, 0)
     @test tt[0] == 1.0
-    TaylorSeries.pow!(tt, 1.0+t, 1.5, 0, 0)
+    TaylorSeries.pow!(tt, 1.0+t, 1.5, 0)
     @test tt[0] == 1.0
+    TaylorSeries.pow!(tt, 0.0*t, 1.5, 0)
+    @test tt[0] == 0.0
+    TaylorSeries.pow!(tt, 0.0+t, 18, 0)
+    @test tt[0] == 0.0
     TaylorSeries.pow!(tt, 1.0+t, 1.5, 0)
     @test tt[0] == 1.0
     TaylorSeries.pow!(tt, 1.0+t, 0.5, 1)