diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 886724bd2..5d0505462 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -23,7 +23,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.3'
+          - '1.7'
           - '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.8' && '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 5cda14603..21ddd9a74 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 8aa9f1b89..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,16 +49,19 @@ quantile_bisect(d::ContinuousUnivariateDistribution, p::Real) =
 #   Distribution, with Application to the Inverse Gaussian Distribution
 #   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)
+function newton((f,df), xs::T=mode(d), xrtol::Real=1e-12) where {T}
+    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)
+    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)
+    x = xs - Δ(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)
-        end
-        return x
+        return newton((f, df), T(xs), xrtol)
     elseif p == 0
         return T(minimum(d))
     elseif p == 1
@@ -66,16 +71,15 @@ 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)
-    x = xs + (ccdf(d, xs)-p) / pdf(d, xs)
+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.
+    Δ(x) = f(x)/df(x)
+    x = xs - Δ(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)
-        end
-        return x
+        return newton((f, df), T(xs), xrtol)
     elseif p == 1
         return T(minimum(d))
     elseif p == 0
@@ -85,22 +89,17 @@ 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)))
+    df(x::T) where {T} = T(1)
     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)))
-            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
+            return newton((f_ver0,df), T(xs), xrtol)
         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
-                x0 = x
-                x = x0 + exp(lp - logpdf(d,x0) + log1mexp(min(logcdf(d,x0)-lp,0)))
-            end
+            return newton((f_ver1,df), T(xs), xrtol)
         end
         return x
     elseif lp == -Inf
@@ -112,22 +111,17 @@ 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)))
+    df(x::T) where {T} = T(1)
     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)))
-            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
+            return newton((f_ver0,df), T(xs), xrtol)
         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
-                x0 = x
-                x = x0 - exp(lp - logpdf(d,x0) + log1mexp(min(logccdf(d,x0)-lp,0)))
-            end
+            return newton((f_ver1,df), T(xs), xrtol)
         end
         return x
     elseif lp == -Inf
diff --git a/test/aqua.jl b/test/aqua.jl
index 4dd8c0706..17b06f424 100644
--- a/test/aqua.jl
+++ b/test/aqua.jl
@@ -9,13 +9,10 @@ import Aqua
     Aqua.test_all(
         Distributions;
         ambiguities = false,
-        # On older Julia versions, installed dependencies are quite old
-        # Thus unbound type parameters show up that are fixed in newer versions
-        unbound_args = VERSION >= v"1.6",
     )
     # Tests are not reliable on older Julia versions and
     # show ambiguities in loaded packages
-    if VERSION >= v"1.6"
+    if VERSION >= v"1.9"
         Aqua.test_ambiguities(Distributions)
     end
 end
diff --git a/test/testutils.jl b/test/testutils.jl
index 8acd37853..d1e215745 100644
--- a/test/testutils.jl
+++ b/test/testutils.jl
@@ -50,10 +50,8 @@ function test_cgf(dist, ts)
     d(f) = Base.Fix1(ForwardDiff.derivative, f)
     κ₁ = d(Base.Fix1(cgf, dist))(0)
     @test κ₁ ≈ mean(dist)
-    if VERSION >= v"1.4"
-        κ₂ = d(d(Base.Fix1(cgf, dist)))(0)
-        @test κ₂ ≈ var(dist)
-    end
+    κ₂ = d(d(Base.Fix1(cgf, dist)))(0)
+    @test κ₂ ≈ var(dist)
     for t in ts
         val = @inferred cgf(dist, t)
         @test isfinite(val)