Division is not properly defined for some cases of Taylor1{TaylorN{Float64}}

[PerezHz]

I'm working with Taylor1{TaylorN{Float64}} variables in Julia 0.5 using TaylorSeries latest master, but are having some trouble with two cases for the division:

• Division of two Taylor1{TaylorN{Float64}}'s
• Division of a Taylor1{TaylorN{Float64}} by a Float64

julia> using TaylorSeries

julia> x = 1.0 + TaylorN(Float64,1,order=5)
 1.0 + 1.0 x₁ + 𝒪(‖x‖⁶)

julia> y = Taylor1(x, 16)
  1.0 + 1.0 x₁ + 𝒪(‖x‖⁶) + 𝒪(t¹⁷)

julia> y/y # division of Taylor1{TaylorN{Float64}}'s
ERROR: MethodError: no method matching isinf(::TaylorSeries.TaylorN{Float64})
Closest candidates are:
  isinf(::BigFloat) at mpfr.jl:792
  isinf(::Float16) at float16.jl:118
  isinf(::Real) at float.jl:362
  ...
 in divfactorization(::TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}}, ::TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}}) at /Users/Jorge/.julia/v0.5/TaylorSeries/src/arithmetic.jl:430
 in /(::TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}}, ::TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}}) at /Users/Jorge/.julia/v0.5/TaylorSeries/src/arithmetic.jl:382

julia> y/1.0 # divide a Taylor1{TaylorN{Float64}} by a Float64
Segmentation fault: 11

It seems that in order to properly define division of two Taylor1{TaylorN{Float64}}'s, either a method isinf(::TaylorSeries.TaylorN{Float64}) or divfactorization{T<:Real}(a1::Taylor1{TaylorN{T}}, b1::Taylor1{TaylorN{T}}) must be defined. Also, in order to define the division of a Taylor1{TaylorN{Float64}} by a Float64, a specialized method /{T<:Real}(a::Taylor1{TaylorN{T}}, b::T) must be defined. Currently, what I'm doing to work around this is:

using TaylorSeries

import Base./

#specialized method of /
/{T<:Real}(a::Taylor1{TaylorN{T}}, b::T) = a * inv(b)

import TaylorSeries.divfactorization

#unsafe, provisional specialized method of divfactorization
#in order to work, isnan, isinf must be defined for TaylorN's
#modified from the currently implemented method in TaylorSeries
function divfactorization{T<:Real}(a1::Taylor1{TaylorN{T}}, b1::Taylor1{TaylorN{T}})
    # order of first factorized term; a1 and b1 assumed to be of the same order
    a1nz = TaylorSeries.findfirst(a1)
    b1nz = TaylorSeries.findfirst(b1)
    orddivfact = min(a1nz, b1nz)
    if orddivfact < 0
        orddivfact = a1.order
    end
    cdivfact = a1.coeffs[orddivfact+1] / b1.coeffs[orddivfact+1]

    return orddivfact, cdivfact
end

And then (although the divfactorization method lacks some checks), everything works smoothly:

julia> x = 1.0 + TaylorN(Float64,1,order=5)
 1.0 + 1.0 x₁ + 𝒪(‖x‖⁶)

julia> y = Taylor1(x, 16)
  1.0 + 1.0 x₁ + 𝒪(‖x‖⁶) + 𝒪(t¹⁷)

julia> y/1.0
  1.0 + 1.0 x₁ + 𝒪(‖x‖⁶) + 𝒪(t¹⁷)

julia> y/y
  1.0 + 𝒪(‖x‖⁶) + 𝒪(t¹⁷)

With this workaround, some more elaborate cases work. As stated above, it seems that the division of two Taylor1{TaylorN{Float64}}'s may also be fixed by defining isinf(::TaylorSeries.TaylorN{Float64}).

[lbenet]

A solution for this is proposed in #93; Once I get the green lights I'll merge.

[PerezHz]

Thanks for looking into this! I checked out the isinf_isnan branch locally; and got the following results:

julia> using TaylorSeries

julia> x = 1.0 + TaylorN(Float64,1,order=5)
 1.0 + 1.0 x₁ + 𝒪(‖x‖⁶)

julia> y = Taylor1(x, 16)
  1.0 + 1.0 x₁ + 𝒪(‖x‖⁶) + 𝒪(t¹⁷)

julia> y/y
  1.0 + 𝒪(‖x‖⁶) + 𝒪(t¹⁷)

julia> y/1.0
Segmentation fault: 11

So the division of two Taylor1{TaylorN{Float64}}'s is solved, but the division of a Taylor1{TaylorN{Float64}} by a Float64 still gives a segmentation fault; I think the problem with the latter is related to promotions. I would suggest adding

/{T<:Real}(a::Taylor1{TaylorN{T}}, b::T) = a * inv(b)

to solve this; what do you think?

[lbenet]

I've just pushed another commit which solves the second case. There is still a problem though, and it is indeed due to promotion. I'm checking this.

[lbenet]

I've merged #93, and I'll close this. The general problem persists; see #94