Dual for 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

Start by reproducing max(Dual(3,1), Dual(3,-1)) and reading Julia Base's max behavior described in the issue. Determine the intended result when values agree but partials differ, then add coverage for the agreed behavior and verify the existing Dual max cases still pass.