fix some floating-point logic

[MikaelSlevinsky]

According to the Julia Base docs, == for floating-point types should fail for NaNs but not for zeros with opposite signs and isequal should have the opposite behaviour for these cases.

Currently, that behaviour is also replicated by complex numbers but not for dual numbers. To me, that seems inconsistent hence the pull request.

The corresponding change is made to isnan, isinf, and isfinite so that isnan(x) == true and isfinite(x) == false iff x != x, closing https://github.com/JuliaDiff/DualNumbers.jl/issues/15.

[KristofferC]

I think the goal is to make the dual numbers take the same code path as the a real with the same real value would.

[MikaelSlevinsky]

Ya but sometimes isfinite is used to check whether an algorithm should exit a loop (because NaN math is so much more expensive and further arithmetic may be irrelevant if returning NaNs).

[KristofferC]

In my view, it would be surprising if you got a different number of iterations until the algorithm exits when using reals or dual numbers.

[MikaelSlevinsky]

Dual numbers carry derivative information, so it's not inconceivable in a Maclaurin series that a derivative might overflow/underflow/NaN before the function does.

This PR pulls dual numbers in line with complex numbers.

[MikaelSlevinsky]

This also fixes inv and / for "large" floating-point numbers.

[MikaelSlevinsky]

Why is this better than the status quo? Well, what Base function currently allows me to check whether all components of a dual number are finite?

To be sure, the Number interface (if there is one) is extremely narrow, since the set of complex numbers with modulus 2 cannot support one or zero (also the set of integers shifted by a half).

However, not having a finite check from Base on field extensions prevents generic code from working across all such extensions. Of course it's possible to adapt an algorithm for field extensions which need a bespoke finite check, but that defeats the purpose of generic code. I want to write an algorithm, test it with complex numbers (as part of Base) as a field extension, and have it Just Work^TM for any other field extension that a user loads.

[dlfivefifty]

Seems like a big change, but feels "right": for all ε we have x + ε ≠ x + 2ε, even as ε -> 0:

julia> 1.0 == nextfloat(1.0)
false

So why would the current behaviour dual(x,ε) == dual(x,2ε) make sense?

Though isinf(dual(1.0,∞)) == true is questionable as ε*∞ is sort of like 0*∞ so not well-defined. Note:

julia> isinf(NaN)
false

so I would say the correct definition is:

isinf(x::Dual) = isinf(value(x)) & isfinite(epsilon(x))

[MikaelSlevinsky]

That's an interesting take on these checks -- ascribing the meaning of the unit dual from the limit definition. This proposal is akin to component-wise checking the matrix definition, any(isinf.([1 Inf; 0 1])), whereas the current code is a check on the diagonal.

Your argument would suggest that:

isnan(x::Dual) = isnan(value(x)) | !isfinite(epsilon(x))

since an Inf*ε would act like a NaN.

It also brings up the following issue with equality:

julia> 1.0+ Inf*ε == 1.0+Inf*ε

If Inf*ε is a NaN, then they shouldn't be equal.