From cbd429bdfcc843d7c6498bd0091c262059629023 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?David=20M=C3=BCller-Widmann?=
 <devmotion@users.noreply.github.com>
Date: Tue, 10 Nov 2026 02:28:02 +0100
Subject: [PATCH] Accurate derivative of tanh for large |x|
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

DiffRules' `1 - tanh(x)^2` returns 0 once `tanh(x)` rounds to ±1
(JuliaDiff/DiffRules.jl#102). Define `tanh(::Dual)` by hand: reuse the
primal value for `|x| < 1` and use `sech(x)^2` otherwise, as in
JuliaDiff/ChainRules.jl#864.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/dual.jl      | 14 +++++++++++++-
 test/DualTest.jl | 22 +++++++++++++++++++++-
 2 files changed, 34 insertions(+), 2 deletions(-)

diff --git a/src/dual.jl b/src/dual.jl
index 48d5d0e3..8d4f0bb3 100644
--- a/src/dual.jl
+++ b/src/dual.jl
@@ -474,7 +474,7 @@ Base.float(d::Dual) = convert(float(typeof(d)), d)
 ###################################
 
 for (M, f, arity) in DiffRules.diffrules(filter_modules = nothing)
-    if (M, f) in ((:Base, :^), (:NaNMath, :pow), (:Base, :/), (:Base, :+), (:Base, :-), (:Base, :sin), (:Base, :cos))
+    if (M, f) in ((:Base, :^), (:NaNMath, :pow), (:Base, :/), (:Base, :+), (:Base, :-), (:Base, :sin), (:Base, :cos), (:Base, :tanh))
         continue  # Skip methods which we define elsewhere.
     elseif !(isdefined(@__MODULE__, M) && isdefined(getfield(@__MODULE__, M), f))
         continue  # Skip rules for methods not defined in the current scope
@@ -721,6 +721,18 @@ end
     return (Dual{T}(sd, cd * partials(d)), Dual{T}(cd, -sd * partials(d)))
 end
 
+# tanh #
+#------#
+
+# `1 - tanh(x)^2` reuses the primal value but suffers from catastrophic cancellation for
+# large `|x|` (JuliaDiff/DiffRules.jl#102), so use `sech(x)^2` there
+@inline function Base.tanh(d::Dual{T}) where T
+    x = value(d)
+    val = tanh(x)
+    deriv = abs(x) < 1 ? 1 - val^2 : sech(x)^2
+    return Dual{T}(val, deriv * partials(d))
+end
+
 # sincospi #
 #----------#
 
diff --git a/test/DualTest.jl b/test/DualTest.jl
index e7429be5..ba151207 100644
--- a/test/DualTest.jl
+++ b/test/DualTest.jl
@@ -559,7 +559,12 @@ ForwardDiff.:≺(::Type{OuterTestTag}, ::Type{TestTag}) = false
                     if actualval isa Real
                         @test dx isa Dual{TestTag}
                         @test value(dx) == actualval
-                        @test partials(dx, 1) == $deriv
+                        if $(M === :Base && f === :tanh)
+                            # defined by hand with an equivalent formula that is faster for |x| < 1
+                            @test partials(dx, 1) ≈ $deriv
+                        else
+                            @test partials(dx, 1) == $deriv
+                        end
                     else
                         @test dx isa Complex{<:Dual{TestTag}}
                         @test value(real(dx)) == real(actualval)
@@ -775,4 +780,19 @@ end
     end
 end
 
+@testset "tanh: no catastrophic cancellation for large |x|" begin # JuliaDiff/DiffRules.jl#102
+    for T in (Float32, Float64), x in (0.5, -0.99, 1.0, -5.0, 10.0, -20.0, 30.0, 300.0)
+        x = T(x)
+        expected = T(setprecision(() -> sech(big(x))^2, BigFloat, 256))
+        d = tanh(Dual{TestTag}(x, one(x)))
+        @test value(d) === tanh(x)
+        @test partials(d, 1) isa T
+        @test partials(d, 1) ≈ expected rtol = 4 * eps(T)
+        # nested duals take the same branch
+        dd = tanh(Dual{TestTag}(Dual{TestTag}(x, one(x)), one(x)))
+        @test value(value(dd)) === tanh(x)
+        @test value(partials(dd, 1)) ≈ expected rtol = 4 * eps(T)
+    end
+end
+
 end # module