Enable NaN-safe mode by default #857

I'd like to propose enabling NaN-safe mode by default, and making the current behaviour opt-in via the nansafe_mode preference.

This was discussed in #719 (starting at #719 (comment)), but since switching the default is not breaking, the question was not resolved there and #719 was closed with the 1.0 release. So I open a separate issue for it, as suggested in #719 (comment).

IMO NaN-safe mode is not primarily about NaNs but about the definition of dual numbers: a dual number with a zero partial is unperturbed in that direction, so it must not contribute to the derivative in that direction, regardless of the derivative of the primal function at that point. The default mode violates this whenever the derivative is not finite (0 * Inf and 0 * NaN are NaN). Our documentation even says so ("oops, this NaN should be 0.0").

As a consequence, the default mode returns incorrect derivatives:

julia> f(v) = sqrt(v[1]) + v[2];

julia> ForwardDiff.gradient(f, [0.0, 1.0])  # ∂f/∂v[2] = 1
2-element Vector{Float64}:
 Inf
 NaN

Moreover, results depend on implementation details that should be irrelevant:

• whether an unperturbed argument is a Real or a Dual with zero partials (see Make x^y with Dual arguments with zero partials consistent with Real arguments #856),
• the chunk size: with the existing all-or-nothing check in the Dual^Real branch of ^, v[1]^0.5 + v[2] returns [Inf, 1.0] with Chunk{1} and [Inf, NaN] with Chunk{2}.

In NaN-safe mode, all of these examples return [Inf, 1.0] for any chunk size (since #777 the check is applied to each partial separately).

IMO NaN-safe mode should be the default regardless of its performance impact. Correctness should be the default, and faster but incorrect results should be opt-in (as with, e.g., @fastmath) - one can always be fast by returning nonsense. The documentation states a performance impact of ~5-10%, but that's outdated: in the Rosenbrock benchmark in #777 the gradient took 853 μs with NaN-safe mode vs 597 μs without (median), i.e., about 40% more. That's worth documenting for users who want to opt out, but it shouldn't decide the default. Also note that since #777 NaN-safe mode does not turn the partials multiplication into branchy code anymore (one of the concerns in #719 (comment)): _mul_partials compiles to SIMD instructions with an additional comparison and mask per partial (see the native code in #777).