Make x^y with Dual arguments with zero partials consistent with Real arguments

[devmotion]

Mathematically, a dual number with zero partials is just a regular number: Dual(v, 0) and v represent the same quantity. Hence derivative rules should return the same result regardless of whether an unperturbed argument is a Real or a Dual with zero partials.

Currently, this is not the case for ^ and NaNMath.pow if both arguments are Duals with the same tag:

julia> using ForwardDiff: Dual

julia> Dual(0.0, 0.0)^Dual(0.5, 1.0)  # same as `0.0^Dual(0.5, 1.0) == Dual(0.0, 0.0)`
Dual{Nothing}(0.0,NaN)

julia> Dual(0.0, 1.0)^Dual(0.0, 0.0)  # same as `Dual(0.0, 1.0)^0.0 == Dual(1.0, 0.0)`
Dual{Nothing}(1.0,NaN)

The second example returns NaN also in NaN-safe mode: there the NaN is not caused by a zero partial but by powval = 0 * 0^(-1).

This PR fixes both cases by setting powval to zero if iszero(vy) || isconstant(x), i.e., the same check as in the existing Dual^Real branch, which makes the Dual^Dual branch consistent with both the Dual^Real and the Real^Dual branch. For consistency, the check in the Dual^Real branch is now also written as iszero(y) || isconstant(x) instead of y == zero(y) || iszero(partials(x)) (which is equivalent). The new tests check that replacing an argument with a zero-partials Dual does not change the result.

With this PR, the derivatives at the special points agree with the Real case, SymPy, and WolframAlpha (default / NaN-safe mode where they differ):

	SymPy	WolframAlpha	master, Real	master, zero-partials Dual	this PR, zero-partials Dual
∂x^y/∂x at (0, 0)	0	0	0	NaN / NaN	0
∂x^y/∂x at (0, 0.5)	∞	∞	Inf	Inf	Inf
∂x^y/∂x at (0, -0.5)	-∞	-∞	-Inf	-Inf	-Inf
∂x^y/∂x at (2, 0)	0	0	0	0	0
∂x^y/∂y at (0, 0)	-∞	-∞	-Inf	NaN / -Inf	-Inf
∂x^y/∂y at (0, 0.5)	0	0	0	NaN / 0	0
∂x^y/∂y at (0, 1.5)	0	0	0	0	0
∂x^y/∂y at (0, -0.5)	-∞	-∞	-Inf	NaN / -Inf	-Inf

For x = 0, the SymPy and WolframAlpha values are the limits x -> 0+ of y * x^(y - 1) and x^y * log(x) (at (0, 0), ∂x^y/∂x is d/dx x^0 = 0). At (0, 0), ∂x^y/∂y = -∞ is also the limit of the one-sided difference quotients of y -> 0^y.

Note that this PR only makes the rules consistent for Duals whose partials are all zero. A Dual with only some zero partials still propagates NaNs in these directions in the default mode, so results depend on the chunk size, also with the existing check in the Dual^Real branch:

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

julia> ForwardDiff.gradient(f, [0.0, 1.0], ForwardDiff.GradientConfig(f, [0.0, 1.0], ForwardDiff.Chunk{1}()))
2-element Vector{Float64}:
 Inf
   1.0

julia> ForwardDiff.gradient(f, [0.0, 1.0], ForwardDiff.GradientConfig(f, [0.0, 1.0], ForwardDiff.Chunk{2}()))
2-element Vector{Float64}:
 Inf
 NaN

That's not specific to ^ (e.g., sqrt(v[1]) + v[2] returns [Inf, NaN] for any chunk size), and only NaN-safe mode handles it consistently. See #857.

🤖 Generated with Claude Code

[codecov]

Codecov Report

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

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #856      +/-   ##
==========================================
+ Coverage   91.23%   91.24%   +0.01%     
==========================================
  Files          11       11              
  Lines        1072     1074       +2     
==========================================
+ Hits          978      980       +2     
  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.