Enable NaN-safe mode by default

[devmotion]

Fixes #857.

This PR enables NaN-safe mode by default. The previous behaviour can still be selected by setting the nansafe_mode preference to false. The documentation section is rewritten accordingly, and the outdated "~5%-10%" performance claim is replaced with the numbers below. Since the CI matrix already sets the preference explicitly to true and false, both modes remain tested.

With this PR, by default:

julia> ForwardDiff.gradient(v -> sqrt(v[1]) + v[2], [0.0, 1.0])
2-element Vector{Float64}:
 Inf
  1.0

julia> log(ForwardDiff.Dual{:tag}(0.0, 0.0))
Dual{:tag}(-Inf,0.0)

(previously [Inf, NaN] and Dual{:tag}(-Inf,NaN)).

I bumped the version to 1.5.0 since the default behaviour changes. Happy to revert that if you'd prefer a patch release.

Timings

Median times (Chairmarks @be, seconds=3). Each mode ran in a fresh process, and the runs alternated twice (false, true, false, true). Both runs agree within a few percent, so the table shows the second run.

Benchmark	NaN-safe disabled	NaN-safe enabled	Ratio
gradient(rosenbrock, rand(1000))	532 μs	770 μs	1.45
hessian(rosenbrock, rand(10))	3.42 μs	5.35 μs	1.56
gradient(x -> sum(exp, x), rand(100))	8.03 μs	8.13 μs	1.01
jacobian(cumprod, rand(10))	624 ns	614 ns	0.98
derivative(x -> sin(x) * exp(x) / x, 1.3)	16.2 ns	16.6 ns	1.03

Overhead shows up in code that does little more than multiply and add (Rosenbrock), where the compare and select are a sizeable fraction of the work. It is within noise once transcendental functions or memory traffic dominate.

Benchmark script

using ForwardDiff, Chairmarks, Random

function rosenbrock(x::AbstractVector)
    a = one(eltype(x))
    b = 100 * a
    result = zero(eltype(x))
    for i in 1:length(x)-1
        result += (a - x[i])^2 + b*(x[i+1] - x[i]^2)^2
    end
    return result
end
sumexp(x) = sum(exp, x)
scalar(x) = sin(x) * exp(x) / x

Random.seed!(1)
x1000 = rand(1000); x100 = rand(100); x10 = rand(10)

@be ForwardDiff.gradient($rosenbrock, $x1000) seconds=3
@be ForwardDiff.gradient($sumexp, $x100) seconds=3
@be ForwardDiff.jacobian($cumprod, $x10) seconds=3
@be ForwardDiff.hessian($rosenbrock, $x10) seconds=3
@be Ref(1.3) ForwardDiff.derivative($scalar, _[]) seconds=3

versioninfo()

Julia Version 1.13.1
Commit 96ca370cf0e (2026-09-25 19:34 UTC)
Platform Info:
  OS: macOS (arm64-apple-darwin27.0.0)
  CPU: 10 × Apple M2 Pro
  WORD_SIZE: 64
  LLVM: libLLVM-20.1.8 (ORCJIT, apple-m2)
Threads: 1 default, 1 interactive, 1 GC (on 10 virtual cores)

Generated code

In NaN-safe mode, each 2-wide vector of partials gets one extra compare and one extra mask instruction (fcmeq + bic on aarch64, fcmp une + select in LLVM IR). The code remains branch-free and fully vectorized.

@code_native debuginfo=:none ForwardDiff._mul_partials(p, p, 2.0, 1.0) with p::Partials{4,Float64}

NaN-safe mode disabled:

	stp	x29, x30, [sp, #-16]!
	mov	x29, sp
	ldp	q2, q3, [x0]
	fmul.2d	v2, v2, v0[0]
	ldp	q4, q5, [x1]
	fmul.2d	v4, v4, v1[0]
	fadd.2d	v2, v2, v4
	fmul.2d	v0, v3, v0[0]
	fmul.2d	v1, v5, v1[0]
	fadd.2d	v0, v0, v1
	stp	q2, q0, [x8]
	ldp	x29, x30, [sp], #16
	ret

NaN-safe mode enabled:

	stp	x29, x30, [sp, #-16]!
	mov	x29, sp
	ldp	q2, q3, [x0]
	fmul.2d	v4, v2, v0[0]
	fcmeq.2d	v2, v2, #0.0
	bic.16b	v2, v4, v2
	ldp	q4, q5, [x1]
	fmul.2d	v6, v4, v1[0]
	fcmeq.2d	v4, v4, #0.0
	bic.16b	v4, v6, v4
	fadd.2d	v2, v2, v4
	fmul.2d	v0, v3, v0[0]
	fcmeq.2d	v3, v3, #0.0
	bic.16b	v0, v0, v3
	fmul.2d	v1, v5, v1[0]
	fcmeq.2d	v3, v5, #0.0
	bic.16b	v1, v1, v3
	fadd.2d	v0, v0, v1
	stp	q2, q0, [x8]
	ldp	x29, x30, [sp], #16
	ret

@code_llvm debuginfo=:none x * x with x::Dual{Nothing,Float64,4} (function body)

NaN-safe mode disabled:

  %1 = load <2 x double>, ptr %"x::Dual.partials_ptr", align 8
  %2 = insertelement <2 x double> poison, double %"y::Dual.unbox", i64 0
  %3 = shufflevector <2 x double> %2, <2 x double> poison, <2 x i32> zeroinitializer
  %4 = fmul <2 x double> %3, %1
  %5 = load <2 x double>, ptr %"y::Dual.partials_ptr", align 8
  %6 = insertelement <2 x double> poison, double %"x::Dual.unbox", i64 0
  %7 = shufflevector <2 x double> %6, <2 x double> poison, <2 x i32> zeroinitializer
  %8 = fmul <2 x double> %7, %5
  %9 = fadd <2 x double> %4, %8
  store <2 x double> %9, ptr %"new::Dual.sroa.2.0.sret_return.sroa_idx", align 8
  ; (same for partials 3 and 4)

NaN-safe mode enabled:

  %1 = load <2 x double>, ptr %"x::Dual.partials_ptr", align 8
  %2 = insertelement <2 x double> poison, double %"y::Dual.unbox", i64 0
  %3 = shufflevector <2 x double> %2, <2 x double> poison, <2 x i32> zeroinitializer
  %4 = fmul <2 x double> %3, %1
  %5 = fcmp une <2 x double> %1, zeroinitializer
  %6 = select <2 x i1> %5, <2 x double> %4, <2 x double> zeroinitializer
  %7 = load <2 x double>, ptr %"y::Dual.partials_ptr", align 8
  %8 = insertelement <2 x double> poison, double %"x::Dual.unbox", i64 0
  %9 = shufflevector <2 x double> %8, <2 x double> poison, <2 x i32> zeroinitializer
  %10 = fmul <2 x double> %9, %7
  %11 = fcmp une <2 x double> %7, zeroinitializer
  %12 = select <2 x i1> %11, <2 x double> %10, <2 x double> zeroinitializer
  %13 = fadd <2 x double> %6, %12
  store <2 x double> %13, ptr %"new::Dual.sroa.2.0.sret_return.sroa_idx", align 8
  ; (same for partials 3 and 4)

@code_llvm debuginfo=:none exp(x) with x::Dual{Nothing,Float64,4} (function body)

NaN-safe mode disabled:

  %0 = call double @j_exp_2207(double %"d::Dual.unbox")
  ...
  %1 = load <2 x double>, ptr %"d::Dual.partials_ptr", align 8
  %2 = insertelement <2 x double> poison, double %0, i64 0
  %3 = shufflevector <2 x double> %2, <2 x double> poison, <2 x i32> zeroinitializer
  %4 = fmul <2 x double> %3, %1
  store <2 x double> %4, ptr %"new::Dual.sroa.2.0.sret_return.sroa_idx", align 8
  ; (same for partials 3 and 4)

NaN-safe mode enabled:

  %0 = call double @j_exp_2244(double %"d::Dual.unbox")
  ...
  %1 = load <2 x double>, ptr %"d::Dual.partials_ptr", align 8
  %2 = insertelement <2 x double> poison, double %0, i64 0
  %3 = shufflevector <2 x double> %2, <2 x double> poison, <2 x i32> zeroinitializer
  %4 = fmul <2 x double> %3, %1
  %5 = fcmp une <2 x double> %1, zeroinitializer
  %6 = select <2 x i1> %5, <2 x double> %4, <2 x double> zeroinitializer
  store <2 x double> %6, ptr %"new::Dual.sroa.2.0.sret_return.sroa_idx", align 8
  ; (same for partials 3 and 4)

Tests

Pkg.test() passes locally in both modes: 9723 tests with the new default and 9641 with nansafe_mode = false.

🤖 Generated with Claude Code

[KristofferC]

Check your benchmarks with #859 perhaps

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 91.23%. Comparing base (a3c0f4f) to head (da4575c).

Additional details and impacted files

@@           Coverage Diff           @@
##           master     #858   +/-   ##
=======================================
  Coverage   91.23%   91.23%           
=======================================
  Files          11       11           
  Lines        1072     1072           
=======================================
  Hits          978      978           
  Misses         94       94           

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[KristofferC]

I don't think we should merge this until something like #859 is in. I didn't fully get the objection to that PR, it was a very long response but feel free to modify as you see fit.

[devmotion]

I benchmarked #859 on more examples and studied a few other alternatives, but none of them was able to consistently outperform this PR here.

Regarding #859 specifically, I think master/this PR is better since

• Conceptually, it implements exactly the underlying math (a zero partial means unperturbed, so its contribution is zero) for every number type, without relying on properties of the factor's type or implementation details of its arithmetic.
• Its lowering is predictable: The per-partial mask is branch-free and SIMD-friendly (fcmeq + bic on aarch64). There's no branch that LLVM could if-convert, speculate, or keep, and no additional code that could push operations over the inlining budget. Its cost is moderate and consistent: In none of the benchmarks above it's much slower than NaN-safe mode disabled (at most 1.59×), whereas the gains of this PR (ie Reduce the cost of NaN-safe mode by checking the scalar factor once #859) come with regressions of up to 3.8× that depend on N, the call site, and LLVM's heuristics.

[devmotion]

As I mentioned in #857, my strong opinion is that the correct mode should be the default, regardless of its performance impact. Unsafe performance optimization must be an opt-in, similar to @fastmath etc. One can always be fast by returning incorrect results - as ForwardDiff does by default. IMO this PR is a bugfix.

Any user can still opt into the faster but unsafe version by setting the preference.

[KristofferC]

regardless of its performance impact.

This is obviously false. You can test by putting a sleep(100) in the operation here and see what the effect will be on the issue tracker.

[KristofferC]

Conceptually,

Who cares about "conceptually", though? The point is to get the result and to get it quickly.

come with regressions of up to 3.8×

This is from benchmarking a completely useless a * x which says absolutely nothing about what the actual performance change will be in real scenarios.

[devmotion]

I disagree with both points.

There's one thing we want to achieve: zero partials out for zero partials in. The approach on master clearly enforces this, and it's easy to check and validate for any reader. The approach in #859 is indirect and does not work in general, hence ::Real fallbacks are needed, division fast path is limited to single-level Duals of IEEEFloats and multiplication fast path to possibly nested Duals of IEEEFloats, the correctness of which are crucially relying on the implementation detail of * and / for Duals. Hence my assessment that the master approach implements the mathematical property directly, without relying on properties of specific number types and implementation details.

Moreover, on my computer I just couldn't see a consistent advantage of #859 in the benchmarks. Neither in microbenchmarks nor in full derivative/gradient calls. Based on my experiments, I just can't reproduce that it is clearly superior.

On the other hand, what my experiments and study of some alternatives revealed was that the compiler much more consistently and reliably vectorizes the approach on master whereas the approach on #859 is compiled in a much less predictable way, even in the alternatives I tried. This seemed to explain also some of the performance regressions of #859.

All in all, in my experiments #859 doesn't seem to be a clear winner but come with a few caveats and disadvantages.