From 53c46efbac2947884822d15af9f3119541b0697f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Thu, 5 Nov 2026 15:02:54 +0100 Subject: [PATCH 1/7] 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 The `tanh` rule `1 - Ω ^ 2` loses all accuracy for large `|x|` since `Ω = tanh(x)` rounds to ±1 (see JuliaDiff/DiffRules.jl#102); use `sech(x) ^ 2` for `|real(x)| ≥ 1`. Rules involving `1 - x ^ 2` or `x ^ 2 - 1` (asin, acos, atanh, acoth, asech, asind, acosd, and the real variants of asec, acsc, asecd, acscd) lose accuracy near `|x| = 1`; use the factorized forms `(1 - x) * (1 + x)` and `(x - 1) * (x + 1)` instead. Co-Authored-By: Claude Opus 5.5 --- Project.toml | 2 +- src/rulesets/Base/base.jl | 19 ++++++++++--------- src/rulesets/Base/fastmath_able.jl | 7 ++++--- test/rulesets/Base/base.jl | 22 ++++++++++++++++++++++ test/rulesets/Base/fastmath_able.jl | 15 +++++++++++++++ 5 files changed, 52 insertions(+), 13 deletions(-) diff --git a/Project.toml b/Project.toml index 952b846d4..65c7177b9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ChainRules" uuid = "082447d4-558c-5d27-93f4-14fc19e9eca2" -version = "1.74.10" +version = "1.74.11" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/src/rulesets/Base/base.jl b/src/rulesets/Base/base.jl index 357a3cd5e..25940c3ac 100644 --- a/src/rulesets/Base/base.jl +++ b/src/rulesets/Base/base.jl @@ -168,21 +168,22 @@ end # Can't multiply though sqrt in acosh because of negative complex case for x @scalar_rule acosh(x) inv(sqrt(x - 1) * sqrt(x + 1)) -@scalar_rule acoth(x) inv(1 - x ^ 2) +# Factorized `(1 - x) * (1 + x)` avoids cancellation for `|x|` close to 1 +@scalar_rule acoth(x) inv((1 - x) * (1 + x)) @scalar_rule acsch(x) -(inv(x ^ 2 * sqrt(1 + x ^ -2))) @scalar_rule acsch(x::Real) -(inv(abs(x) * sqrt(1 + x ^ 2))) -@scalar_rule asech(x) -(inv(x * sqrt(1 - x ^ 2))) +@scalar_rule asech(x) -inv(x * sqrt((1 - x) * (1 + x))) @scalar_rule asinh(x) inv(sqrt(x ^ 2 + 1)) -@scalar_rule atanh(x) inv(1 - x ^ 2) +@scalar_rule atanh(x) inv((1 - x) * (1 + x)) -@scalar_rule acosd(x) -inv(deg2rad(sqrt(1 - x ^ 2))) +@scalar_rule acosd(x) -inv(deg2rad(sqrt((1 - x) * (1 + x)))) @scalar_rule acotd(x) -inv(deg2rad(1 + x ^ 2)) @scalar_rule acscd(x) -inv(deg2rad(x^2 * sqrt(1 - x ^ -2))) -@scalar_rule acscd(x::Real) -inv(deg2rad(abs(x) * sqrt(x ^ 2 - 1))) +@scalar_rule acscd(x::Real) -inv(deg2rad(abs(x) * sqrt((x - 1) * (x + 1)))) @scalar_rule asecd(x) inv(deg2rad(x ^ 2 * sqrt(1 - x ^ -2))) -@scalar_rule asecd(x::Real) inv(deg2rad(abs(x) * sqrt(x ^ 2 - 1))) -@scalar_rule asind(x) inv(deg2rad(sqrt(1 - x ^ 2))) +@scalar_rule asecd(x::Real) inv(deg2rad(abs(x) * sqrt((x - 1) * (x + 1)))) +@scalar_rule asind(x) inv(deg2rad(sqrt((1 - x) * (1 + x)))) @scalar_rule atand(x) inv(deg2rad(1 + x ^ 2)) @scalar_rule cot(x) -((1 + Ω ^ 2)) @@ -197,9 +198,9 @@ end @scalar_rule acot(x) -(inv(1 + x ^ 2)) @scalar_rule acsc(x) -(inv(x ^ 2 * sqrt(1 - x ^ -2))) -@scalar_rule acsc(x::Real) -(inv(abs(x) * sqrt(x ^ 2 - 1))) +@scalar_rule acsc(x::Real) -inv(abs(x) * sqrt((x - 1) * (x + 1))) @scalar_rule asec(x) inv(x ^ 2 * sqrt(1 - x ^ -2)) -@scalar_rule asec(x::Real) inv(abs(x) * sqrt(x ^ 2 - 1)) +@scalar_rule asec(x::Real) inv(abs(x) * sqrt((x - 1) * (x + 1))) @scalar_rule cosd(x) -deg2rad(sind(x)) @scalar_rule cospi(x) -π * sinpi(x) diff --git a/src/rulesets/Base/fastmath_able.jl b/src/rulesets/Base/fastmath_able.jl index e8473661f..44b77bf9f 100644 --- a/src/rulesets/Base/fastmath_able.jl +++ b/src/rulesets/Base/fastmath_able.jl @@ -38,11 +38,12 @@ let # Trig-Hyperbolic @scalar_rule cosh(x) sinh(x) @scalar_rule sinh(x) cosh(x) - @scalar_rule tanh(x) 1 - Ω ^ 2 + # `1 - Ω ^ 2` cancels catastrophically for large `|x|` (DiffRules.jl#102) + @scalar_rule tanh(x) abs(real(x)) < 1 ? 1 - Ω ^ 2 : sech(x) ^ 2 # Trig- Inverses - @scalar_rule acos(x) -(inv(sqrt(1 - x ^ 2))) - @scalar_rule asin(x) inv(sqrt(1 - x ^ 2)) + @scalar_rule acos(x) -inv(sqrt((1 - x) * (1 + x))) + @scalar_rule asin(x) inv(sqrt((1 - x) * (1 + x))) @scalar_rule atan(x) inv(1 + x ^ 2) # Trig-Multivariate diff --git a/test/rulesets/Base/base.jl b/test/rulesets/Base/base.jl index e9d0fa3ef..89eefb4d5 100644 --- a/test/rulesets/Base/base.jl +++ b/test/rulesets/Base/base.jl @@ -67,6 +67,28 @@ end test_scalar(acoth, x + 1) end + # finite differencing is not accurate enough to detect cancellation near |x| = 1 + @testset "Inverse accuracy near |x| = 1" begin + setprecision(BigFloat, 256) do + for (f, df, xs) in ( + (atanh, x -> inv(1 - x^2), (1 - 1e-8, -1 + 1e-8, ComplexF64(1, 1e-9))), + (acoth, x -> inv(1 - x^2), (1 + 1e-8, -1 - 1e-8, ComplexF64(1, 1e-9))), + (asech, x -> -inv(x * sqrt(1 - x^2)), (1 - 1e-8, ComplexF64(1, 1e-9))), + (asind, x -> inv(deg2rad(sqrt(1 - x^2))), (1 - 1e-8, -1 + 1e-8)), + (acosd, x -> -inv(deg2rad(sqrt(1 - x^2))), (1 - 1e-8, -1 + 1e-8)), + (asecd, x -> inv(deg2rad(abs(x) * sqrt(x^2 - 1))), (1 + 1e-8, -1 - 1e-8)), + (acscd, x -> -inv(deg2rad(abs(x) * sqrt(x^2 - 1))), (1 + 1e-8, -1 - 1e-8)), + (asec, x -> inv(abs(x) * sqrt(x^2 - 1)), (1 + 1e-8, -1 - 1e-8)), + (acsc, x -> -inv(abs(x) * sqrt(x^2 - 1)), (1 + 1e-8, -1 - 1e-8)), + ) + @testset "$f($x)" for x in xs + _, ẏ = frule((NoTangent(), one(x)), f, x) + @test ẏ ≈ df(big(x)) rtol = 4eps() + end + end + end + end + @testset "Inverse degrees" for x = (0.5, Complex(0.5, 0.25)) test_scalar(asind, x) test_scalar(acosd, x) diff --git a/test/rulesets/Base/fastmath_able.jl b/test/rulesets/Base/fastmath_able.jl index c411811a2..2c5cf48ae 100644 --- a/test/rulesets/Base/fastmath_able.jl +++ b/test/rulesets/Base/fastmath_able.jl @@ -59,6 +59,21 @@ const FASTABLE_AST = quote test_scalar(acos, x) test_scalar(atan, x) end + # finite differencing is not accurate enough to detect cancellation + @testset "Accuracy" begin + setprecision(BigFloat, 256) do + for (f, df, xs) in ( + (tanh, x -> sech(x)^2, (0.5, 4.0, 14.0, 25.0, -25.0, ComplexF64(14, 0.3), ComplexF64(30, 1))), + (asin, x -> inv(sqrt(1 - x^2)), (1 - 1e-8, -1 + 1e-8, ComplexF64(1, 1e-9))), + (acos, x -> -inv(sqrt(1 - x^2)), (1 - 1e-8, -1 + 1e-8, ComplexF64(1, 1e-9))), + ) + @testset "$f($x)" for x in xs + _, ẏ = frule((NoTangent(), one(x)), f, x) + @test ẏ ≈ df(big(x)) rtol = 4eps() + end + end + end + end @testset "Multivariate" begin @testset "sincos(x::$T)" for T in (Float64, ComplexF64) Δz = Tangent{Tuple{T,T}}(randn(T), randn(T)) From 6cd5e3e3365c57631651ed11ca75531eece6d44e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Thu, 5 Nov 2026 15:25:29 +0100 Subject: [PATCH 2/7] Apply batched suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- test/rulesets/Base/base.jl | 12 ++++++++++-- test/rulesets/Base/fastmath_able.jl | 26 +++++++++++++++++++++++--- 2 files changed, 33 insertions(+), 5 deletions(-) diff --git a/test/rulesets/Base/base.jl b/test/rulesets/Base/base.jl index 89eefb4d5..8f7317fa9 100644 --- a/test/rulesets/Base/base.jl +++ b/test/rulesets/Base/base.jl @@ -76,8 +76,16 @@ end (asech, x -> -inv(x * sqrt(1 - x^2)), (1 - 1e-8, ComplexF64(1, 1e-9))), (asind, x -> inv(deg2rad(sqrt(1 - x^2))), (1 - 1e-8, -1 + 1e-8)), (acosd, x -> -inv(deg2rad(sqrt(1 - x^2))), (1 - 1e-8, -1 + 1e-8)), - (asecd, x -> inv(deg2rad(abs(x) * sqrt(x^2 - 1))), (1 + 1e-8, -1 - 1e-8)), - (acscd, x -> -inv(deg2rad(abs(x) * sqrt(x^2 - 1))), (1 + 1e-8, -1 - 1e-8)), + ( + asecd, + x -> inv(deg2rad(abs(x) * sqrt(x^2 - 1))), + (1 + 1e-8, -1 - 1e-8), + ), + ( + acscd, + x -> -inv(deg2rad(abs(x) * sqrt(x^2 - 1))), + (1 + 1e-8, -1 - 1e-8), + ), (asec, x -> inv(abs(x) * sqrt(x^2 - 1)), (1 + 1e-8, -1 - 1e-8)), (acsc, x -> -inv(abs(x) * sqrt(x^2 - 1)), (1 + 1e-8, -1 - 1e-8)), ) diff --git a/test/rulesets/Base/fastmath_able.jl b/test/rulesets/Base/fastmath_able.jl index 2c5cf48ae..797d91ef2 100644 --- a/test/rulesets/Base/fastmath_able.jl +++ b/test/rulesets/Base/fastmath_able.jl @@ -63,9 +63,29 @@ const FASTABLE_AST = quote @testset "Accuracy" begin setprecision(BigFloat, 256) do for (f, df, xs) in ( - (tanh, x -> sech(x)^2, (0.5, 4.0, 14.0, 25.0, -25.0, ComplexF64(14, 0.3), ComplexF64(30, 1))), - (asin, x -> inv(sqrt(1 - x^2)), (1 - 1e-8, -1 + 1e-8, ComplexF64(1, 1e-9))), - (acos, x -> -inv(sqrt(1 - x^2)), (1 - 1e-8, -1 + 1e-8, ComplexF64(1, 1e-9))), + ( + tanh, + x -> sech(x)^2, + ( + 0.5, + 4.0, + 14.0, + 25.0, + -25.0, + ComplexF64(14, 0.3), + ComplexF64(30, 1), + ), + ), + ( + asin, + x -> inv(sqrt(1 - x^2)), + (1 - 1e-8, -1 + 1e-8, ComplexF64(1, 1e-9)), + ), + ( + acos, + x -> -inv(sqrt(1 - x^2)), + (1 - 1e-8, -1 + 1e-8, ComplexF64(1, 1e-9)), + ), ) @testset "$f($x)" for x in xs _, ẏ = frule((NoTangent(), one(x)), f, x) From f31f108e4420b1bf06f1f9e3eeec6bbc06bf9b85 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Fri, 6 Nov 2026 00:17:41 +0100 Subject: [PATCH 3/7] Cite reference for the `(1 - x) * (1 + x)` factorization Co-Authored-By: Claude Opus 5.5 --- src/rulesets/Base/base.jl | 3 ++- src/rulesets/Base/fastmath_able.jl | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/rulesets/Base/base.jl b/src/rulesets/Base/base.jl index 25940c3ac..20cdd7bb6 100644 --- a/src/rulesets/Base/base.jl +++ b/src/rulesets/Base/base.jl @@ -168,7 +168,8 @@ end # Can't multiply though sqrt in acosh because of negative complex case for x @scalar_rule acosh(x) inv(sqrt(x - 1) * sqrt(x + 1)) -# Factorized `(1 - x) * (1 + x)` avoids cancellation for `|x|` close to 1 +# `(1 - x) * (1 + x)` avoids cancellation in `1 - x ^ 2` for `|x|` close to 1 +# (Goldberg (1991), "What Every Computer Scientist Should Know About Floating-Point Arithmetic") @scalar_rule acoth(x) inv((1 - x) * (1 + x)) @scalar_rule acsch(x) -(inv(x ^ 2 * sqrt(1 + x ^ -2))) @scalar_rule acsch(x::Real) -(inv(abs(x) * sqrt(1 + x ^ 2))) diff --git a/src/rulesets/Base/fastmath_able.jl b/src/rulesets/Base/fastmath_able.jl index 44b77bf9f..4044747d7 100644 --- a/src/rulesets/Base/fastmath_able.jl +++ b/src/rulesets/Base/fastmath_able.jl @@ -42,6 +42,7 @@ let @scalar_rule tanh(x) abs(real(x)) < 1 ? 1 - Ω ^ 2 : sech(x) ^ 2 # Trig- Inverses + # `(1 - x) * (1 + x)` avoids cancellation in `1 - x ^ 2` for `|x|` close to 1 (see `acoth` in base.jl) @scalar_rule acos(x) -inv(sqrt((1 - x) * (1 + x))) @scalar_rule asin(x) inv(sqrt((1 - x) * (1 + x))) @scalar_rule atan(x) inv(1 + x ^ 2) From e59d0bf1840e727ba1bceb330795a14b3879fa64 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Fri, 6 Nov 2026 01:06:14 +0100 Subject: [PATCH 4/7] =?UTF-8?q?Avoid=20cancellation=20near=20=C2=B1i=20and?= =?UTF-8?q?=20overflow=20in=20asinh,=20atan,=20acot,=20atand,=20acotd?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Complex `1 + x ^ 2` is factorized as `(1 + im * x) * (1 - im * x)` - Real `asinh` uses `hypot(x, 1)`, `atan(y, x)` uses `hypot(x, y)` Co-Authored-By: Claude Opus 5.5 --- src/rulesets/Base/base.jl | 9 +++++++++ src/rulesets/Base/fastmath_able.jl | 5 ++++- test/rulesets/Base/base.jl | 18 ++++++++++++++++++ test/rulesets/Base/fastmath_able.jl | 18 ++++++++++++++++++ 4 files changed, 49 insertions(+), 1 deletion(-) diff --git a/src/rulesets/Base/base.jl b/src/rulesets/Base/base.jl index 20cdd7bb6..ff0012430 100644 --- a/src/rulesets/Base/base.jl +++ b/src/rulesets/Base/base.jl @@ -175,17 +175,25 @@ end @scalar_rule acsch(x::Real) -(inv(abs(x) * sqrt(1 + x ^ 2))) @scalar_rule asech(x) -inv(x * sqrt((1 - x) * (1 + x))) @scalar_rule asinh(x) inv(sqrt(x ^ 2 + 1)) +# `hypot` avoids overflow of `x ^ 2` for large `|x|` +@scalar_rule asinh(x::Real) inv(hypot(x, 1)) +# `(1 + ix) * (1 - ix)` avoids cancellation in `1 + x ^ 2` near `x = ±i`; +# `ix` is constructed explicitly since `im * x` loses the sign of zero +@scalar_rule asinh(x::Complex) @setup(ix = Complex(-imag(x), real(x))) inv(sqrt((1 + ix) * (1 - ix))) @scalar_rule atanh(x) inv((1 - x) * (1 + x)) @scalar_rule acosd(x) -inv(deg2rad(sqrt((1 - x) * (1 + x)))) @scalar_rule acotd(x) -inv(deg2rad(1 + x ^ 2)) +# `(1 + im * x) * (1 - im * x)` avoids cancellation in `1 + x ^ 2` near `x = ±i` +@scalar_rule acotd(x::Complex) -inv(deg2rad((1 + im * x) * (1 - im * x))) @scalar_rule acscd(x) -inv(deg2rad(x^2 * sqrt(1 - x ^ -2))) @scalar_rule acscd(x::Real) -inv(deg2rad(abs(x) * sqrt((x - 1) * (x + 1)))) @scalar_rule asecd(x) inv(deg2rad(x ^ 2 * sqrt(1 - x ^ -2))) @scalar_rule asecd(x::Real) inv(deg2rad(abs(x) * sqrt((x - 1) * (x + 1)))) @scalar_rule asind(x) inv(deg2rad(sqrt((1 - x) * (1 + x)))) @scalar_rule atand(x) inv(deg2rad(1 + x ^ 2)) +@scalar_rule atand(x::Complex) inv(deg2rad((1 + im * x) * (1 - im * x))) @scalar_rule cot(x) -((1 + Ω ^ 2)) @scalar_rule coth(x) -(csch(x) ^ 2) @@ -198,6 +206,7 @@ end @scalar_rule sech(x) -(tanh(x)) * Ω @scalar_rule acot(x) -(inv(1 + x ^ 2)) +@scalar_rule acot(x::Complex) -inv((1 + im * x) * (1 - im * x)) @scalar_rule acsc(x) -(inv(x ^ 2 * sqrt(1 - x ^ -2))) @scalar_rule acsc(x::Real) -inv(abs(x) * sqrt((x - 1) * (x + 1))) @scalar_rule asec(x) inv(x ^ 2 * sqrt(1 - x ^ -2)) diff --git a/src/rulesets/Base/fastmath_able.jl b/src/rulesets/Base/fastmath_able.jl index 4044747d7..a28f87de4 100644 --- a/src/rulesets/Base/fastmath_able.jl +++ b/src/rulesets/Base/fastmath_able.jl @@ -46,9 +46,12 @@ let @scalar_rule acos(x) -inv(sqrt((1 - x) * (1 + x))) @scalar_rule asin(x) inv(sqrt((1 - x) * (1 + x))) @scalar_rule atan(x) inv(1 + x ^ 2) + # `(1 + im * x) * (1 - im * x)` avoids cancellation in `1 + x ^ 2` near `x = ±i` + @scalar_rule atan(x::Complex) inv((1 + im * x) * (1 - im * x)) # Trig-Multivariate - @scalar_rule atan(y, x) @setup(u = x ^ 2 + y ^ 2) (x / u, -y / u) + # `hypot` avoids overflow and underflow of `x ^ 2 + y ^ 2` + @scalar_rule atan(y, x) @setup(h = hypot(x, y)) (x / h / h, -y / h / h) @scalar_rule sincos(x) @setup((sinx, cosx) = Ω) cosx -sinx # exponents diff --git a/test/rulesets/Base/base.jl b/test/rulesets/Base/base.jl index 8f7317fa9..8c2e807eb 100644 --- a/test/rulesets/Base/base.jl +++ b/test/rulesets/Base/base.jl @@ -97,6 +97,24 @@ end end end + # finite differencing is not accurate enough to detect cancellation near x = ±i + @testset "Inverse accuracy near x = ±i and for large |x|" begin + setprecision(BigFloat, 256) do + zs = (ComplexF64(1e-9, 1 - 1e-8), ComplexF64(-1e-9, -1 + 1e-8)) + for (f, df, xs) in ( + (asinh, x -> inv(sqrt(1 + x^2)), (zs..., 1e200, -1e200)), + (acot, x -> -inv(1 + x^2), zs), + (atand, x -> inv(deg2rad(1 + x^2)), zs), + (acotd, x -> -inv(deg2rad(1 + x^2)), zs), + ) + @testset "$f($x)" for x in xs + _, ẏ = frule((NoTangent(), one(x)), f, x) + @test ẏ ≈ df(big(x)) rtol = 4eps() + end + end + end + end + @testset "Inverse degrees" for x = (0.5, Complex(0.5, 0.25)) test_scalar(asind, x) test_scalar(acosd, x) diff --git a/test/rulesets/Base/fastmath_able.jl b/test/rulesets/Base/fastmath_able.jl index 797d91ef2..bcf850593 100644 --- a/test/rulesets/Base/fastmath_able.jl +++ b/test/rulesets/Base/fastmath_able.jl @@ -86,6 +86,11 @@ const FASTABLE_AST = quote x -> -inv(sqrt(1 - x^2)), (1 - 1e-8, -1 + 1e-8, ComplexF64(1, 1e-9)), ), + ( + atan, + x -> inv(1 + x^2), + (ComplexF64(1e-9, 1 - 1e-8), ComplexF64(-1e-9, -1 + 1e-8)), + ), ) @testset "$f($x)" for x in xs _, ẏ = frule((NoTangent(), one(x)), f, x) @@ -334,3 +339,16 @@ end @eval @testset "fastmath_able FastMath functions" begin $(Base.FastMath.make_fastmath(FASTABLE_AST)) end + +# Not part of `FASTABLE_AST`: `hypot_fast` does not guard against overflow and underflow +@testset "atan(y, x) for large and small arguments" begin + setprecision(BigFloat, 256) do + @testset "atan($y, $x)" for (y, x) in ((1e-200, 1e-200), (1e200, 3e200), (-1e-170, 3e-170)) + _, pb = rrule(atan, y, x) + _, ȳ, x̄ = pb(1.0) + u = big(x)^2 + big(y)^2 + @test ȳ ≈ big(x) / u rtol = 4eps() + @test x̄ ≈ -big(y) / u rtol = 4eps() + end + end +end From 20403eb8d18d4036d10455fee5ed173b5305f2d6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Fri, 6 Nov 2026 01:10:23 +0100 Subject: [PATCH 5/7] Apply batched suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/rulesets/Base/base.jl | 4 +++- src/rulesets/Base/fastmath_able.jl | 4 ++-- test/rulesets/Base/fastmath_able.jl | 3 ++- 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/rulesets/Base/base.jl b/src/rulesets/Base/base.jl index ff0012430..d8204d8bd 100644 --- a/src/rulesets/Base/base.jl +++ b/src/rulesets/Base/base.jl @@ -179,7 +179,9 @@ end @scalar_rule asinh(x::Real) inv(hypot(x, 1)) # `(1 + ix) * (1 - ix)` avoids cancellation in `1 + x ^ 2` near `x = ±i`; # `ix` is constructed explicitly since `im * x` loses the sign of zero -@scalar_rule asinh(x::Complex) @setup(ix = Complex(-imag(x), real(x))) inv(sqrt((1 + ix) * (1 - ix))) +@scalar_rule asinh(x::Complex) @setup(ix=Complex(-imag(x), real(x))) inv( + sqrt((1 + ix) * (1 - ix)) +) @scalar_rule atanh(x) inv((1 - x) * (1 + x)) diff --git a/src/rulesets/Base/fastmath_able.jl b/src/rulesets/Base/fastmath_able.jl index a28f87de4..9f04e8844 100644 --- a/src/rulesets/Base/fastmath_able.jl +++ b/src/rulesets/Base/fastmath_able.jl @@ -51,8 +51,8 @@ let # Trig-Multivariate # `hypot` avoids overflow and underflow of `x ^ 2 + y ^ 2` - @scalar_rule atan(y, x) @setup(h = hypot(x, y)) (x / h / h, -y / h / h) - @scalar_rule sincos(x) @setup((sinx, cosx) = Ω) cosx -sinx + @scalar_rule atan(y, x) @setup(h=hypot(x, y)) (x / h / h, -y / h / h) + @scalar_rule sincos(x) @setup((sinx, cosx)=Ω) cosx -sinx # exponents @scalar_rule cbrt(x) inv(3 * Ω ^ 2) diff --git a/test/rulesets/Base/fastmath_able.jl b/test/rulesets/Base/fastmath_able.jl index bcf850593..5cf2a9b1a 100644 --- a/test/rulesets/Base/fastmath_able.jl +++ b/test/rulesets/Base/fastmath_able.jl @@ -343,7 +343,8 @@ end # Not part of `FASTABLE_AST`: `hypot_fast` does not guard against overflow and underflow @testset "atan(y, x) for large and small arguments" begin setprecision(BigFloat, 256) do - @testset "atan($y, $x)" for (y, x) in ((1e-200, 1e-200), (1e200, 3e200), (-1e-170, 3e-170)) + @testset "atan($y, $x)" for (y, x) in + ((1e-200, 1e-200), (1e200, 3e200), (-1e-170, 3e-170)) _, pb = rrule(atan, y, x) _, ȳ, x̄ = pb(1.0) u = big(x)^2 + big(y)^2 From 26ba16f26a6cdd11038c808e7f644bcbbebe176e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Fri, 6 Nov 2026 12:25:15 +0100 Subject: [PATCH 6/7] Use cheap overflow guards instead of hypot in asinh and atan(y, x) `hypot` prevents SIMD and made the derivatives ~10x slower. `asinh` now falls back to `1 / |x|` only if `x ^ 2` overflows, and `atan(y, x)` scales by `max(|x|, |y|)`, which also fixes `atan_fast`. Co-Authored-By: Claude Opus 5.5 --- src/rulesets/Base/base.jl | 4 ++-- src/rulesets/Base/fastmath_able.jl | 6 ++++-- test/rulesets/Base/fastmath_able.jl | 8 ++++---- 3 files changed, 10 insertions(+), 8 deletions(-) diff --git a/src/rulesets/Base/base.jl b/src/rulesets/Base/base.jl index d8204d8bd..8a50a7a05 100644 --- a/src/rulesets/Base/base.jl +++ b/src/rulesets/Base/base.jl @@ -175,8 +175,8 @@ end @scalar_rule acsch(x::Real) -(inv(abs(x) * sqrt(1 + x ^ 2))) @scalar_rule asech(x) -inv(x * sqrt((1 - x) * (1 + x))) @scalar_rule asinh(x) inv(sqrt(x ^ 2 + 1)) -# `hypot` avoids overflow of `x ^ 2` for large `|x|` -@scalar_rule asinh(x::Real) inv(hypot(x, 1)) +# if `x ^ 2` overflows, `1 / |x|` is exact up to rounding +@scalar_rule asinh(x::Real) (x2=x ^ 2; isfinite(x2) ? inv(sqrt(x2 + 1)) : inv(abs(x))) # `(1 + ix) * (1 - ix)` avoids cancellation in `1 + x ^ 2` near `x = ±i`; # `ix` is constructed explicitly since `im * x` loses the sign of zero @scalar_rule asinh(x::Complex) @setup(ix=Complex(-imag(x), real(x))) inv( diff --git a/src/rulesets/Base/fastmath_able.jl b/src/rulesets/Base/fastmath_able.jl index 9f04e8844..03cc7324e 100644 --- a/src/rulesets/Base/fastmath_able.jl +++ b/src/rulesets/Base/fastmath_able.jl @@ -50,8 +50,10 @@ let @scalar_rule atan(x::Complex) inv((1 + im * x) * (1 - im * x)) # Trig-Multivariate - # `hypot` avoids overflow and underflow of `x ^ 2 + y ^ 2` - @scalar_rule atan(y, x) @setup(h=hypot(x, y)) (x / h / h, -y / h / h) + # scaling by `max(|x|, |y|)` avoids overflow and underflow of `x ^ 2 + y ^ 2` + @scalar_rule atan(y, x) @setup( + s = max(abs(x), abs(y)), xs = x / s, ys = y / s, v = xs ^ 2 + ys ^ 2 + ) (xs / v / s, -ys / v / s) @scalar_rule sincos(x) @setup((sinx, cosx)=Ω) cosx -sinx # exponents diff --git a/test/rulesets/Base/fastmath_able.jl b/test/rulesets/Base/fastmath_able.jl index 5cf2a9b1a..070d3eb49 100644 --- a/test/rulesets/Base/fastmath_able.jl +++ b/test/rulesets/Base/fastmath_able.jl @@ -340,12 +340,12 @@ end $(Base.FastMath.make_fastmath(FASTABLE_AST)) end -# Not part of `FASTABLE_AST`: `hypot_fast` does not guard against overflow and underflow @testset "atan(y, x) for large and small arguments" begin setprecision(BigFloat, 256) do - @testset "atan($y, $x)" for (y, x) in - ((1e-200, 1e-200), (1e200, 3e200), (-1e-170, 3e-170)) - _, pb = rrule(atan, y, x) + @testset "$f($y, $x)" for f in (atan, Base.FastMath.atan_fast), + (y, x) in ((1e-200, 1e-200), (1e200, 3e200), (-1e-170, 3e-170)) + + _, pb = rrule(f, y, x) _, ȳ, x̄ = pb(1.0) u = big(x)^2 + big(y)^2 @test ȳ ≈ big(x) / u rtol = 4eps() From a2e18f7cac2c40b90e09c879a7732481e339ada0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Fri, 6 Nov 2026 14:27:06 +0100 Subject: [PATCH 7/7] Update src/rulesets/Base/fastmath_able.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/rulesets/Base/fastmath_able.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/rulesets/Base/fastmath_able.jl b/src/rulesets/Base/fastmath_able.jl index 03cc7324e..443943074 100644 --- a/src/rulesets/Base/fastmath_able.jl +++ b/src/rulesets/Base/fastmath_able.jl @@ -52,7 +52,7 @@ let # Trig-Multivariate # scaling by `max(|x|, |y|)` avoids overflow and underflow of `x ^ 2 + y ^ 2` @scalar_rule atan(y, x) @setup( - s = max(abs(x), abs(y)), xs = x / s, ys = y / s, v = xs ^ 2 + ys ^ 2 + s=max(abs(x), abs(y)), xs=x / s, ys=y / s, v=xs ^ 2 + ys ^ 2 ) (xs / v / s, -ys / v / s) @scalar_rule sincos(x) @setup((sinx, cosx)=Ω) cosx -sinx