From 7a6fa443a273483cf597993d9e45a7fc9bb0d16c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Tue, 10 Nov 2026 02:27:27 +0100 Subject: [PATCH 1/3] Avoid catastrophic cancellation in tanh and inverse trig/hyperbolic derivatives MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `1 - tanh(x)^2` returns 0 once `tanh(x)` rounds to ±1 (#102); use the branch-free `sech(x)^2`. Rules involving `1 - x^2` or `x^2 - 1` lose accuracy near `|x| = 1`; use the factorized forms `(1 - x) * (1 + x)` and `(x - 1) * (x + 1)` instead, as in JuliaDiff/ChainRules.jl#864. Co-Authored-By: Claude Opus 5.5 --- src/rules.jl | 38 +++++++++++++++++++------------------- test/runtests.jl | 38 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 57 insertions(+), 19 deletions(-) diff --git a/src/rules.jl b/src/rules.jl index caa3078..8a69428 100644 --- a/src/rules.jl +++ b/src/rules.jl @@ -12,8 +12,8 @@ @define_diffrule Base.abs2(x) = :( $x + $x ) @define_diffrule Base.inv(x) = :( -abs2(inv($x)) ) @define_diffrule Base.log(x) = :( inv($x) ) -@define_diffrule Base.log10(x) = :( inv($x) / $logten ) -@define_diffrule Base.log2(x) = :( inv($x) / $logtwo ) +@define_diffrule Base.log10(x) = :( inv($logten * $x) ) +@define_diffrule Base.log2(x) = :( inv($logtwo * $x) ) @define_diffrule Base.log1p(x) = :( inv($x + 1) ) @define_diffrule Base.exp(x) = :( exp($x) ) @define_diffrule Base.exp2(x) = :( exp2($x) * $logtwo ) @@ -33,30 +33,30 @@ @define_diffrule Base.cotd(x) = :( - deg2rad(1 + cotd($x)^2) ) @define_diffrule Base.sinpi(x) = :( π * cospi($x) ) @define_diffrule Base.cospi(x) = :( -(π * sinpi($x)) ) -@define_diffrule Base.asin(x) = :( inv(sqrt(1 - $x^2)) ) -@define_diffrule Base.acos(x) = :( -inv(sqrt(1 - $x^2)) ) +@define_diffrule Base.asin(x) = :( inv(sqrt((1 - $x) * (1 + $x))) ) +@define_diffrule Base.acos(x) = :( -inv(sqrt((1 - $x) * (1 + $x))) ) @define_diffrule Base.atan(x) = :( inv(1 + $x^2) ) -@define_diffrule Base.asec(x) = :( inv(abs($x) * sqrt($x^2 - 1)) ) -@define_diffrule Base.acsc(x) = :( -inv(abs($x) * sqrt($x^2 - 1)) ) +@define_diffrule Base.asec(x) = :( inv(abs($x) * sqrt(($x - 1) * ($x + 1))) ) +@define_diffrule Base.acsc(x) = :( -inv(abs($x) * sqrt(($x - 1) * ($x + 1))) ) @define_diffrule Base.acot(x) = :( -inv(1 + $x^2) ) -@define_diffrule Base.asind(x) = :( inv(deg2rad(sqrt(1 - $x^2))) ) -@define_diffrule Base.acosd(x) = :( -inv(deg2rad(sqrt(1 - $x^2))) ) +@define_diffrule Base.asind(x) = :( inv(deg2rad(sqrt((1 - $x) * (1 + $x)))) ) +@define_diffrule Base.acosd(x) = :( -inv(deg2rad(sqrt((1 - $x) * (1 + $x)))) ) @define_diffrule Base.atand(x) = :( inv(deg2rad(1 + $x^2)) ) -@define_diffrule Base.asecd(x) = :( inv(deg2rad(abs($x) * sqrt($x^2 - 1))) ) -@define_diffrule Base.acscd(x) = :( -inv(deg2rad(abs($x) * sqrt($x^2 - 1))) ) +@define_diffrule Base.asecd(x) = :( inv(deg2rad(abs($x) * sqrt(($x - 1) * ($x + 1)))) ) +@define_diffrule Base.acscd(x) = :( -inv(deg2rad(abs($x) * sqrt(($x - 1) * ($x + 1)))) ) @define_diffrule Base.acotd(x) = :( -inv(deg2rad(1 + $x^2)) ) @define_diffrule Base.sinh(x) = :( cosh($x) ) @define_diffrule Base.cosh(x) = :( sinh($x) ) -@define_diffrule Base.tanh(x) = :( 1 - tanh($x)^2 ) +@define_diffrule Base.tanh(x) = :( sech($x)^2 ) @define_diffrule Base.sech(x) = :( -tanh($x) * sech($x) ) @define_diffrule Base.csch(x) = :( -coth($x) * csch($x) ) @define_diffrule Base.coth(x) = :( -(csch($x)^2) ) @define_diffrule Base.asinh(x) = :( inv(sqrt($x^2 + 1)) ) -@define_diffrule Base.acosh(x) = :( inv(sqrt($x^2 - 1)) ) -@define_diffrule Base.atanh(x) = :( inv(1 - $x^2) ) -@define_diffrule Base.asech(x) = :( -inv($x * sqrt(1 - $x^2)) ) +@define_diffrule Base.acosh(x) = :( inv(sqrt($x - 1) * sqrt($x + 1)) ) +@define_diffrule Base.atanh(x) = :( inv((1 - $x) * (1 + $x)) ) +@define_diffrule Base.asech(x) = :( -inv($x * sqrt((1 - $x) * (1 + $x))) ) @define_diffrule Base.acsch(x) = :( -inv(abs($x) * sqrt(1 + $x^2)) ) -@define_diffrule Base.acoth(x) = :( inv(1 - $x^2) ) +@define_diffrule Base.acoth(x) = :( inv((1 - $x) * (1 + $x)) ) @define_diffrule Base.sinc(x) = :( cosc($x) ) @define_diffrule Base.deg2rad(x) = :( deg2rad(one($x)) ) @define_diffrule Base.mod2pi(x) = :( isinteger($x / $twoπ) ? oftype(float($x), NaN) : one(float($x)) ) @@ -265,10 +265,10 @@ _abs_deriv(x) = signbit(x) ? -one(x) : one(x) @define_diffrule NaNMath.sin(x) = :( NaNMath.cos($x) ) @define_diffrule NaNMath.cos(x) = :( -NaNMath.sin($x) ) @define_diffrule NaNMath.tan(x) = :( 1 + NaNMath.pow(NaNMath.tan($x), 2) ) -@define_diffrule NaNMath.asin(x) = :( inv(NaNMath.sqrt(1 - NaNMath.pow($x, 2))) ) -@define_diffrule NaNMath.acos(x) = :( -inv(NaNMath.sqrt(1 - NaNMath.pow($x, 2))) ) -@define_diffrule NaNMath.acosh(x) = :( inv(NaNMath.sqrt(NaNMath.pow($x, 2) - 1)) ) -@define_diffrule NaNMath.atanh(x) = :( inv(1 - NaNMath.pow($x, 2)) ) +@define_diffrule NaNMath.asin(x) = :( inv(NaNMath.sqrt((1 - $x) * (1 + $x))) ) +@define_diffrule NaNMath.acos(x) = :( -inv(NaNMath.sqrt((1 - $x) * (1 + $x))) ) +@define_diffrule NaNMath.acosh(x) = :( inv(NaNMath.sqrt(($x - 1) * ($x + 1))) ) +@define_diffrule NaNMath.atanh(x) = :( inv((1 - $x) * (1 + $x)) ) @define_diffrule NaNMath.log(x) = :( inv($x) ) @define_diffrule NaNMath.log2(x) = :( inv($logtwo * $x) ) @define_diffrule NaNMath.log10(x) = :( inv($logten * $x) ) diff --git a/test/runtests.jl b/test/runtests.jl index f1fa041..4487639 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -186,6 +186,44 @@ for xytype in [:Float32, :Float64, :BigFloat] end end +# Rules must not suffer from catastrophic cancellation, e.g., `1 - x^2` near `|x| = 1` or +# `1 - tanh(x)^2` for large `|x|` (#102): compare with the rule evaluated in high precision +below_one(T) = (1 - 4 * eps(T), -(1 - 64 * eps(T))) # 1 ± x for |x| < 1 near 1 +above_one(T) = (1 + 4 * eps(T), -(1 + 64 * eps(T))) # x ± 1 for |x| > 1 near 1 +@testset "accuracy: ($M, $f)" for (M, f, xs) in [ + (:Base, :asin, below_one), + (:Base, :acos, below_one), + (:Base, :asind, below_one), + (:Base, :acosd, below_one), + (:Base, :atanh, below_one), + (:Base, :asech, T -> abs.(below_one(T))), + (:Base, :acoth, above_one), + (:Base, :acosh, T -> abs.(above_one(T))), + (:Base, :asec, above_one), + (:Base, :acsc, above_one), + (:Base, :asecd, above_one), + (:Base, :acscd, above_one), + (:Base, :tanh, T -> T.((0.5, -5, 10, -20, 30, 300))), + (:NaNMath, :asin, below_one), + (:NaNMath, :acos, below_one), + (:NaNMath, :atanh, below_one), + (:NaNMath, :acosh, T -> abs.(above_one(T))), +] + deriv = DiffRules.diffrule(M, f, :x) + @eval let + for T in (Float32, Float64), x in $xs(T) + # `x` is exact in `BigFloat`, so this is the derivative at `x` up to ~1e-70 + expected = setprecision(BigFloat, 256) do + let x = big(x) + T($deriv) + end + end + @test $deriv isa T + @test $deriv ≈ expected rtol = 8 * eps(T) + end + end +end + end @testset "diffrules" begin From 798d17f8d56c6f20975dfc690c4817c6f9117b69 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Tue, 10 Nov 2026 09:04:03 +0100 Subject: [PATCH 2/3] Avoid overflow in the derivative of NaNMath.acosh Use `sqrt(x - 1) * sqrt(x + 1)` as for `Base.acosh`: `(x - 1) * (x + 1)` overflows for large `x`, and the derivative is now NaN for `x < 1`, like the primal. Test both rules at a point where `x^2` overflows. Co-Authored-By: Claude Opus 5.5 --- src/rules.jl | 2 +- test/runtests.jl | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/rules.jl b/src/rules.jl index 8a69428..35d0c25 100644 --- a/src/rules.jl +++ b/src/rules.jl @@ -267,7 +267,7 @@ _abs_deriv(x) = signbit(x) ? -one(x) : one(x) @define_diffrule NaNMath.tan(x) = :( 1 + NaNMath.pow(NaNMath.tan($x), 2) ) @define_diffrule NaNMath.asin(x) = :( inv(NaNMath.sqrt((1 - $x) * (1 + $x))) ) @define_diffrule NaNMath.acos(x) = :( -inv(NaNMath.sqrt((1 - $x) * (1 + $x))) ) -@define_diffrule NaNMath.acosh(x) = :( inv(NaNMath.sqrt(($x - 1) * ($x + 1))) ) +@define_diffrule NaNMath.acosh(x) = :( inv(NaNMath.sqrt($x - 1) * NaNMath.sqrt($x + 1)) ) @define_diffrule NaNMath.atanh(x) = :( inv((1 - $x) * (1 + $x)) ) @define_diffrule NaNMath.log(x) = :( inv($x) ) @define_diffrule NaNMath.log2(x) = :( inv($logtwo * $x) ) diff --git a/test/runtests.jl b/test/runtests.jl index 4487639..8fc8994 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -190,6 +190,7 @@ end # `1 - tanh(x)^2` for large `|x|` (#102): compare with the rule evaluated in high precision below_one(T) = (1 - 4 * eps(T), -(1 - 64 * eps(T))) # 1 ± x for |x| < 1 near 1 above_one(T) = (1 + 4 * eps(T), -(1 + 64 * eps(T))) # x ± 1 for |x| > 1 near 1 +acosh_xs(T) = (abs.(above_one(T))..., 2 * sqrt(floatmax(T))) # x^2 overflows for the last @testset "accuracy: ($M, $f)" for (M, f, xs) in [ (:Base, :asin, below_one), (:Base, :acos, below_one), @@ -198,7 +199,7 @@ above_one(T) = (1 + 4 * eps(T), -(1 + 64 * eps(T))) # x ± 1 for |x| > 1 near (:Base, :atanh, below_one), (:Base, :asech, T -> abs.(below_one(T))), (:Base, :acoth, above_one), - (:Base, :acosh, T -> abs.(above_one(T))), + (:Base, :acosh, acosh_xs), (:Base, :asec, above_one), (:Base, :acsc, above_one), (:Base, :asecd, above_one), @@ -207,7 +208,7 @@ above_one(T) = (1 + 4 * eps(T), -(1 + 64 * eps(T))) # x ± 1 for |x| > 1 near (:NaNMath, :asin, below_one), (:NaNMath, :acos, below_one), (:NaNMath, :atanh, below_one), - (:NaNMath, :acosh, T -> abs.(above_one(T))), + (:NaNMath, :acosh, acosh_xs), ] deriv = DiffRules.diffrule(M, f, :x) @eval let From cf65762e7d491c95d6f286706faa6a68d178a0e7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Tue, 10 Nov 2026 09:12:46 +0100 Subject: [PATCH 3/3] Bump version to 1.17.0 Co-Authored-By: Claude Opus 5.5 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e30e835..c8d40c5 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "DiffRules" uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" -version = "1.17.0-dev" +version = "1.17.0" [deps] IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6"