Use dualpart to determine value on branch cut?

[dlfivefifty]

I'm proposing that functions with branch cuts be overriden to use the dualpart to determine direction of approach, that is realpart(f(dual(a,b))) should always return lim_h->0 f(a+h*b).

For example,

    f(z)=sqrt(z)
    realpart(f(dual(-1.,im)))     # correctly returns 1.im
    realpart(f(dual(-1.,-im)))     # should return -1.im but returns 1.im

While one can accomplish this example using -1.+0.im and -1.-0.im,this breaks if the branch cut is rotated, where using dual numbers to determine the angle could still work:

    f(z)=sqrt((1.+1.im)*z)
    z=-1/(1.+1.im)
    realpart(f(dual(z,-1.im)))  # would return -1im which is lim_{h->0} f(z-h*im)

[dlfivefifty]

@MikaelSlevinsky

[dlfivefifty]

One issue that might make this too difficult is that there would need to be a tolerance to determine whether one is on the branch cut or not.

[MikaelSlevinsky]

It's only sensible to take the limit if the limiting value is on the same branch, so yeah it makes sense to me!

I suppose we should look at all the overridden functions and see which have branch cuts. Branch cuts in log and sqrt (and cbrt) imply that all the inverse trig and hyperbolic functions have branch cuts. The erf and airy families are entire, gamma has poles that lgamma inherits as branch points, and besselj0/besselj1 are entire but bessely0/bessely1 have branch cuts along (-infty,0].

Not considering corner cases, it could be as simple as:

julia> using DualNumbers

julia> function Base.sqrt{T<:Complex}(z::Dual{T})
           d = Dual(sqrt(realpart(z)),epsilon(z)*(1/2/sqrt(realpart(z))))
           if imag(realpart(z)) == 0 && real(realpart(z)) < 0
               return d*sign(imag(epsilon(z)))
           else
               return d
           end
       end
sqrt (generic function with 17 methods)

julia> zd = dual(-1.0,0.0+1.0im)
-1.0 + 0.0im + 0.0ɛ + 1.0imɛ

julia> sqrt(zd)
0.0 + 1.0im + 0.5ɛ + 0.0imɛ

julia> zd = dual(-1.0,0.0-1.0im)
-1.0 + 0.0im + 0.0ɛ - 1.0imɛ

julia> sqrt(zd)
-0.0 - 1.0im + 0.5ɛ + 0.0imɛ

I think tolerances are for other methods that would depend on the DualNumbers package, but I don't feel too strongly either way about this.

[MikaelSlevinsky]

EDIT: fixed the redefinition so the dual part is positive on both sides