From c443cbddc7eaad2a6b9846ea41f537b06b4503b7 Mon Sep 17 00:00:00 2001
From: Roman Lebedev <lebedev.ri@gmail.com>
Date: Thu, 12 Sep 2024 05:10:10 +0300
Subject: [PATCH 01/17] `quantile_newton()`: extract `df()`

---
 src/quantilealgs.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index 8aa9f1b89a..0cf41064ff 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -48,13 +48,14 @@ quantile_bisect(d::ContinuousUnivariateDistribution, p::Real) =
 #   http://www.statsci.org/smyth/pubs/qinvgaussPreprint.pdf
 
 function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12)
-    x = xs + (p - cdf(d, xs)) / pdf(d, xs)
+    f(x) = (p - cdf(d, x)) / pdf(d, x)
+    x = xs + f(xs)
     T = typeof(x)
     if 0 < p < 1
         x0 = T(xs)
         while abs(x-x0) > max(abs(x),abs(x0)) * tol
             x0 = x
-            x = x0 + (p - cdf(d, x0)) / pdf(d, x0)
+            x = x0 + f(x0)
         end
         return x
     elseif p == 0

From 4a40f790fe82beb5e35fa02edc02514bf338dff6 Mon Sep 17 00:00:00 2001
From: Roman Lebedev <lebedev.ri@gmail.com>
Date: Thu, 12 Sep 2024 05:37:44 +0300
Subject: [PATCH 02/17] `cquantile_newton()`: extract `df()`

---
 src/quantilealgs.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index 0cf41064ff..a1235df1fc 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -68,13 +68,14 @@ function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=
 end
 
 function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12)
-    x = xs + (ccdf(d, xs)-p) / pdf(d, xs)
+    f(x) = (ccdf(d, x)-p) / pdf(d, x)
+    x = xs + f(xs)
     T = typeof(x)
     if 0 < p < 1
         x0 = T(xs)
         while abs(x-x0) > max(abs(x),abs(x0)) * tol
             x0 = x
-            x = x0 + (ccdf(d, x0)-p) / pdf(d, x0)
+            x = x0 + f(x0)
         end
         return x
     elseif p == 1

From c02e6f7f5d7c0f3d45f940a74081bc7e592ef0a8 Mon Sep 17 00:00:00 2001
From: Roman Lebedev <lebedev.ri@gmail.com>
Date: Thu, 12 Sep 2024 05:38:16 +0300
Subject: [PATCH 03/17] `quantile_newton()`: extract `newton()`

---
 src/quantilealgs.jl | 19 +++++++++++++------
 1 file changed, 13 insertions(+), 6 deletions(-)

diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index a1235df1fc..79bdc4be7f 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -47,17 +47,24 @@ quantile_bisect(d::ContinuousUnivariateDistribution, p::Real) =
 #   Distribution, with Application to the Inverse Gaussian Distribution
 #   http://www.statsci.org/smyth/pubs/qinvgaussPreprint.pdf
 
+function newton(f, xs::T=mode(d), tol::Real=1e-12) where {T}
+    x = xs + f(xs)
+    @assert typeof(x) === T
+    x0 = T(xs)
+    while abs(x-x0) > max(abs(x),abs(x0)) * tol
+        x0 = x
+        x = x0 + f(x0)
+    end
+    return x
+end
+
 function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12)
     f(x) = (p - cdf(d, x)) / pdf(d, x)
+    # FIXME: can this be expressed via `promote_type()`? Test coverage missing.
     x = xs + f(xs)
     T = typeof(x)
     if 0 < p < 1
-        x0 = T(xs)
-        while abs(x-x0) > max(abs(x),abs(x0)) * tol
-            x0 = x
-            x = x0 + f(x0)
-        end
-        return x
+        return newton(f, T(xs), tol)
     elseif p == 0
         return T(minimum(d))
     elseif p == 1

From e0928ba3c14d8a8c0ec061ff43598a9906626b57 Mon Sep 17 00:00:00 2001
From: Roman Lebedev <lebedev.ri@gmail.com>
Date: Thu, 12 Sep 2024 05:38:33 +0300
Subject: [PATCH 04/17] `cquantile_newton()`: refactor using `newton()`

---
 src/quantilealgs.jl | 8 ++------
 1 file changed, 2 insertions(+), 6 deletions(-)

diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index 79bdc4be7f..60e5034c50 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -76,15 +76,11 @@ end
 
 function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12)
     f(x) = (ccdf(d, x)-p) / pdf(d, x)
+    # FIXME: can this be expressed via `promote_type()`? Test coverage missing.
     x = xs + f(xs)
     T = typeof(x)
     if 0 < p < 1
-        x0 = T(xs)
-        while abs(x-x0) > max(abs(x),abs(x0)) * tol
-            x0 = x
-            x = x0 + f(x0)
-        end
-        return x
+        return newton(f, T(xs), tol)
     elseif p == 1
         return T(minimum(d))
     elseif p == 0

From 8b1267e96950cb238c614b57bc93012a8fca8840 Mon Sep 17 00:00:00 2001
From: Roman Lebedev <lebedev.ri@gmail.com>
Date: Thu, 12 Sep 2024 06:20:51 +0300
Subject: [PATCH 05/17] `invlogcdf_newton()`/`invlogccdf_newton()`: use strict
 inequality

This is what the `quantile_newton()`/`cquantile_newton()` does,
because otherwise they were able to end up in an endless loops,
when the initial point and the mode are the same, see #666.

I'm not sure this is needed here, but the next change
is going to refactor them to use general `newton()`,
which would make this change anyway,
so unless we need the current behaviour,
let's do this change explicitly.
---
 src/quantilealgs.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index 60e5034c50..b1d1de3494 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -96,13 +96,13 @@ function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Rea
         x0 = T(xs)
         if lp < logcdf(d,x0)
             x = x0 - exp(lp - logpdf(d,x0) + logexpm1(max(logcdf(d,x0)-lp,0)))
-            while abs(x-x0) >= max(abs(x),abs(x0)) * tol
+            while abs(x-x0) > max(abs(x),abs(x0)) * tol
                 x0 = x
                 x = x0 - exp(lp - logpdf(d,x0) + logexpm1(max(logcdf(d,x0)-lp,0)))
             end
         else
             x = x0 + exp(lp - logpdf(d,x0) + log1mexp(min(logcdf(d,x0)-lp,0)))
-            while abs(x-x0) >= max(abs(x),abs(x0))*tol
+            while abs(x-x0) > max(abs(x),abs(x0))*tol
                 x0 = x
                 x = x0 + exp(lp - logpdf(d,x0) + log1mexp(min(logcdf(d,x0)-lp,0)))
             end
@@ -123,13 +123,13 @@ function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Re
         x0 = T(xs)
         if lp < logccdf(d,x0)
             x = x0 + exp(lp - logpdf(d,x0) + logexpm1(max(logccdf(d,x0)-lp,0)))
-            while abs(x-x0) >= max(abs(x),abs(x0)) * tol
+            while abs(x-x0) > max(abs(x),abs(x0)) * tol
                 x0 = x
                 x = x0 + exp(lp - logpdf(d,x0) + logexpm1(max(logccdf(d,x0)-lp,0)))
             end
         else
             x = x0 - exp(lp - logpdf(d,x0) + log1mexp(min(logccdf(d,x0)-lp,0)))
-            while abs(x-x0) >= max(abs(x),abs(x0)) * tol
+            while abs(x-x0) > max(abs(x),abs(x0)) * tol
                 x0 = x
                 x = x0 - exp(lp - logpdf(d,x0) + log1mexp(min(logccdf(d,x0)-lp,0)))
             end

From 905075a0bf43f9fa946a3ea213e003ab3cf09816 Mon Sep 17 00:00:00 2001
From: Roman Lebedev <lebedev.ri@gmail.com>
Date: Thu, 12 Sep 2024 06:49:46 +0300
Subject: [PATCH 06/17] `invlogcdf_newton()`: extract `df(x)`

---
 src/quantilealgs.jl | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index b1d1de3494..d46ba013e2 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -92,19 +92,21 @@ end
 
 function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12)
     T = typeof(lp - logpdf(d,xs))
+    f_a(x) = -exp(lp - logpdf(d,x) + logexpm1(max(logcdf(d,x)-lp,0)))
+    f_b(x) = exp(lp - logpdf(d,x) + log1mexp(min(logcdf(d,x)-lp,0)))
     if -Inf < lp < 0
         x0 = T(xs)
         if lp < logcdf(d,x0)
-            x = x0 - exp(lp - logpdf(d,x0) + logexpm1(max(logcdf(d,x0)-lp,0)))
+            x = x0 + f_a(x0)
             while abs(x-x0) > max(abs(x),abs(x0)) * tol
                 x0 = x
-                x = x0 - exp(lp - logpdf(d,x0) + logexpm1(max(logcdf(d,x0)-lp,0)))
+                x = x0 + f_a(x0)
             end
         else
-            x = x0 + exp(lp - logpdf(d,x0) + log1mexp(min(logcdf(d,x0)-lp,0)))
+            x = x0 + f_b(x0)
             while abs(x-x0) > max(abs(x),abs(x0))*tol
                 x0 = x
-                x = x0 + exp(lp - logpdf(d,x0) + log1mexp(min(logcdf(d,x0)-lp,0)))
+                x = x0 + f_b(x0)
             end
         end
         return x

From 8d0ebe2b79b331f1e2ecc21531aa42e0bf650c52 Mon Sep 17 00:00:00 2001
From: Roman Lebedev <lebedev.ri@gmail.com>
Date: Thu, 12 Sep 2024 06:56:31 +0300
Subject: [PATCH 07/17] `invlogccdf_newton()`: extract `df(x)`

---
 src/quantilealgs.jl | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index d46ba013e2..e7995b0746 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -121,19 +121,21 @@ end
 
 function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12)
     T = typeof(lp - logpdf(d,xs))
+    f_a(x) = exp(lp - logpdf(d,x) + logexpm1(max(logccdf(d,x)-lp,0)))
+    f_b(x) = -exp(lp - logpdf(d,x) + log1mexp(min(logccdf(d,x)-lp,0)))
     if -Inf < lp < 0
         x0 = T(xs)
         if lp < logccdf(d,x0)
-            x = x0 + exp(lp - logpdf(d,x0) + logexpm1(max(logccdf(d,x0)-lp,0)))
+            x = x0 + f_a(x0)
             while abs(x-x0) > max(abs(x),abs(x0)) * tol
                 x0 = x
-                x = x0 + exp(lp - logpdf(d,x0) + logexpm1(max(logccdf(d,x0)-lp,0)))
+                x = x0 + f_a(x0)
             end
         else
-            x = x0 - exp(lp - logpdf(d,x0) + log1mexp(min(logccdf(d,x0)-lp,0)))
+            x = x0 + f_b(x0)
             while abs(x-x0) > max(abs(x),abs(x0)) * tol
                 x0 = x
-                x = x0 - exp(lp - logpdf(d,x0) + log1mexp(min(logccdf(d,x0)-lp,0)))
+                x = x0 + f_b(x0)
             end
         end
         return x

From efe20ec75af2aa83cae94c47ce6c93eb8b32a0f8 Mon Sep 17 00:00:00 2001
From: Roman Lebedev <lebedev.ri@gmail.com>
Date: Thu, 12 Sep 2024 07:25:36 +0300
Subject: [PATCH 08/17] `invlogcdf_newton()`: refactor using `newton()`

---
 src/quantilealgs.jl | 12 ++----------
 1 file changed, 2 insertions(+), 10 deletions(-)

diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index e7995b0746..59c294f4f7 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -97,17 +97,9 @@ function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Rea
     if -Inf < lp < 0
         x0 = T(xs)
         if lp < logcdf(d,x0)
-            x = x0 + f_a(x0)
-            while abs(x-x0) > max(abs(x),abs(x0)) * tol
-                x0 = x
-                x = x0 + f_a(x0)
-            end
+            return newton(f_a, T(xs), tol)
         else
-            x = x0 + f_b(x0)
-            while abs(x-x0) > max(abs(x),abs(x0))*tol
-                x0 = x
-                x = x0 + f_b(x0)
-            end
+            return newton(f_b, T(xs), tol)
         end
         return x
     elseif lp == -Inf

From db6bf0aee348ab94c67ff3e032d300784f192406 Mon Sep 17 00:00:00 2001
From: Roman Lebedev <lebedev.ri@gmail.com>
Date: Thu, 12 Sep 2024 07:25:43 +0300
Subject: [PATCH 09/17] `invlogccdf_newton()`: refactor using `newton()`

---
 src/quantilealgs.jl | 12 ++----------
 1 file changed, 2 insertions(+), 10 deletions(-)

diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index 59c294f4f7..b987b39aec 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -118,17 +118,9 @@ function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Re
     if -Inf < lp < 0
         x0 = T(xs)
         if lp < logccdf(d,x0)
-            x = x0 + f_a(x0)
-            while abs(x-x0) > max(abs(x),abs(x0)) * tol
-                x0 = x
-                x = x0 + f_a(x0)
-            end
+            return newton(f_a, T(xs), tol)
         else
-            x = x0 + f_b(x0)
-            while abs(x-x0) > max(abs(x),abs(x0)) * tol
-                x0 = x
-                x = x0 + f_b(x0)
-            end
+            return newton(f_b, T(xs), tol)
         end
         return x
     elseif lp == -Inf

From b97bfcaa1fbf32421fa41d6044c9f559445b5439 Mon Sep 17 00:00:00 2001
From: Roman Lebedev <lebedev.ri@gmail.com>
Date: Fri, 13 Sep 2024 02:43:02 +0300
Subject: [PATCH 10/17] `newton()`: more closely match text-book method, which
 performs subtraction,

---
 src/quantilealgs.jl | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index b987b39aec..465a58281f 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -48,20 +48,20 @@ quantile_bisect(d::ContinuousUnivariateDistribution, p::Real) =
 #   http://www.statsci.org/smyth/pubs/qinvgaussPreprint.pdf
 
 function newton(f, xs::T=mode(d), tol::Real=1e-12) where {T}
-    x = xs + f(xs)
+    x = xs - f(xs)
     @assert typeof(x) === T
     x0 = T(xs)
     while abs(x-x0) > max(abs(x),abs(x0)) * tol
         x0 = x
-        x = x0 + f(x0)
+        x = x0 - f(x0)
     end
     return x
 end
 
 function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12)
-    f(x) = (p - cdf(d, x)) / pdf(d, x)
+    f(x) = -((p - cdf(d, x)) / pdf(d, x))
     # FIXME: can this be expressed via `promote_type()`? Test coverage missing.
-    x = xs + f(xs)
+    x = xs - f(xs)
     T = typeof(x)
     if 0 < p < 1
         return newton(f, T(xs), tol)
@@ -75,9 +75,9 @@ function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=
 end
 
 function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12)
-    f(x) = (ccdf(d, x)-p) / pdf(d, x)
+    f(x) = -((ccdf(d, x)-p) / pdf(d, x))
     # FIXME: can this be expressed via `promote_type()`? Test coverage missing.
-    x = xs + f(xs)
+    x = xs - f(xs)
     T = typeof(x)
     if 0 < p < 1
         return newton(f, T(xs), tol)
@@ -92,8 +92,8 @@ end
 
 function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12)
     T = typeof(lp - logpdf(d,xs))
-    f_a(x) = -exp(lp - logpdf(d,x) + logexpm1(max(logcdf(d,x)-lp,0)))
-    f_b(x) = exp(lp - logpdf(d,x) + log1mexp(min(logcdf(d,x)-lp,0)))
+    f_a(x) = exp(lp - logpdf(d,x) + logexpm1(max(logcdf(d,x)-lp,0)))
+    f_b(x) = -exp(lp - logpdf(d,x) + log1mexp(min(logcdf(d,x)-lp,0)))
     if -Inf < lp < 0
         x0 = T(xs)
         if lp < logcdf(d,x0)
@@ -113,8 +113,8 @@ end
 
 function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12)
     T = typeof(lp - logpdf(d,xs))
-    f_a(x) = exp(lp - logpdf(d,x) + logexpm1(max(logccdf(d,x)-lp,0)))
-    f_b(x) = -exp(lp - logpdf(d,x) + log1mexp(min(logccdf(d,x)-lp,0)))
+    f_a(x) = -exp(lp - logpdf(d,x) + logexpm1(max(logccdf(d,x)-lp,0)))
+    f_b(x) = exp(lp - logpdf(d,x) + log1mexp(min(logccdf(d,x)-lp,0)))
     if -Inf < lp < 0
         x0 = T(xs)
         if lp < logccdf(d,x0)

From fd6017dfd50d677157854798305fb28c63fdf26d Mon Sep 17 00:00:00 2001
From: Roman Lebedev <lebedev.ri@gmail.com>
Date: Fri, 13 Sep 2024 02:49:15 +0300
Subject: [PATCH 11/17] `newton()`: `f()` is more of a function to calculate
 the actual delta

---
 src/quantilealgs.jl | 34 +++++++++++++++++-----------------
 1 file changed, 17 insertions(+), 17 deletions(-)

diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index 465a58281f..e039a62471 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -47,24 +47,24 @@ quantile_bisect(d::ContinuousUnivariateDistribution, p::Real) =
 #   Distribution, with Application to the Inverse Gaussian Distribution
 #   http://www.statsci.org/smyth/pubs/qinvgaussPreprint.pdf
 
-function newton(f, xs::T=mode(d), tol::Real=1e-12) where {T}
-    x = xs - f(xs)
+function newton(Δ, xs::T=mode(d), tol::Real=1e-12) where {T}
+    x = xs - Δ(xs)
     @assert typeof(x) === T
     x0 = T(xs)
     while abs(x-x0) > max(abs(x),abs(x0)) * tol
         x0 = x
-        x = x0 - f(x0)
+        x = x0 - Δ(x0)
     end
     return x
 end
 
 function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12)
-    f(x) = -((p - cdf(d, x)) / pdf(d, x))
+    Δ(x) = -((p - cdf(d, x)) / pdf(d, x))
     # FIXME: can this be expressed via `promote_type()`? Test coverage missing.
-    x = xs - f(xs)
+    x = xs - Δ(xs)
     T = typeof(x)
     if 0 < p < 1
-        return newton(f, T(xs), tol)
+        return newton(Δ, T(xs), tol)
     elseif p == 0
         return T(minimum(d))
     elseif p == 1
@@ -75,12 +75,12 @@ function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=
 end
 
 function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12)
-    f(x) = -((ccdf(d, x)-p) / pdf(d, x))
+    Δ(x) = -((ccdf(d, x)-p) / pdf(d, x))
     # FIXME: can this be expressed via `promote_type()`? Test coverage missing.
-    x = xs - f(xs)
+    x = xs - Δ(xs)
     T = typeof(x)
     if 0 < p < 1
-        return newton(f, T(xs), tol)
+        return newton(Δ, T(xs), tol)
     elseif p == 1
         return T(minimum(d))
     elseif p == 0
@@ -92,14 +92,14 @@ end
 
 function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12)
     T = typeof(lp - logpdf(d,xs))
-    f_a(x) = exp(lp - logpdf(d,x) + logexpm1(max(logcdf(d,x)-lp,0)))
-    f_b(x) = -exp(lp - logpdf(d,x) + log1mexp(min(logcdf(d,x)-lp,0)))
+    Δ_ver0(x) = exp(lp - logpdf(d,x) + logexpm1(max(logcdf(d,x)-lp,0)))
+    Δ_ver1(x) = -exp(lp - logpdf(d,x) + log1mexp(min(logcdf(d,x)-lp,0)))
     if -Inf < lp < 0
         x0 = T(xs)
         if lp < logcdf(d,x0)
-            return newton(f_a, T(xs), tol)
+            return newton(Δ_ver0, T(xs), tol)
         else
-            return newton(f_b, T(xs), tol)
+            return newton(Δ_ver1, T(xs), tol)
         end
         return x
     elseif lp == -Inf
@@ -113,14 +113,14 @@ end
 
 function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12)
     T = typeof(lp - logpdf(d,xs))
-    f_a(x) = -exp(lp - logpdf(d,x) + logexpm1(max(logccdf(d,x)-lp,0)))
-    f_b(x) = exp(lp - logpdf(d,x) + log1mexp(min(logccdf(d,x)-lp,0)))
+    Δ_ver0(x) = -exp(lp - logpdf(d,x) + logexpm1(max(logccdf(d,x)-lp,0)))
+    Δ_ver1(x) = exp(lp - logpdf(d,x) + log1mexp(min(logccdf(d,x)-lp,0)))
     if -Inf < lp < 0
         x0 = T(xs)
         if lp < logccdf(d,x0)
-            return newton(f_a, T(xs), tol)
+            return newton(Δ_ver0, T(xs), tol)
         else
-            return newton(f_b, T(xs), tol)
+            return newton(Δ_ver1, T(xs), tol)
         end
         return x
     elseif lp == -Inf

From f3533d14eb4f14ae5addea95c6cde93686a2429d Mon Sep 17 00:00:00 2001
From: Roman Lebedev <lebedev.ri@gmail.com>
Date: Fri, 13 Sep 2024 02:52:51 +0300
Subject: [PATCH 12/17] `[c]quantile_newton()`: explicitly pass function and
 derivative into the `newton()`

---
 src/quantilealgs.jl | 27 ++++++++++++++++++---------
 1 file changed, 18 insertions(+), 9 deletions(-)

diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index e039a62471..362609356b 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -47,7 +47,7 @@ quantile_bisect(d::ContinuousUnivariateDistribution, p::Real) =
 #   Distribution, with Application to the Inverse Gaussian Distribution
 #   http://www.statsci.org/smyth/pubs/qinvgaussPreprint.pdf
 
-function newton(Δ, xs::T=mode(d), tol::Real=1e-12) where {T}
+function newton_impl(Δ, xs::T=mode(d), tol::Real=1e-12) where {T}
     x = xs - Δ(xs)
     @assert typeof(x) === T
     x0 = T(xs)
@@ -58,13 +58,20 @@ function newton(Δ, xs::T=mode(d), tol::Real=1e-12) where {T}
     return x
 end
 
+function newton((f,df), xs::T=mode(d), tol::Real=1e-12) where {T}
+    Δ(x) = f(x)/df(x)
+    return newton_impl(Δ, xs, tol)
+end
+
 function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12)
-    Δ(x) = -((p - cdf(d, x)) / pdf(d, x))
+    f(x) = -(p - cdf(d, x))
+    df(x) = pdf(d, x)
     # FIXME: can this be expressed via `promote_type()`? Test coverage missing.
+    Δ(x) = f(x)/df(x)
     x = xs - Δ(xs)
     T = typeof(x)
     if 0 < p < 1
-        return newton(Δ, T(xs), tol)
+        return newton((f, df), T(xs), tol)
     elseif p == 0
         return T(minimum(d))
     elseif p == 1
@@ -75,12 +82,14 @@ function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=
 end
 
 function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12)
-    Δ(x) = -((ccdf(d, x)-p) / pdf(d, x))
+    f(x) = -(ccdf(d, x)-p)
+    df(x) = pdf(d, x)
     # FIXME: can this be expressed via `promote_type()`? Test coverage missing.
+    Δ(x) = f(x)/df(x)
     x = xs - Δ(xs)
     T = typeof(x)
     if 0 < p < 1
-        return newton(Δ, T(xs), tol)
+        return newton((f, df), T(xs), tol)
     elseif p == 1
         return T(minimum(d))
     elseif p == 0
@@ -97,9 +106,9 @@ function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Rea
     if -Inf < lp < 0
         x0 = T(xs)
         if lp < logcdf(d,x0)
-            return newton(Δ_ver0, T(xs), tol)
+            return newton_impl(Δ_ver0, T(xs), tol)
         else
-            return newton(Δ_ver1, T(xs), tol)
+            return newton_impl(Δ_ver1, T(xs), tol)
         end
         return x
     elseif lp == -Inf
@@ -118,9 +127,9 @@ function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Re
     if -Inf < lp < 0
         x0 = T(xs)
         if lp < logccdf(d,x0)
-            return newton(Δ_ver0, T(xs), tol)
+            return newton_impl(Δ_ver0, T(xs), tol)
         else
-            return newton(Δ_ver1, T(xs), tol)
+            return newton_impl(Δ_ver1, T(xs), tol)
         end
         return x
     elseif lp == -Inf

From 46f4a106dd259e8fb7acbd662bc94efa52fe11b6 Mon Sep 17 00:00:00 2001
From: Roman Lebedev <lebedev.ri@gmail.com>
Date: Fri, 13 Sep 2024 18:42:55 +0300
Subject: [PATCH 13/17] `invlog[c]cdf_newton()`: explicitly pass function and
 derivative into the `newton()
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

These have `df(x)=1`, at least that is what the `df(x)` is if we solve
`Δ(x)=f(x)/f'(x)` as an equation for `f'(x)`.
---
 src/quantilealgs.jl | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index 362609356b..dc368298e6 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -101,14 +101,15 @@ end
 
 function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12)
     T = typeof(lp - logpdf(d,xs))
-    Δ_ver0(x) = exp(lp - logpdf(d,x) + logexpm1(max(logcdf(d,x)-lp,0)))
-    Δ_ver1(x) = -exp(lp - logpdf(d,x) + log1mexp(min(logcdf(d,x)-lp,0)))
+    f_ver0(x) = exp(lp - logpdf(d,x) + logexpm1(max(logcdf(d,x)-lp,0)))
+    f_ver1(x) = -exp(lp - logpdf(d,x) + log1mexp(min(logcdf(d,x)-lp,0)))
+    df(x::T) where {T} = T(1)
     if -Inf < lp < 0
         x0 = T(xs)
         if lp < logcdf(d,x0)
-            return newton_impl(Δ_ver0, T(xs), tol)
+            return newton((f_ver0,df), T(xs), tol)
         else
-            return newton_impl(Δ_ver1, T(xs), tol)
+            return newton((f_ver1,df), T(xs), tol)
         end
         return x
     elseif lp == -Inf
@@ -122,14 +123,15 @@ end
 
 function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12)
     T = typeof(lp - logpdf(d,xs))
-    Δ_ver0(x) = -exp(lp - logpdf(d,x) + logexpm1(max(logccdf(d,x)-lp,0)))
-    Δ_ver1(x) = exp(lp - logpdf(d,x) + log1mexp(min(logccdf(d,x)-lp,0)))
+    f_ver0(x) = -exp(lp - logpdf(d,x) + logexpm1(max(logccdf(d,x)-lp,0)))
+    f_ver1(x) = exp(lp - logpdf(d,x) + log1mexp(min(logccdf(d,x)-lp,0)))
+    df(x::T) where {T} = T(1)
     if -Inf < lp < 0
         x0 = T(xs)
         if lp < logccdf(d,x0)
-            return newton_impl(Δ_ver0, T(xs), tol)
+            return newton((f_ver0,df), T(xs), tol)
         else
-            return newton_impl(Δ_ver1, T(xs), tol)
+            return newton((f_ver1,df), T(xs), tol)
         end
         return x
     elseif lp == -Inf

From 4d021a701411bd058f9076c919190bf4ef3acfe2 Mon Sep 17 00:00:00 2001
From: Roman Lebedev <lebedev.ri@gmail.com>
Date: Fri, 13 Sep 2024 02:56:27 +0300
Subject: [PATCH 14/17] `[c]quantile_newton()`: sink negation into `f(x)`

---
 src/quantilealgs.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index dc368298e6..543a16f72a 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -64,7 +64,7 @@ function newton((f,df), xs::T=mode(d), tol::Real=1e-12) where {T}
 end
 
 function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12)
-    f(x) = -(p - cdf(d, x))
+    f(x) = cdf(d, x) - p
     df(x) = pdf(d, x)
     # FIXME: can this be expressed via `promote_type()`? Test coverage missing.
     Δ(x) = f(x)/df(x)
@@ -82,7 +82,7 @@ function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=
 end
 
 function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12)
-    f(x) = -(ccdf(d, x)-p)
+    f(x) = p - ccdf(d, x)
     df(x) = pdf(d, x)
     # FIXME: can this be expressed via `promote_type()`? Test coverage missing.
     Δ(x) = f(x)/df(x)

From 83992243fd8f18df2872b6bbe8497f22b76a1d51 Mon Sep 17 00:00:00 2001
From: Roman Lebedev <lebedev.ri@gmail.com>
Date: Fri, 13 Sep 2024 03:21:52 +0300
Subject: [PATCH 15/17] `newton_impl()`: use `isapprox()`

---
 src/quantilealgs.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index 543a16f72a..94608f5753 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -51,7 +51,7 @@ function newton_impl(Δ, xs::T=mode(d), tol::Real=1e-12) where {T}
     x = xs - Δ(xs)
     @assert typeof(x) === T
     x0 = T(xs)
-    while abs(x-x0) > max(abs(x),abs(x0)) * tol
+    while !isapprox(x, x0, atol=0, rtol=tol)
         x0 = x
         x = x0 - Δ(x0)
     end

From 4a66424939d6eec05afba853b4ec8169981ddeda Mon Sep 17 00:00:00 2001
From: Roman Lebedev <lebedev.ri@gmail.com>
Date: Fri, 13 Sep 2024 03:22:50 +0300
Subject: [PATCH 16/17] Be more specific that `tol` is really a `xrtol`

---
 src/quantilealgs.jl | 28 ++++++++++++++--------------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index 94608f5753..b6eb35cbd2 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -47,23 +47,23 @@ quantile_bisect(d::ContinuousUnivariateDistribution, p::Real) =
 #   Distribution, with Application to the Inverse Gaussian Distribution
 #   http://www.statsci.org/smyth/pubs/qinvgaussPreprint.pdf
 
-function newton_impl(Δ, xs::T=mode(d), tol::Real=1e-12) where {T}
+function newton_impl(Δ, xs::T=mode(d), xrtol::Real=1e-12) where {T}
     x = xs - Δ(xs)
     @assert typeof(x) === T
     x0 = T(xs)
-    while !isapprox(x, x0, atol=0, rtol=tol)
+    while !isapprox(x, x0, atol=0, rtol=xrtol)
         x0 = x
         x = x0 - Δ(x0)
     end
     return x
 end
 
-function newton((f,df), xs::T=mode(d), tol::Real=1e-12) where {T}
+function newton((f,df), xs::T=mode(d), xrtol::Real=1e-12) where {T}
     Δ(x) = f(x)/df(x)
-    return newton_impl(Δ, xs, tol)
+    return newton_impl(Δ, xs, xrtol)
 end
 
-function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12)
+function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), xrtol::Real=1e-12)
     f(x) = cdf(d, x) - p
     df(x) = pdf(d, x)
     # FIXME: can this be expressed via `promote_type()`? Test coverage missing.
@@ -71,7 +71,7 @@ function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=
     x = xs - Δ(xs)
     T = typeof(x)
     if 0 < p < 1
-        return newton((f, df), T(xs), tol)
+        return newton((f, df), T(xs), xrtol)
     elseif p == 0
         return T(minimum(d))
     elseif p == 1
@@ -81,7 +81,7 @@ function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=
     end
 end
 
-function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12)
+function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), xrtol::Real=1e-12)
     f(x) = p - ccdf(d, x)
     df(x) = pdf(d, x)
     # FIXME: can this be expressed via `promote_type()`? Test coverage missing.
@@ -89,7 +89,7 @@ function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real
     x = xs - Δ(xs)
     T = typeof(x)
     if 0 < p < 1
-        return newton((f, df), T(xs), tol)
+        return newton((f, df), T(xs), xrtol)
     elseif p == 1
         return T(minimum(d))
     elseif p == 0
@@ -99,7 +99,7 @@ function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real
     end
 end
 
-function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12)
+function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), xrtol::Real=1e-12)
     T = typeof(lp - logpdf(d,xs))
     f_ver0(x) = exp(lp - logpdf(d,x) + logexpm1(max(logcdf(d,x)-lp,0)))
     f_ver1(x) = -exp(lp - logpdf(d,x) + log1mexp(min(logcdf(d,x)-lp,0)))
@@ -107,9 +107,9 @@ function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Rea
     if -Inf < lp < 0
         x0 = T(xs)
         if lp < logcdf(d,x0)
-            return newton((f_ver0,df), T(xs), tol)
+            return newton((f_ver0,df), T(xs), xrtol)
         else
-            return newton((f_ver1,df), T(xs), tol)
+            return newton((f_ver1,df), T(xs), xrtol)
         end
         return x
     elseif lp == -Inf
@@ -121,7 +121,7 @@ function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Rea
     end
 end
 
-function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12)
+function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), xrtol::Real=1e-12)
     T = typeof(lp - logpdf(d,xs))
     f_ver0(x) = -exp(lp - logpdf(d,x) + logexpm1(max(logccdf(d,x)-lp,0)))
     f_ver1(x) = exp(lp - logpdf(d,x) + log1mexp(min(logccdf(d,x)-lp,0)))
@@ -129,9 +129,9 @@ function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Re
     if -Inf < lp < 0
         x0 = T(xs)
         if lp < logccdf(d,x0)
-            return newton((f_ver0,df), T(xs), tol)
+            return newton((f_ver0,df), T(xs), xrtol)
         else
-            return newton((f_ver1,df), T(xs), tol)
+            return newton((f_ver1,df), T(xs), xrtol)
         end
         return x
     elseif lp == -Inf

From 53c217192aa583118e58725118cbaac8321bd224 Mon Sep 17 00:00:00 2001
From: Roman Lebedev <lebedev.ri@gmail.com>
Date: Fri, 13 Sep 2024 02:32:55 +0300
Subject: [PATCH 17/17] `newton()`: defer to `Roots.jl` implementation

Note that we really do need `Roots.jl` 2.2.0-or-newer,
because of https://github.com/JuliaMath/Roots.jl/pull/445
---
 Project.toml        |  2 ++
 src/quantilealgs.jl | 16 +++-------------
 2 files changed, 5 insertions(+), 13 deletions(-)

diff --git a/Project.toml b/Project.toml
index b3ae9b837a..db0c7de3b9 100644
--- a/Project.toml
+++ b/Project.toml
@@ -13,6 +13,7 @@ PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
@@ -48,6 +49,7 @@ PDMats = "0.10, 0.11"
 Printf = "<0.0.1, 1"
 QuadGK = "2"
 Random = "<0.0.1, 1"
+Roots = "2.2"
 SparseArrays = "<0.0.1, 1"
 SpecialFunctions = "1.2, 2"
 StableRNGs = "1"
diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index b6eb35cbd2..4733e3dc0d 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -1,3 +1,5 @@
+using Roots
+
 # Various algorithms for computing quantile
 
 function quantile_bisect(d::ContinuousUnivariateDistribution, p::Real, lx::T, rx::T) where {T<:Real}
@@ -47,20 +49,8 @@ quantile_bisect(d::ContinuousUnivariateDistribution, p::Real) =
 #   Distribution, with Application to the Inverse Gaussian Distribution
 #   http://www.statsci.org/smyth/pubs/qinvgaussPreprint.pdf
 
-function newton_impl(Δ, xs::T=mode(d), xrtol::Real=1e-12) where {T}
-    x = xs - Δ(xs)
-    @assert typeof(x) === T
-    x0 = T(xs)
-    while !isapprox(x, x0, atol=0, rtol=xrtol)
-        x0 = x
-        x = x0 - Δ(x0)
-    end
-    return x
-end
-
 function newton((f,df), xs::T=mode(d), xrtol::Real=1e-12) where {T}
-    Δ(x) = f(x)/df(x)
-    return newton_impl(Δ, xs, xrtol)
+    return find_zero((f,df), xs, Roots.Newton(), xatol=0, xrtol=xrtol, atol=0, rtol=eps(float(T)), maxiters=typemax(Int))
 end
 
 function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), xrtol::Real=1e-12)