Correct the returned series when factorizing, with tests

.travis.yml

@@ -11,6 +11,10 @@ julia:
 notifications:
     email: false
 
+matrix:
+  allow_failures:
+  - julia: nightly
+
 script:
  - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
  - julia --check-bounds=yes -e 'Pkg.clone(pwd()); Pkg.build("TaylorSeries"); Pkg.test("TaylorSeries"; coverage=true)'

src/arithmetic.jl

@@ -414,7 +414,7 @@ function /(a::Taylor1{T}, b::Taylor1{T}) where {T<:Number}
     ordfact, cdivfact = divfactorization(a, b)
 
     c = Taylor1(cdivfact, a.order)
-    for ord = 1:a.order-ordfact
+    for ord = 1:a.order
         div!(c, a, b, ord, ordfact) # updates c[ord]
     end
 
@@ -486,9 +486,14 @@ term of the denominator.
     end
 
     @inbounds for i = 0:k-1
+        k+ordfact-i > b.order && continue
         c[k] += c[i] * b[k+ordfact-i]
     end
-    @inbounds c[k] = (a[k+ordfact]-c[k]) / b[ordfact]
+    if k+ordfact ≤ b.order
+        @inbounds c[k] = (a[k+ordfact]-c[k]) / b[ordfact]
+    else
+        @inbounds c[k] = - c[k] / b[ordfact]
+    end
     return nothing
 end
 

src/power.jl

@@ -88,7 +88,7 @@ function ^(a::Taylor1, r::S) where {S<:Real}
     k0 = lnull+l0nz
     c = Taylor1( zero(aux), a.order)
     @inbounds c[lnull] = aux
-    for k = k0+1:a.order
+    for k = k0+1:a.order+l0nz
         pow!(c, a, r, k, l0nz)
     end
 
@@ -149,6 +149,7 @@ coefficient of `a`.
     end
 
     for i = 0:k-l0-1
+        k-i > a.order && continue
         aux = r*(k-i) - i
         @inbounds c[k-l0] += aux * a[k-i] * c[i]
     end
@@ -322,7 +323,7 @@ function sqrt(a::Taylor1)
 
     c = Taylor1( zero(T), a.order )
     @inbounds c[lnull] = aux
-    for k = lnull+1:a.order-l0nz
+    for k = lnull+1:a.order
         sqrt!(c, a, k, lnull)
     end
 
@@ -377,9 +378,14 @@ coefficient, which must be even.
     kodd = (k - k0)%2
     kend = div(k - k0 - 2 + kodd, 2)
     @inbounds for i = k0+1:k0+kend
+        (k+k0-i > a.order) || (i > a.order) && continue
         c[k] += c[i] * c[k+k0-i]
     end
-    @inbounds aux = a[k+k0] - 2*c[k]
+    if k+k0 ≤ a.order
+        @inbounds aux = a[k+k0] - 2*c[k]
+    else
+        @inbounds aux = - 2*c[k]
+    end
     if kodd == 0
         @inbounds aux = aux - (c[kend+k0+1])^2
     end

test/onevariable.jl

@@ -144,6 +144,13 @@ end
     @test Taylor1(BigFloat,5)/(6*Taylor1(3)) == 1/BigInt(6)
     @test Taylor1(BigFloat,5)/(6im*Taylor1(3)) == -1im/BigInt(6)
 
+    # These tests involve some sort of factorization
+    @test t/(t+t^2) == 1/(1+t)
+    @test sqrt(t^2+t^3) == t*sqrt(1+t)
+    @test (t^3+t^4)^(1/3) ≈ t*(1+t)^(1/3)
+    @test norm((t^3+t^4)^(1/3) - t*(1+t)^(1/3), Inf) < eps()
+    @test ((t^3+t^4)^(1/3))[15] ≈ -8617640/1162261467
+
     trational = ta(0//1)
     @inferred ta(0//1) == Taylor1{Rational{Int}}
     @test eltype(trational) == Rational{Int}