diff --git a/Project.toml b/Project.toml
index 092dfae81..72e8e1695 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Distributions"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
 authors = ["JuliaStats"]
-version = "0.25.127"
+version = "0.25.128"
 
 [deps]
 AliasTables = "66dad0bd-aa9a-41b7-9441-69ab47430ed8"
@@ -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"
@@ -51,6 +52,7 @@ PDMats = "0.11.35"
 Printf = "<0.0.1, 1"
 QuadGK = "2"
 Random = "<0.0.1, 1"
+Roots = "3"
 SparseArrays = "<0.0.1, 1"
 SparseConnectivityTracer = "1"
 SpecialFunctions = "1.2, 2"
diff --git a/src/Distributions.jl b/src/Distributions.jl
index ed2a2cdde..d00daf121 100644
--- a/src/Distributions.jl
+++ b/src/Distributions.jl
@@ -4,6 +4,7 @@ using StatsBase, PDMats, StatsFuns, Statistics
 using StatsFuns: logtwo, invsqrt2, invsqrt2π
 
 import QuadGK: quadgk
+using Roots: find_zero, ITP, Newton
 import Base: size, length, convert, show, getindex, rand, vec, inv
 import Base: sum, maximum, minimum, extrema, +, -, *, ==
 import Base.Math: @horner
diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index 8aa9f1b89..6ab209208 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -11,27 +11,11 @@ function quantile_bisect(d::ContinuousUnivariateDistribution, p::Real, lx::T, rx
         return middle(lx)
     end
 
-    # base tolerance on types to support e.g. `Float32` (avoids an infinite loop)
-    # ≈ 3.7e-11 for Float64
-    # ≈ 2.4e-5 for Float32
-    tol = cbrt(eps(float(T)))^2
-    # find quantile using bisect algorithm
-    cl = cdf(d, lx)
-    cr = cdf(d, rx)
-    cl <= p <= cr ||
-        throw(ArgumentError("[$lx, $rx] is not a valid bracketing interval for `quantile(d, $p)`"))
-    while rx - lx > tol
-        m = (lx + rx)/2
-        c = cdf(d, m)
-        if p > c
-            cl = c
-            lx = m
-        else
-            cr = c
-            rx = m
-        end
-    end
-    return (lx + rx)/2
+    # Rely on Roots' default tolerances for the bracketing solver. Unlike the previous
+    # hand-rolled bisection, ITP uses a relative `xrtol = eps`, so it converges even for
+    # brackets far from zero where an absolute tolerance is below the floating-point spacing
+    # (#1611).
+    return find_zero(x -> cdf(d, x) - p, (lx, rx), ITP())
 end
 
 function quantile_bisect(d::ContinuousUnivariateDistribution, p::Real, lx::Real, rx::Real)
@@ -47,16 +31,16 @@ 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)
+function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d))
     x = xs + (p - cdf(d, xs)) / pdf(d, xs)
     T = typeof(x)
+    F(x) = cdf(d, x) - p
+    f(x) = pdf(d, 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
+        # Newton's method with an ITP bracketing fallback: Roots switches to ITP whenever the
+        # iteration brackets the root (sign change of `F`), which resolves the oscillation/stalling
+        # observed for extreme quantiles (#2061, #1898). Roots' default tolerances are used.
+        return find_zero((F, f), x, Newton(), ITP())
     elseif p == 0
         return T(minimum(d))
     elseif p == 1
@@ -66,16 +50,13 @@ 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))
     x = xs + (ccdf(d, xs)-p) / pdf(d, xs)
     T = typeof(x)
+    F(x) = ccdf(d, x) - p
+    f(x) = -pdf(d, 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 find_zero((F, f), x, Newton(), ITP())
     elseif p == 1
         return T(minimum(d))
     elseif p == 0
@@ -85,24 +66,16 @@ 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))
     T = typeof(lp - logpdf(d,xs))
+    F(x) = logcdf(d, x) - lp
+    f(x) = exp(logpdf(d, x) - logcdf(d, x))
     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
-        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
-        end
-        return x
+        x = lp < logcdf(d,x0) ?
+            x0 - exp(lp - logpdf(d,x0) + logexpm1(max(logcdf(d,x0)-lp,0))) :
+            x0 + exp(lp - logpdf(d,x0) + log1mexp(min(logcdf(d,x0)-lp,0)))
+        return find_zero((F, f), x, Newton(), ITP())
     elseif lp == -Inf
         return T(minimum(d))
     elseif lp == 0
@@ -112,24 +85,16 @@ 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))
     T = typeof(lp - logpdf(d,xs))
+    F(x) = logccdf(d, x) - lp
+    f(x) = -exp(logpdf(d, x) - logccdf(d, x))
     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
-        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
-        end
-        return x
+        x = lp < logccdf(d,x0) ?
+            x0 + exp(lp - logpdf(d,x0) + logexpm1(max(logccdf(d,x0)-lp,0))) :
+            x0 - exp(lp - logpdf(d,x0) + log1mexp(min(logccdf(d,x0)-lp,0)))
+        return find_zero((F, f), x, Newton(), ITP())
     elseif lp == -Inf
         return T(maximum(d))
     elseif lp == 0
diff --git a/src/univariate/continuous/chernoff.jl b/src/univariate/continuous/chernoff.jl
index 3d319db25..a06dad7e8 100644
--- a/src/univariate/continuous/chernoff.jl
+++ b/src/univariate/continuous/chernoff.jl
@@ -189,14 +189,7 @@ function quantile(d::Chernoff, tau::Real)
     # one good approximation of the quantiles can be computed using Normal(0.0, stdapprox) with stdapprox = 0.52
     stdapprox = 0.52
     dnorm = Normal(0.0, 1.0)
-    if tau < 0.001
-        return -newton(x -> tau - ChernoffComputations._cdfbar(x), ChernoffComputations._pdf, quantile(dnorm, 1.0 - tau)*stdapprox)
-
-    end
-    if tau > 0.999
-        return newton(x -> 1.0 - tau - ChernoffComputations._cdfbar(x), ChernoffComputations._pdf, quantile(dnorm, tau)*stdapprox)
-    end
-    return newton(x -> ChernoffComputations._cdf(x) - tau, ChernoffComputations._pdf, quantile(dnorm, tau)*stdapprox)   # should consider replacing x-> construct for speed
+    return quantile_newton(d, tau, quantile(dnorm, tau)*stdapprox)
 end
 
 minimum(d::Chernoff) = -Inf
diff --git a/test/mixture.jl b/test/mixture.jl
index d2b2742eb..92493f1a0 100644
--- a/test/mixture.jl
+++ b/test/mixture.jl
@@ -290,4 +290,28 @@ end
             end
         end
     end
+
+    # issue #1611: bracket of adjacent floating point numbers far from zero
+    @testset "adjacent floating point bracket (#1611)" begin
+        d = MixtureModel(
+            [
+                Uniform{Float64}(-0.0001, 0.0001),
+                LogNormal{Float64}(11.174347445936371, 1.6086247197750911),
+            ],
+            [0.8832, 0.1168],
+        )
+        p = 0.99
+        x = quantile(d, p)
+        @test isfinite(x)
+        @test cdf(d, x) ≈ p
+    end
+
+    # issue #1807: high quantile of a mixture of exponentials with very different scales
+    @testset "high quantile of exponential mixture (#1807)" begin
+        d = MixtureModel(Exponential, [10_000, 1_000_000], [0.5, 0.5])
+        p = 0.999
+        x = quantile(d, p)
+        @test isfinite(x)
+        @test cdf(d, x) ≈ p
+    end
 end
diff --git a/test/quantile_newton.jl b/test/quantile_newton.jl
index da9291b0b..bb9bcc168 100644
--- a/test/quantile_newton.jl
+++ b/test/quantile_newton.jl
@@ -5,3 +5,43 @@ using Test
 d = Normal()
 @test Distributions.quantile_newton(d, 0.5) == quantile(d, 0.5)
 @test Distributions.cquantile_newton(d, 0.5) == cquantile(d, 0.5)
+
+# issue #1571
+@testset "InverseGaussian quantile convergence (#1571)" begin
+    d = InverseGaussian(1.0, 0.25)
+    p = 0.999996
+    x = quantile(d, p)
+    @test isfinite(x)
+    @test cdf(d, x) ≈ p
+end
+
+# issue #2061
+@testset "InverseGaussian quantile convergence (#2061)" begin
+    d = InverseGaussian(2.8853900817779268)
+    p = 0.9999996485182184
+    x = quantile(d, p)
+    @test isfinite(x)
+    @test cdf(d, x) ≈ p
+end
+
+# issue #1898: large-σ InverseGaussian, `p` far in the tail. `2 * cdf(Normal(), 5) - 1` is the
+# value `erf(5 / sqrt(2))` from the original report.
+@testset "InverseGaussian quantile convergence (#1898)" begin
+    d = InverseGaussian(1.187997687788096, 60.467382225458564)
+    p = 2 * cdf(Normal(), 5) - 1
+    x = quantile(d, p)
+    @test isfinite(x)
+    @test cdf(d, x) ≈ p
+end
+
+# the Roots defaults are type-aware: `Float32` must converge (and preserve the element type)
+# without the absolute tolerance floor that the hand-rolled loops needed to avoid an infinite loop.
+@testset "Float32 quantile round-trip" begin
+    d = InverseGaussian(2.0f0, 3.0f0)
+    for p in (0.1f0, 0.5f0, 0.9f0, 0.999f0)
+        x = quantile(d, p)
+        @test x isa Float32
+        @test isfinite(x)
+        @test cdf(d, x) ≈ p
+    end
+end
diff --git a/test/univariate/continuous/chernoff.jl b/test/univariate/continuous/chernoff.jl
index a5f28d02b..c72a3c049 100644
--- a/test/univariate/continuous/chernoff.jl
+++ b/test/univariate/continuous/chernoff.jl
@@ -155,4 +155,15 @@
     for i=1:size(pdftest, 1)
         @test isapprox(pdf(d, pdftest[i, 1]), pdftest[i, 2] ; atol = 1e-6)
     end
+
+    # issue #1999
+    @testset "truncated quantile for non-precomputed probability (#1999)" begin
+        td = truncated(Chernoff(); lower=0.1)
+        p = 0.75
+        x = quantile(td, p)
+
+        @test isfinite(x)
+        @test x >= 0.1
+        @test cdf(td, x) ≈ p
+    end
 end