Fix zero outer-reverse second derivatives of saturated tanh by hoisting the tangent out of its JVP rule

[kodlan]

Fixes #40315.

For f(v) = tanh(u) * u * u with u = -1e8 * v, the second derivative at v = 1 is -2e16 and jacfwd(jacfwd) computes it exactly, but grad(grad) and any other combination with reverse mode outermost returned -0.0. tanh(u) rounds to exactly -1 there, and the correct value flows through terms that survive that rounding, so this was a real forward/reverse asymmetry, not expected precision loss. Notably clip and erf in the same program shape were fine; only tanh misbehaved.

The cause is the exact form of the default tanh JVP rule, the Herbie rewrite from 2020: (g + gans) * (1 - ans). It distributes the tangent g inside the (1 + ans) factor. Evaluated forward, g(1 + ans) collapses to an exact zero early (ans is exactly -1) and everything stays small. The transpose, which is what an outer grad evaluates, multiplies the cotangent by the large residuals first (intermediates around 1e24 for the repro), so the true answer term at 2e8 scale is absorbed by rounding before the matching pair cancels to zero. The transpose is mathematically exact; the evaluation order it induces is catastrophic.

The fix hoists g out: g * ((1 - ans) * (1 + ans)). The tangent path becomes a single multiplication by a factor computed purely from the primal, so the transpose sees the exact zero immediately, like the logistic rule next to it already does. Same operation count, and the compiled HLO for a tanh-MLP gradient is identical op for op (XLA fuses both forms the same way). The AccuracyMode.HIGHEST variant is untouched.

On first-derivative numerics, being fully upfront for anyone with golden gradient values: forward-mode JVPs with unit tangents are bitwise identical to before (checked over 200k points in f32 and f64), but reverse-mode grad(tanh) transposes to a differently associated expression and shifts in the last bit for many inputs. The shifts stay inside each form's error band, and the worst case strictly improves: the old g + g*ans add cancels at g's magnitude, with unbounded relative error deep in the negative saturation tail (observed errors up to about 1), while the new form stays within about 2 ulp of the value implied by the rounded primal everywhere we sampled, since 1 + ans is exact by the Sterbenz lemma for ans in [-1, -0.5].

The regression test checks grad(grad), grad(jacfwd) and jacfwd(jacfwd) agree on the exact second derivative at a saturated point, eager and jit, at a float32-safe scale; it fails before this change (-0.0) and passes after.