Avoid overflow and underflow of y^2 in the derivatives of x / y

Project.toml

@@ -1,6 +1,6 @@
 name = "ForwardDiff"
 uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "1.4.6"
+version = "1.4.7"
 
 [deps]
 CommonSubexpressions = "bbf7d656-a473-5ed7-a52c-81e309532950"

src/partials.jl

@@ -96,8 +96,9 @@ end
     return Partials(mul_tuples(a.values, b.values, x_a, x_b))
 end
 
+# Divide by `bval` instead of multiplying by `inv(bval)` or dividing by `bval^2`, which may overflow or underflow
 @inline function _div_partials(a::Partials, b::Partials, aval, bval)
-    return _mul_partials(a, b, inv(bval), -(aval / (bval*bval)))
+    return (a - (aval / bval) * b) / bval
 end
 
 # NaN/Inf-safe methods #

test/DerivativeTest.jl

@@ -123,6 +123,19 @@ end
     end
 end
 
+@testset "division by tiny and huge numbers" begin
+    for x0 in (1e-200, 1e-310, 1e200)
+        @test iszero(ForwardDiff.derivative(x -> 2x / x, x0))
+        @test iszero(ForwardDiff.derivative(x -> zero(x) / x, x0))
+    end
+
+    # Third derivative of x / (x^2 + x^3) = 1 / (x + x^2) (#866)
+    d3(f, x) = ForwardDiff.derivative(x -> ForwardDiff.derivative(x -> ForwardDiff.derivative(f, x), x), x)
+    for x0 in (1e-20, 1e-30, 1e-40)
+        @test d3(x -> x / (x^2 + x^3), x0) ≈ d3(x -> 1 / (x + x^2), x0)
+    end
+end
+
 @testset "Givens rotations: Derivatives" begin
     # Test different branches in `LinearAlgebra.givensAlgorithm`
     for f in [randexp(), -randexp()], g in [0.0, f / 2, 2f, -f / 2, -2f], i in 1:3

test/PartialsTest.jl

@@ -132,7 +132,7 @@ samerng() = MersenneTwister(1)
         # The following properties should always be satisfied, regardless of whether NaN-safe mode is enabled or disabled
         # We use `isequal` for comparisons in the presence of `NaN`s
         for p1 in (PARTIALS, ALLZERO, FIRSTZERO), p2 in (PARTIALS2, ALLZERO, FIRSTZERO), v1 in (X, NaN, Inf), v2 in (Y, NaN, Inf)
-            @test isequal(ForwardDiff._div_partials(p1, p2, v1, v2), ForwardDiff._mul_partials(p1, p2, inv(v2), -v1/(v2^2)))
+            @test isequal(ForwardDiff._div_partials(p1, p2, v1, v2), (p1 - (v1 / v2) * p2) / v2)
             @test isequal(ForwardDiff._mul_partials(p1, p2, v1, v2), v1 * p1 + v2 * p2)
         end
         for v1 in (X, NaN, Inf), v2 in (Y, NaN, Inf)