[devmotion]

Fixes #857.

This PR enables NaN-safe mode by default. With this PR:

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.

unsafe_mode preference

Since NaN-safe mode is now the default, the preference is renamed to make clear that turning it off trades correctness for speed: NaN-safe mode is disabled by setting the unsafe_mode preference to true (instead of setting nansafe_mode to false). The internal constant NANSAFE_MODE_ENABLED is renamed to UNSAFE_MODE_ENABLED accordingly, and the CI matrix sets unsafe_mode, so both modes remain tested.

The nansafe_mode preference is deprecated. If unsafe_mode is not set, nansafe_mode = false still disables NaN-safe mode. If nansafe_mode is set at all, loading ForwardDiff emits a deprecation warning with force = true, i.e. even with --depwarn=no:

┌ Warning: The `nansafe_mode` preference of ForwardDiff is deprecated. NaN-safe mode is enabled by default since ForwardDiff 1.5.0, and the preference to disable it is now called `unsafe_mode`. For now, `nansafe_mode = false` still disables NaN-safe mode unless `unsafe_mode` is set. Please remove the `nansafe_mode` preference with `Preferences.delete_preferences!(ForwardDiff, "nansafe_mode"; force = true)` and restart Julia.
│   caller = run_module_init(mod::Module, i::Int64) at loading.jl:1470
└ @ Base ./loading.jl:1470

Note that with --depwarn=error this warning becomes an error, so loading ForwardDiff fails until the old preference is removed.

The documentation section is rewritten: disabling NaN-safe mode is described as trading correctness for performance, similar to @fastmath, and the outdated "~5%-10%" claim is removed since the cost depends on the function and the architecture (see the numbers below). The deprecated preference is mentioned in a compat admonition.

Reducing the cost of NaN-safe mode

This PR includes #859 (Kristoffer's commit is kept as is), which checks the scalar factor once instead of masking every partial, plus follow-ups that make its fast path stricter:

• Following ChainRulesCore's strong_mul (Avoid NaN-propagation in scalar rules ChainRulesCore.jl#551), the masking helpers are called _strong_mul/_strong_div. _mul_is_strong(V, x)/_div_is_strong(V, x) decide whether plain *// already treat zero partials of type V as strong zeros. With NaN-safe mode disabled they are constant true, so nothing changes there.
• Plain *// is only used if both the partials' element type and the factor are known Base number types (IEEEFloat, BitInteger, Bool, BigInt, BigFloat, Rationals of these, AbstractIrrational) or dual numbers of those. Any other Real type always takes the masked path.
• For *, the factor must be finite (in all components).
• For /, the divisor must be non-zero and not NaN; for dual numbers, its value must satisfy this (recursively) and its partials must be finite. ±Inf does not force the masked path.
• _mul_partials keeps Reduce the cost of NaN-safe mode by checking the scalar factor once #859's single check for both factors: choosing plain or masked separately for each factor made the first hessian call compile 2–3× longer, for at most 4–6% runtime.

Division of dual numbers

_div_partials (used by Dual / Dual) multiplied by inv(bval) and divided by bval * bval, which overflow or underflow for tiny and huge bval. For example, on master the derivative of x -> 2x / x is -Inf at 1e-200 and 2.0e-200 at 1e200, and dividing a dual zero by a dual number with value 1e-200 yields NaN partials (which would also break the fast path for dual divisors above). Following the quotient rule, _div_partials now computes (a - (aval / bval) * b) / bval, avoiding inv as Dual / Real and Real / Dual already do. Dividing each partial is slower than multiplying by inv(bval) in division-heavy code (see below), but correct.

Timings

Minimum times (Chairmarks @be, seconds=3), each variant in a fresh process, two alternating runs. "NaN-safe disabled" is this PR with unsafe_mode = true; "master" and "#859" are with nansafe_mode = true:

Benchmark	NaN-safe disabled	master	#859	this PR
gradient(rosenbrock, rand(1000))	506–520 μs	737–758 μs	606 μs	591–607 μs
hessian(rosenbrock, rand(10))	3.28–3.33 μs	5.10–5.20 μs	3.04–3.05 μs	2.99 μs
hessian(rosenbrock, rand(20))	17.6–18.0 μs	33.9–34.0 μs	16.3–16.4 μs	15.3–15.5 μs
hessian(logreg_nll, randn(8)), 500 × 8	131–134 μs	205–211 μs	173 μs	160–164 μs
gradient(x -> sum(exp, x), rand(100))	7.60–7.81 μs	7.78–7.93 μs	6.74–6.85 μs	6.69–6.82 μs
jacobian(cumprod, rand(10))*	611–617 ns	609–611 ns	588–590 ns	596–597 ns
derivative(x -> sin(x) * exp(x) / x, 1.3)	16.3–16.7 ns	16.1–16.5 ns	17.3–17.7 ns	18.6–19.1 ns

* Measured on its own in a fresh process. Within the full script below, this PR consistently takes 733–761 ns for it (#859: 565–587 ns), i.e. it depends on the benchmarks that ran before.

Division change, NaN-safe mode enabled, with ratios(x) = sum(x[i] / x[i+1] for i in 1:length(x)-1) and rosen_div(x) = sum((1 - x[i])^2 / (1 + x[i+1]^2) for i in 1:length(x)-1) on rand(n) .+ 0.5:

Benchmark	inv(bval), aval / (bval * bval)	(a - (aval / bval) * b) / bval
gradient(rosen_div, 1000)	490 μs	585–589 μs
hessian(rosen_div, 20)	15.1–15.3 μs	17.5–17.8 μs
hessian(ratios, 10)	2.61–2.65 μs	2.76–2.79 μs
gradient(ratios, 1000)	646–663 μs	653–709 μs
derivative(x -> x / (1 + x^2), 1.3)	2.3 ns	2.3 ns

Benchmark script

using ForwardDiff, Chairmarks, LinearAlgebra, 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

# Negative log-likelihood of a logistic regression model
function logreg_nll(β, X, y)
    result = zero(eltype(β))
    for i in axes(X, 1)
        η = dot(view(X, i, :), β)
        result += log1p(exp(η)) - y[i] * η
    end
    return result
end

Random.seed!(1)
x1000 = rand(1000); x100 = rand(100); x10 = rand(10); x20 = rand(20)
const X = randn(500, 8); const y = rand(Bool, 500); β = randn(8) ./ 4
nll(β) = logreg_nll(β, X, y)

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

The division benchmark uses the same setup with ratios and rosen_div as defined above, on rand(n) .+ 0.5.

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)

Tests

test/PreferencesTest.jl checks the deprecated preference in subprocesses (resulting mode and warning text). test/PartialsTest.jl checks strong zeros for unknown number types, non-finite dual factors and divisors, and tiny nested divisors; test/DerivativeTest.jl checks division by tiny and huge numbers. Pkg.test() passes locally both with the new default (9784 tests) and with unsafe_mode = true (9664 tests).

🤖 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.42%. Comparing base (a3c0f4f) to head (33f771b).

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #858      +/-   ##
==========================================
+ Coverage   91.23%   91.42%   +0.19%     
==========================================
  Files          11       11              
  Lines        1072     1096      +24     
==========================================
+ Hits          978     1002      +24     
  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.

[devmotion]

I put together a larger set of hopefully somewhat realistic comparisons of NaN-safe mode on master (which #858 enables by default) and in #859. I deliberately included cases where I expected #859 to be slower.

All workloads are hand-written without dependencies besides ForwardDiff and use the public API with preallocated configs and default chunk sizes: Rosenbrock; logistic regression, Weibull survival regression, sigmoid Emax PK-PD model, MLP, polynomial and Michaelis-Menten least squares; ODE right-hand side Jacobians (ROBER, 2D Brusselator, 1D Fisher-KPP); gradients through RK4 solves of PK models; Newton's method for Kepler's equation; and the benchmarks from #858 and #859.

Results

Times are relative to master with NaN-safe mode disabled; in the last column, < 1 means #859 is faster than master with NaN-safe mode enabled (with bootstrap 95% CI).

Workload	NaN-safe off, master	NaN-safe on, master	NaN-safe on, #859	#859 / master (NaN-safe on)
Hessians
Rosenbrock, n = 10	2.33 μs	1.91×	0.92×	0.48× (0.46–0.50)
Rosenbrock, n = 20	17.40 μs	2.03×	0.91×	0.45× (0.42–0.46)
logistic regression NLL, 500 × 8	115.3 μs	1.70×	1.73×	1.02× (0.97–1.04)
Weibull regression NLL, 300 × 7	95.90 μs	3.07×	1.34×	0.44× (0.42–0.44)
sigmoid Emax PK-PD NLL, 200 obs	60.27 μs	2.04×	1.30×	0.64× (0.60–0.65)
Gradients, Jacobians, derivatives
gradient: Rosenbrock, n = 1000	552.5 μs	1.45×	1.15×	0.79× (0.75–0.81)
gradient: logistic regression NLL (X * β), 500 × 8	13.42 μs	1.14×	1.04×	0.91× (0.87–0.94)
gradient: logistic regression NLL (loops), 500 × 8	11.62 μs	1.27×	1.51×	1.19× (1.16–1.21)
gradient: Weibull regression NLL, 300 × 7	22.69 μs	1.03×	1.04×	1.01× (0.98–1.02)
gradient: MLP 4-8-1 loss, 49 params, batch 64	43.46 μs	1.18×	1.00×	0.85× (0.79–0.91)
gradient: sigmoid Emax PK-PD NLL, 200 obs	23.52 μs	1.05×	1.03×	0.98× (0.94–1.01)
gradient: quadratic least squares (broadcast), 1000 points	5.57 μs	1.61×	1.53×	0.95× (0.93–0.97)
gradient: Michaelis-Menten least squares (broadcast), 1000 points	5.46 μs	1.62×	2.14×	1.32× (1.28–1.39)
gradient: RK4 PK 1-compartment, 3 params	31.27 μs	1.14×	1.07×	0.93× (0.88–1.01)
gradient: RK4 PK 2-compartment, 5 params	50.38 μs	1.52×	1.31×	0.87× (0.82–0.91)
gradient: sum(exp, x), n = 100	7.62 μs	1.04×	0.86×	0.83× (0.81–0.87)
gradient: sum(xi -> xi^x[1], x), n = 100	40.56 μs	1.05×	1.35×	1.28× (1.24–1.34)
Jacobian: ROBER RHS, 3 states	24.6 ns	1.14×	1.01×	0.88× (0.84–0.92)
Jacobian: 1D Fisher-KPP RHS, 16 states	539.0 ns	1.08×	1.07×	0.99× (0.96–1.02)
Jacobian: 2D Brusselator RHS, 128 states	45.46 μs	1.09×	1.01×	0.92× (0.88–1.02)
Jacobian: cumprod, n = 10	313.4 ns	1.04×	1.01×	0.97× (0.94–1.01)
derivative: Newton for Kepler's equation, 1000 solves	144.4 μs	1.02×	0.98×	0.97× (0.95–0.98)
derivative: sin(x) * exp(x) / x	16.4 ns	1.06×	1.13×	1.06× (1.04–1.11)

With NaN-safe mode disabled, #859 is within 0.98×–1.02× of master everywhere.

Summary

• Hessians: NaN-safe mode on master costs 1.7×–2.0×, and about 3× for the Weibull Hessian (3.0×–3.2× in repeated runs with other data). Reduce the cost of NaN-safe mode by checking the scalar factor once #859 cuts this to 0.44×–0.64× of master, except for logistic regression (1.02×). The overhead on master comes mostly from masking nested duals, not from branch mispredictions.
• First-order code: Reduce the cost of NaN-safe mode by checking the scalar factor once #859 is mostly faster, but slower for logistic regression written with loops (1.19×), Michaelis-Menten least squares (1.32×), the xi^x[1] gradient (1.28×) and the scalar derivative (1.06×). These are not caused by the inlining failure I described in Reduce the cost of NaN-safe mode by checking the scalar factor once #859. The same logistic regression model is 9% faster with Reduce the cost of NaN-safe mode by checking the scalar factor once #859 when written with X * β and broadcasting, and 19% slower when written with loops.
• Exploit symmetry in the Hessian #844 doesn't change this: in a separate run of the Hessians, both the overhead of NaN-safe mode on master and Reduce the cost of NaN-safe mode by checking the scalar factor once #859's advantage are the same with Exploit symmetry in the Hessian #844.

Overall, #859 is not uniformly faster. It is clearly worthwhile for Hessians, but for first-order code the gains are modest, some realistic workloads regress by up to 1.3×, and whether #859 is faster or slower can flip just based on how the same model is written (X * β vs. loops). Since the gains are concentrated in nested duals, one option might be a hybrid that keeps the per-partial mask for first-order partials and uses the factor check only for nested duals; I haven't tested whether that keeps the Hessian gains without the first-order regressions. Independent of performance, I still think the factor check should be restricted to IEEEFloats (and Duals of them), for the reasons in #859.

Setup and caveats

• Apple M2 Pro, Julia 1.13.1, Chairmarks; aarch64 only, x86 is untested.
• Every configuration ran in a fresh process, 10 interleaved rounds in random order. Per process, each workload was sampled with @be for 1 s (0.5 s for the 7 workloads added in a second run: loop logistic regression, Weibull, Fisher-KPP, quadratic and Michaelis-Menten least squares, Kepler); the table shows the median over rounds of the per-process medians. The Exploit symmetry in the Hessian #844 comparison used 8 rounds at 0.5 s.
• Results agree between master and Reduce the cost of NaN-safe mode by checking the scalar factor once #859 up to signed zeros, except for a last-bit difference in Michaelis-Menten: the per-element terms are identical, but the final sum is vectorized differently (as also happens between NaN-safe mode on and off on master).

Benchmark script (suite.jl)

Run once per configuration in a fresh process, with nansafe_mode set in LocalPreferences.toml: julia --project=<env> suite.jl <outfile> 0.5. It writes, per workload, the median and minimum time and a hash of the result.

# Realistic ForwardDiff workloads for comparing NaN-safe mode implementations.
# Usage: julia --project=<env> suite.jl <outfile> [seconds per benchmark]
# Writes one tab-separated line per benchmark: name, median ns, minimum ns, result hash.

using ForwardDiff, DiffResults, LinearAlgebra, Random, Chairmarks

const SECONDS = length(ARGS) >= 2 ? parse(Float64, ARGS[2]) : 0.5

# Hash of a result, normalizing signed zeros (NaN-safe mode implementations may differ in
# the sign of zero partials)
_normalize(x::Real) = iszero(x) ? zero(x) : x
result_hash(x::Real) = hash(_normalize(x))
result_hash(x::AbstractArray) = hash(map(_normalize, x))
result_hash(x::DiffResults.DiffResult) = hash(map(result_hash, (DiffResults.value(x), DiffResults.derivative(x))))
result_hash(x::Tuple) = hash(map(result_hash, x))

const BENCHMARKS = Pair{String,Any}[]
add!(name, f) = push!(BENCHMARKS, name => f)

# Extended Rosenbrock (Optim/CUTEst test problem), as in #858
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
let x = rand(MersenneTwister(1), 1000), out = similar(x), cfg = ForwardDiff.GradientConfig(rosenbrock, x)
    add!("gradient: Rosenbrock, n = 1000", () -> ForwardDiff.gradient!(out, rosenbrock, x, cfg))
end
for n in (10, 20)
    let x = rand(MersenneTwister(1), n), out = zeros(n, n), cfg = ForwardDiff.HessianConfig(rosenbrock, x)
        add!("Hessian: Rosenbrock, n = $n", () -> ForwardDiff.hessian!(out, rosenbrock, x, cfg))
    end
end

# Logistic regression negative log-likelihood (GLM fitting, Newton/IRLS, Laplace)
function logreg_nll(β, X, y)
    η = X * β
    return sum(log1p.(exp.(η)) .- y .* η)
end
let rng = MersenneTwister(2), X = randn(rng, 500, 8), y = Float64.(rand(rng, Bool, 500)),
    β = 0.1 * randn(rng, 8), f = β -> logreg_nll(β, X, y),
    gout = similar(β), gcfg = ForwardDiff.GradientConfig(f, β),
    hout = zeros(8, 8), hcfg = ForwardDiff.HessianConfig(f, β)
    add!("gradient: logistic regression NLL (`X * β`), 500 × 8", () -> ForwardDiff.gradient!(gout, f, β, gcfg))
    add!("Hessian: logistic regression NLL, 500 × 8", () -> ForwardDiff.hessian!(hout, f, β, hcfg))
end

# Small MLP (4 → 8 → 1, tanh), mean squared error over a batch of 64 samples,
# gradient w.r.t. all 49 parameters
function mlp_loss(θ, X, y)
    W1 = reshape(view(θ, 1:32), 8, 4)
    b1 = view(θ, 33:40)
    w2 = view(θ, 41:48)
    b2 = θ[49]
    h = tanh.(W1 * X .+ b1)
    ŷ = vec(w2' * h) .+ b2
    return sum(abs2, ŷ .- y) / length(y)
end
let rng = MersenneTwister(3), X = randn(rng, 4, 64), y = randn(rng, 64), θ = 0.3 * randn(rng, 49),
    f = θ -> mlp_loss(θ, X, y), out = similar(θ), cfg = ForwardDiff.GradientConfig(f, θ)
    add!("gradient: MLP 4-8-1 loss, 49 params, batch 64", () -> ForwardDiff.gradient!(out, f, θ, cfg))
end

# Sigmoid Emax (Hill) PK-PD model with power error model, negative log-likelihood
# θ = (E0, Emax, EC50, γ, σ); `C^γ` and `EC50^γ` (`Real^Dual`, `Dual^Dual`), `f^0.8` (`Dual^Real`)
function emax_nll(θ, C, y)
    E0, Emax, EC50, γ, σ = θ
    nll = zero(eltype(θ))
    for i in eachindex(C, y)
        Cγ = C[i]^γ
        f = E0 + Emax * Cγ / (EC50^γ + Cγ)
        sd = σ * f^0.8
        nll += log(sd) + ((y[i] - f) / sd)^2 / 2
    end
    return nll
end
let rng = MersenneTwister(4), C = 10 .^ range(-1, 2; length = 200), θ = [1.0, 10.0, 5.0, 1.5, 0.1],
    y = [θ[1] + θ[2] * c^θ[4] / (θ[3]^θ[4] + c^θ[4]) for c in C] .* (1 .+ 0.1 .* randn(rng, 200)),
    f = θ -> emax_nll(θ, C, y), gout = similar(θ), gcfg = ForwardDiff.GradientConfig(f, θ),
    hout = zeros(5, 5), hcfg = ForwardDiff.HessianConfig(f, θ)
    add!("gradient: sigmoid Emax PK-PD NLL, 200 obs", () -> ForwardDiff.gradient!(gout, f, θ, gcfg))
    add!("Hessian: sigmoid Emax PK-PD NLL, 200 obs", () -> ForwardDiff.hessian!(hout, f, θ, hcfg))
end

# ODE right-hand side Jacobians, as computed by stiff solvers
# ROBER (3 states)
function rober!(du, u)
    k1, k2, k3 = 0.04, 3e7, 1e4
    y1, y2, y3 = u
    du[1] = -k1 * y1 + k3 * y2 * y3
    du[2] = k1 * y1 - k2 * y2^2 - k3 * y2 * y3
    du[3] = k2 * y2^2
    return du
end
let u = [1.0, 1e-4, 0.5], du = similar(u), J = zeros(3, 3), cfg = ForwardDiff.JacobianConfig(rober!, du, u)
    add!("Jacobian: ROBER RHS, 3 states", () -> ForwardDiff.jacobian!(J, rober!, du, u, cfg))
end
# 2D Brusselator (SciML benchmark), 8 × 8 grid, 128 states
const BRUSS_N = 8
function brusselator!(du, u)
    N = BRUSS_N
    A, B, α = 3.4, 1.0, 10.0 * (N - 1)^2
    xs = range(0, 1; length = N)
    U = reshape(view(u, 1:N^2), N, N)
    V = reshape(view(u, N^2+1:2N^2), N, N)
    dU = reshape(view(du, 1:N^2), N, N)
    dV = reshape(view(du, N^2+1:2N^2), N, N)
    lim(k) = k == 0 ? N : (k == N + 1 ? 1 : k)
    @inbounds for j in 1:N, i in 1:N
        ip, im, jp, jm = lim(i + 1), lim(i - 1), lim(j + 1), lim(j - 1)
        u2v = U[i, j]^2 * V[i, j]
        forcing = ((xs[i] - 0.3)^2 + (xs[j] - 0.6)^2 <= 0.01) ? 5.0 : 0.0
        dU[i, j] = α * (U[im, j] + U[ip, j] + U[i, jp] + U[i, jm] - 4U[i, j]) + B + u2v - (A + 1) * U[i, j] + forcing
        dV[i, j] = α * (V[im, j] + V[ip, j] + V[i, jp] + V[i, jm] - 4V[i, j]) + A * U[i, j] - u2v
    end
    return du
end
let rng = MersenneTwister(5), u = rand(rng, 2BRUSS_N^2), du = similar(u), J = zeros(length(u), length(u)),
    cfg = ForwardDiff.JacobianConfig(brusselator!, du, u)
    add!("Jacobian: 2D Brusselator RHS, 128 states", () -> ForwardDiff.jacobian!(J, brusselator!, du, u, cfg))
end

# Forward sensitivities through a fixed-step RK4 solve of PK models,
# sum of squared residuals at 12 observation times
@inline function rk4_step(f, u, p, dt)
    k1 = f(u, p)
    k2 = f(u .+ (dt / 2) .* k1, p)
    k3 = f(u .+ (dt / 2) .* k2, p)
    k4 = f(u .+ dt .* k3, p)
    return u .+ (dt / 6) .* (k1 .+ 2 .* k2 .+ 2 .* k3 .+ k4)
end
# One-compartment model with first-order absorption, p = (ka, CL, V)
pk1cmt(u, p) = (-p[1] * u[1], p[1] * u[1] - p[2] / p[3] * u[2])
# Two-compartment model with first-order absorption, p = (ka, CL, V1, Q, V2)
pk2cmt(u, p) = (-p[1] * u[1],
                p[1] * u[1] - (p[2] + p[4]) / p[3] * u[2] + p[4] / p[5] * u[3],
                p[4] / p[3] * u[2] - p[4] / p[5] * u[3])
function pk_sse(f, p, u0, obs_steps, obs, dt, cmt)
    u = map(x -> x * one(eltype(p)), u0)
    sse = zero(eltype(p))
    step = 0
    for (k, n) in enumerate(obs_steps)
        while step < n
            u = rk4_step(f, u, p, dt)
            step += 1
        end
        sse += (obs[k] - u[cmt] / p[3])^2
    end
    return sse
end
let obs_steps = [5, 10, 20, 40, 60, 80, 120, 160, 240, 320, 400, 480], dt = 0.05
    for (name, f, p, u0) in (("1-compartment, 3 params", pk1cmt, [1.2, 3.0, 20.0], (100.0, 0.0)),
                             ("2-compartment, 5 params", pk2cmt, [1.2, 3.0, 20.0, 2.0, 40.0], (100.0, 0.0, 0.0)))
        obs = [5 * exp(-0.01 * n) for n in obs_steps]
        g = p -> pk_sse(f, p, u0, obs_steps, obs, dt, 2)
        out = similar(p)
        cfg = ForwardDiff.GradientConfig(g, p)
        add!("gradient: RK4 PK $name, 480 steps", () -> ForwardDiff.gradient!(out, g, p, cfg))
    end
end

# The benchmarks from #858 and from the review in #859
let x = rand(MersenneTwister(1), 100), out = similar(x), f = x -> sum(exp, x), cfg = ForwardDiff.GradientConfig(f, x)
    add!("gradient: sum(exp, x), n = 100", () -> ForwardDiff.gradient!(out, f, x, cfg))
end
let x = rand(MersenneTwister(1), 10), J = zeros(10, 10), cfg = ForwardDiff.JacobianConfig(cumprod, x)
    add!("Jacobian: cumprod, n = 10", () -> ForwardDiff.jacobian!(J, cumprod, x, cfg))
end
let f = x -> sin(x) * exp(x) / x, x = Ref(1.3)
    add!("derivative: sin(x) * exp(x) / x", () -> ForwardDiff.derivative(f, x[]))
end
let x = rand(MersenneTwister(1), 100), out = similar(x), f = x -> sum(xi -> xi^x[1], x), cfg = ForwardDiff.GradientConfig(f, x)
    add!("gradient: sum(xi -> xi^x[1], x), n = 100", () -> ForwardDiff.gradient!(out, f, x, cfg))
end

# Workloads exercising the code paths where #859 regressed in microbenchmarks

# Logistic regression NLL written with explicit loops (`Float64 * Dual{8}` per term)
function logreg_nll_loop(β, X, y)
    nll = zero(eltype(β))
    for i in axes(X, 1)
        η = zero(eltype(β))
        for j in axes(X, 2)
            η += X[i, j] * β[j]
        end
        nll += log1p(exp(η)) - y[i] * η
    end
    return nll
end
let rng = MersenneTwister(2), X = randn(rng, 500, 8), y = Float64.(rand(rng, Bool, 500)),
    β = 0.1 * randn(rng, 8), f = β -> logreg_nll_loop(β, X, y), out = similar(β), cfg = ForwardDiff.GradientConfig(f, β)
    add!("gradient: logistic regression NLL (loops), 500 × 8", () -> ForwardDiff.gradient!(out, f, β, cfg))
end

# Weibull survival regression NLL, 7 covariates + shape (8 params);
# `(t/λ)^k` is `Dual^Dual`, the linear predictor is `Float64 * Dual{8}`
function weibull_nll(θ, X, t, δ)
    k = θ[end]
    nll = zero(eltype(θ))
    for i in eachindex(t, δ)
        η = zero(eltype(θ))
        for j in axes(X, 2)
            η += X[i, j] * θ[j]
        end
        z = t[i] / exp(η)
        nll -= δ[i] * (log(k) - η + (k - 1) * log(z)) - z^k
    end
    return nll
end
let rng = MersenneTwister(6), X = randn(rng, 300, 7), t = 0.1 .+ 2 .* rand(rng, 300), δ = Float64.(rand(rng, 300) .< 0.7),
    θ = [0.1 * randn(rng, 7); 1.5], f = θ -> weibull_nll(θ, X, t, δ),
    gout = similar(θ), gcfg = ForwardDiff.GradientConfig(f, θ), hout = zeros(8, 8), hcfg = ForwardDiff.HessianConfig(f, θ)
    add!("gradient: Weibull regression NLL, 300 × 7", () -> ForwardDiff.gradient!(gout, f, θ, gcfg))
    add!("Hessian: Weibull regression NLL, 300 × 7", () -> ForwardDiff.hessian!(hout, f, θ, hcfg))
end

# Jacobian of a 1D reaction-diffusion (Fisher-KPP) RHS, 16 states (chunk size 8)
function fisher_kpp!(du, u)
    n = length(u)
    D, r = 0.1 * (n + 1)^2, 1.0
    for i in 1:n
        left = i == 1 ? zero(u[i]) : u[i-1]
        right = i == n ? zero(u[i]) : u[i+1]
        du[i] = D * (left - 2 * u[i] + right) + r * u[i] * (1 - u[i])
    end
    return du
end
let u = rand(MersenneTwister(7), 16), du = similar(u), J = zeros(16, 16), cfg = ForwardDiff.JacobianConfig(fisher_kpp!, du, u)
    add!("Jacobian: 1D Fisher-KPP RHS, 16 states", () -> ForwardDiff.jacobian!(J, fisher_kpp!, du, u, cfg))
end

# Polynomial least squares with broadcasting, 3 params, 1000 points
poly_sse(θ, t, y) = sum(abs2, y .- (θ[1] .+ θ[2] .* t .+ θ[3] .* t .^ 2))
let rng = MersenneTwister(8), t = range(0, 1; length = 1000), y = 1 .+ 2 .* t .- t .^ 2 .+ 0.1 .* randn(rng, 1000),
    θ = [0.5, 1.0, 0.0], f = θ -> poly_sse(θ, t, y), out = similar(θ), cfg = ForwardDiff.GradientConfig(f, θ)
    add!("gradient: quadratic least squares (broadcast), 1000 points", () -> ForwardDiff.gradient!(out, f, θ, cfg))
end

# Michaelis-Menten least squares with broadcasting, 2 params, 1000 points (division-heavy)
mm_sse(θ, S, v) = sum(abs2, v .- θ[1] .* S ./ (θ[2] .+ S))
let rng = MersenneTwister(9), S = range(0.1, 10; length = 1000), v = 2 .* S ./ (1.5 .+ S) .+ 0.05 .* randn(rng, 1000),
    θ = [1.0, 1.0], f = θ -> mm_sse(θ, S, v), out = similar(θ), cfg = ForwardDiff.GradientConfig(f, θ)
    add!("gradient: Michaelis-Menten least squares (broadcast), 1000 points", () -> ForwardDiff.gradient!(out, f, θ, cfg))
end

# Newton's method for Kepler's equation `E - e sin(E) = M`, derivatives via ForwardDiff
function kepler(M, e)
    g(E) = E - e * sin(E) - M
    E = M
    for _ in 1:6
        E -= g(E) / ForwardDiff.derivative(g, E)
    end
    return E
end
let Ms = range(0, 2π; length = 1000)
    add!("derivative: Newton for Kepler's equation, 1000 solves", () -> sum(M -> kepler(M, 0.3), Ms))
end

#######
# Run #
#######

function run_suite(outfile)
    open(outfile, "w") do io
        for (name, f) in BENCHMARKS
            h = result_hash(f())
            b = @be f() seconds = SECONDS
            println(io, name, '\t', Chairmarks.median(b).time * 1e9, '\t', Chairmarks.minimum(b).time * 1e9, '\t', h)
            flush(io)
        end
    end
end

run_suite(ARGS[1])

Driver (run.jl): environments and interleaved rounds

# Runs suite.jl for every configuration in fresh, interleaved processes.
# Usage: julia run.jl <master checkout> <PR checkout> <rounds> [seconds per benchmark]

using Pkg, Random

const SOURCES = Dict("master" => ARGS[1], "pr" => ARGS[2])
const CONFIGS = [(src, nansafe) for src in ("master", "pr") for nansafe in (false, true)]
config_name((src, nansafe)) = "$(src)-nansafe=$(nansafe)"

function setup_env(cfg)
    src, nansafe = cfg
    dir = joinpath(@__DIR__, "envs", config_name(cfg))
    mkpath(dir)
    write(joinpath(dir, "LocalPreferences.toml"), "[ForwardDiff]\nnansafe_mode = $nansafe\n")
    Pkg.activate(dir) do
        Pkg.develop(path = SOURCES[src])
        Pkg.add(["Chairmarks", "DiffResults"])
        Pkg.precompile()
    end
    return dir
end

function main(rounds, seconds)
    envs = Dict(cfg => setup_env(cfg) for cfg in CONFIGS)
    resdir = joinpath(@__DIR__, "results")
    mkpath(resdir)
    rng = MersenneTwister(42)
    for round in 1:rounds, cfg in shuffle(rng, CONFIGS)
        out = joinpath(resdir, "$(config_name(cfg))-round=$(round).tsv")
        run(`$(Base.julia_cmd()) --startup-file=no --project=$(envs[cfg]) $(joinpath(@__DIR__, "suite.jl")) $out $seconds`)
    end
end

main(parse(Int, ARGS[3]), length(ARGS) >= 4 ? ARGS[4] : "0.5")

[devmotion]

That being said, I want to repeat once more what I already wrote repeatedly: In my opinion this PR is a long overdue bug fix and hence the performance implications are a side effect but not relevant for whether the bug should be fixed or not and this PR should be merged or not. I can't count how often I ran into this correctness issue, in e.g. Turing, Distributions and Pumas. At least I know that I could try turning on the correct mode, in contrast to many users who don't. We're only having this discussion since for whatever reason what should have been a simple bugfix resolving a correctness issue ended up as a secondary non-default opt-in mode.

[devmotion]

@KristofferC I assume the windows CI failure on 1.14.0-alpha1 is caused by an upstream issue, maybe Julia or LLVM

[devmotion]

Upstream issue: JuliaLang/julia#64013