max when values agree?

The partial derivative of max currently uses the Julia Base implementation:

max(x,y) = ifelse(y < x, x, y)

Unfortunately, this may cause incorrect results when the Dual values agree but the partials do not:

max(Dual(3,1), Dual(3,-1))

Here, the values are both 3 but the partials have opposite signs. I think the partial here may be undefined.

Maybe something like:

function Base.max(a::Dual, b::Dual)
    if value(a) > value(b)
        return a
    elseif value(a) < value(b)
        return b
    else # value(a) == value(b)
        return ifelse(epsilon(a) == epsilon(b), a, Dual(value(a), NaN))
    end
end

Begin with the Base.max behavior described in issue #53 and reproduce max(Dual(3,1), Dual(3,-1)). Before changing code, establish the expected result for equal values with equal and unequal partials; done means the chosen behavior is covered by tests.