[devmotion]

Fixes #243, fixes #794.

Problem

The generic implementations of norm (LinearAlgebra.generic_norm2/generic_normp and StaticArrays' _norm/_norm_scaled) and hypot with four or more arguments (Base.Math._hypot) rescale by the maximum absolute value to avoid overflow and underflow. Their branches (iszero, isfinite, <) take partials into account for Duals (< since #791), so Duals end up in the wrong branch:

input	before	after
gradient(norm, [0.0, 0.0])	[NaN, NaN]	[0.0, 0.0]
gradient(x -> norm(x)^2, [0.0, 0.0])	[NaN, NaN]	[0.0, 0.0]
gradient(x -> norm(x, 3), [0.0, 0.0])	[NaN, NaN]	[0.0, 0.0]
norm([Dual(1e-200, 1e-100), Dual(0.0, 0.0)])	Dual(0.0, Inf) (wrong primal)	Dual(1e-200, 1e-100)
gradient(norm, [1e-200, 1e-200])	[Inf, Inf]	[0.707, 0.707]
norm(SA[Dual(0, 1), Dual(0, 1), Dual(0, 1)])	Dual(NaN, NaN)	Dual(0.0, 0.0)
norm((Dual(0.0, 1.0), Dual(0.0, 1.0)))	Dual(NaN, NaN)	Dual(0.0, 0.0)
norm([Dual(0.0, 1.0), Dual(1.0, 0.0)], 0)	Dual(2.0, 0.0) (wrong primal)	Dual(1.0, 0.0)
derivative(p -> norm([3.0, 4.0], p), 3.0)	MethodError	-0.304
hypot(Dual(0.0, 1.0), Dual(0.0, 1.0))	Dual(0.0, NaN)	Dual(0.0, 0.0)
hypot(Dual(0.0, 1.0), Dual(0.0, 1.0), Dual(0.0, 1.0), Dual(0.0, 1.0))	Dual(NaN, NaN)	Dual(0.0, 0.0)

All of these fail with NaN-safe mode enabled as well.

Approach

Similar to the rules in ChainRules, the norm of the primal values is computed with the existing (scale-robust) implementations, and the chain rule is applied with weights w_i = |u_i|^(p-1) * sign(u_i), where u_i = x_i / norm(x, p):

• The weights are bounded by 1 in absolute value for p ≥ 1 and p < 0, so they cannot overflow or underflow there. For 0 < p < 1 they are at least 1, as is the true derivative.
• All p use the same code path: the derivative is zero for p = 0, and only the first element with the largest (smallest) absolute value contributes for p = Inf (p = -Inf).
• At norm(x, p) == 0 and for zero elements x_i == 0, the subgradient zero is used. ChainRules uses the same convention (via pinv and sign).
• The derivative factors are computed in the precision of the primal norm (e.g., Float32 inputs give Float32 results).
• With respect to p: norm(x, p) / p * sum(|u_i|^p * log(|u_i|)). This is bounded and avoids the cancellation in sum(|x_i|^p * log(|x_i|)) - norm(x, p) * log(norm(x, p)).

The following are supported:

• complex arrays
• nested differentiation, including Hessians and all tag orderings of x and p
• StaticArrays, Adjoint/Transpose and sparse arrays (new weak dependency/extension)
• tuples, promoted to a common element type

Method ambiguities with LinearAlgebra, StaticArrays and SparseArrays are resolved explicitly.

hypot uses the same convention at zero. With four or more arguments it reuses the rule for norm: Base promotes the arguments and calls Base.Math._hypot with a tuple, and that internal function has the same signature on Julia 1.10 through 1.14-alpha. Additionally, the three-argument hypot no longer drops the partials of outer tags in nested differentiation (it used the tag-agnostic value/partials; this overlaps with #835).

Performance

Dense arrays of Duals with isbits values (including nested Duals) are reinterpreted without copying, so the primal norm of real Duals uses BLAS nrm2 on a strided view. The partials are summed without allocating (a lazy broadcast sum for arrays, sum(map(...)) for tuples and StaticArrays).

	master	this PR
norm(x), x::Vector{<:Dual}, length 100	268 ns	228 ns, no allocations
norm(x, Inf), length 100	185 ns	74 ns, no allocations
norm(x), length 3	8.5 ns	20 ns, no allocations

For very short vectors, master's single fused (but unscaled) pass is faster than the two passes needed here: the scale-robust primal norm, then the bounded weights.

Behavior changes

• norm(x, 1) uses the subgradient zero for zero elements, as for all other p. For example, gradient(x -> norm(x, 1), [1.0, 0.0, -2.0]) returns [1.0, 0.0, -1.0] instead of [1.0, 1.0, -1.0].
• With NaN-safe mode disabled, infinite elements yield NaN partials, e.g., norm([Dual(Inf, 0.0), Dual(1.0, 1.0)]) == Dual(Inf, NaN) instead of Dual(Inf, 0.0): their weight is Inf / Inf. With NaN-safe mode enabled, the result is Dual(Inf, 0.0).

Limitations

• Second-order derivatives through norm at zero remain incorrect (e.g., hessian(x -> norm(x)^2, zeros(2)) returns zeros instead of 2I): Since norm is not differentiable at zero, no choice of partials of norm(x) can recover the correct result. This is documented in the limitations section of the docs; sum(abs2, x) should be used instead.
• Tuples that mix Duals with non-Dual numbers (e.g., norm((Dual(0.0, 1.0), 0.0))) still use the generic implementation: dispatching on "some element is a Dual" would require type piracy.
• Arrays of arrays only benefit through the norm of their elements.

🤖 Generated with Claude Code

[codecov]

Codecov Report

❌ Patch coverage is 96.90722% with 3 lines in your changes missing coverage. Please review.
✅ Project coverage is 92.24%. Comparing base (70ef06c) to head (a390d78).

Files with missing lines	Patch %	Lines
src/dual.jl	96.55%	3 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #872      +/-   ##
==========================================
+ Coverage   91.24%   92.24%   +1.00%     
==========================================
  Files          11       12       +1     
  Lines        1085     1174      +89     
==========================================
+ Hits          990     1083      +93     
+ Misses         95       91       -4     

☔ 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.