From 5a13d0280359392d0671416452df315a217c1a5a 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] `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

This forces minimal Julia version to be 1.6
---
 .github/workflows/CI.yml |  4 ++--
 Project.toml             |  4 +++-
 src/quantilealgs.jl      | 16 +++-------------
 3 files changed, 8 insertions(+), 16 deletions(-)

diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 886724bd2..90f808845 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -23,7 +23,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.3'
+          - '1.6'
           - '1'
           - pre
         os:
@@ -36,7 +36,7 @@ jobs:
         with:
           version: ${{ matrix.version }}
           # ARM64 on macos-latest is neither supported by older Julia versions nor setup-julia
-          arch: ${{ matrix.os == 'macos-latest' && matrix.version != '1.3' && 'aarch64' || 'x64' }}
+          arch: ${{ matrix.os == 'macos-latest' && matrix.version <= '1.7' && 'aarch64' || 'x64' }}
           show-versioninfo: true
       - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
diff --git a/Project.toml b/Project.toml
index af12e4b7a..aab9155d1 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"
@@ -57,7 +59,7 @@ StatsAPI = "1.6"
 StatsBase = "0.32, 0.33, 0.34"
 StatsFuns = "0.9.15, 1"
 Test = "<0.0.1, 1"
-julia = "1.3"
+julia = "1.6"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index b6eb35cbd..4733e3dc0 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)