[devmotion]

Ports the fix for JuliaDiff/DiffRules.jl#102 to ChainRules, and fixes related cancellations.

tanh

1 - Ω ^ 2 with Ω = tanh(x) cancels catastrophically for large |x|. From x ≈ 19.07, Ω rounds to 1 and the derivative is 0. Rearranging Ω, e.g. (1 - Ω) * (1 + Ω), cannot fix this, so the derivative is computed from x: sech(x) ^ 2 for |real(x)| ≥ 1. The cheap 1 - Ω ^ 2 is kept for |real(x)| < 1, where |sech(x)|^2 ≥ 0.42, so it doesn't suffer from cancellation.

Relative error vs. BigFloat:

x	before	after
4	2e-14	1e-16
14	1e-5	8e-17
19	0.77	2e-16
25	1.0	2e-16
30 + 1im	0.42	3e-16

The Hermitian/Symmetric matrix tanh rule calls the scalar frule, so it picks up the fix too.

1 - x ^ 2 / x ^ 2 - 1 near |x| = 1

asin, acos, atanh, acoth, asech, asind, acosd and the Real methods of asec, acsc, asecd, acscd now use (1 - x) * (1 + x) resp. (x - 1) * (x + 1). Near |x| = 1, the factor 1 - x (resp. 1 + x) is computed exactly, which avoids the cancellation in 1 - x ^ 2. Example: atanh at x = 1 - 1e-8 goes from 6e-10 to 1e-16 relative error. On the complex branch cuts (signed-zero imaginary parts), the signs match the previous rules.

The generic (complex) 1 - x ^ -2 forms of asec/acsc/asecd/acscd are not changed. The obvious accurate rewrite, (x - 1) * (x + 1) / x ^ 2, puts signed-zero inputs on (-1, 1) on the wrong side of the branch cut.

References: this is the standard rewrite of x² - y² as (x - y)(x + y), which replaces catastrophic cancellation with a benign one (D. Goldberg, What Every Computer Scientist Should Know About Floating-Point Arithmetic, 1991, section "Cancellation"). For 1/2 ≤ x ≤ 2, 1 - x is computed exactly (ibid., Theorem 11, the Sterbenz lemma). Julia Base uses the same factorization in atanh(::Complex), and sqrt(1-z)*sqrt(1+z) in asin(::Complex). ChainRules' existing acosh rule, inv(sqrt(x - 1) * sqrt(x + 1)), is similar.

Complex 1 + x ^ 2 near x = ±i

asinh, atan, acot, atand and acotd have new Complex methods that use (1 + ix) * (1 - ix), which is the same rewrite as above: near x = ±i, one factor is computed exactly. For example, at x = 1e-9 + (1 - 1e-8)im, the relative error drops from 3e-10 (asinh) and 6e-10 (the others) to below eps(). In asinh, ix is constructed as Complex(-imag(x), real(x)), because im * x loses the sign of a zero real part, which matters on the branch cuts. On a grid of signed zeros, ±1, ±2, ±1e300, ±Inf and NaN in both components, the new rules agree with the old ones up to rounding. Complex acsch (1 + x ^ -2) is not changed, for the same branch-cut reason as complex asec/acsc.

Overflow and underflow of x ^ 2

• asinh(x::Real) now returns inv(abs(x)) if x ^ 2 overflows. This is exact up to rounding, since 1 / |x| matches 1 / sqrt(1 + x ^ 2) to within eps / 2 already for |x| > 1 / sqrt(eps) = 2 ^ 26, long before x ^ 2 overflows at |x| ≈ 1.3e154. Before, asinh at 1e200 gave 0.0 instead of 1e-200.
• atan(y, x) now divides x and y by s = max(|x|, |y|) first: with xs = x / s, ys = y / s and v = xs ^ 2 + ys ^ 2 (which is in [1, 2]), the partials are (xs / v / s, -ys / v / s). Before, (1e-200, 1e-200) gave NaN instead of ±5e199, and (1e200, 3e200) gave 0.0; Int inputs could also overflow silently in x ^ 2 + y ^ 2. The maximum relative error over the full range is about 1.6 eps. For ordinary inputs, results can differ from before in the last bits. This also works for atan_fast.
• hypot would also avoid the overflow, but it is ~10× slower than the old formulas, since its branches and correction step prevent SIMD.

Possible regressions

• Values: none found, apart from atan(y, x) with infinite arguments (next item). Overflow/limit cases (±800, ±Inf, NaN, 800 + 1im, Inf + 1im, near the poles at iπ/2) give the same results as before. Element types (Float16, Float32, Int, ComplexF64) are preserved. Float32 improves, e.g. tanh at 30f0: 0 before, 3.5e-26 now.
• Values of atan(y, x) when exactly one argument is infinite: both partials are now NaN. Before, one of them was ±0.0 (the correct limit) and the other NaN, e.g. (Inf, 1.0) gave (0.0, NaN). All other special values (±0, ±1, ±Inf, NaN in both arguments) give the same results as before, including (NaN, NaN) at (0, 0).
• Performance: for |x| ≥ 1, the tanh derivative now evaluates sech(x) instead of reusing Ω. A broadcast over 10⁶ values with |x| up to about 5 took ~25% longer (5.8 ms → 7.1 ms).
• Performance: real asinh is unchanged (0.45 ms → 0.47 ms for 10⁶ values). atan(y, x) takes three divisions instead of two (0.37 ms → 0.96 ms for 10⁶ values), which is small compared to the primal atan(y, x) (~17 ms).
• The tanh rule now branches on abs(real(x)) < 1, and the real asinh rule on isfinite(x ^ 2), so they require Bool-valued comparisons. Number types whose comparisons return symbolic or traced values would fail here.
• FastMath variants: LLVM may reassociate (1 - x) * (1 + x) into 1 - x * x, so asin_fast etc. may not gain accuracy. The new tests pass for them.

Tests

FiniteDifferences can't detect these errors, so the new testsets compare frule against the analytic derivative evaluated in BigFloat, with rtol = 4eps().

🤖 Generated with Claude Code

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 96.65%. Comparing base (8f6cea7) to head (a2e18f7).

Additional details and impacted files

@@           Coverage Diff           @@
##             main     #864   +/-   ##
=======================================
  Coverage   96.64%   96.65%           
=======================================
  Files          28       28           
  Lines        3609     3614    +5     
=======================================
+ Hits         3488     3493    +5     
  Misses        121      121           

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