Simplify quantile root-finding to Roots defaults

[devmotion]

First, thanks to @Cedriq1astaken for the work in #2066 — diagnosing the quantile_newton infinite loop, introducing Roots.find_zero, and assembling the regression cases. This PR builds directly on that branch (its commits are preserved here) and pushes the idea a little further.

I started looking into this more closely after hitting the same class of bug from a different direction in #2069 (adding cdf/ccdf/quantile to SkewNormal). There, the hand-rolled Newton iteration hangs whenever the root is near zero: SkewNormal's support straddles 0, and the relative-only stopping rule abs(x - x0) > max(abs(x), abs(x0)) * tol shrinks below the cdf/pdf evaluation noise as the root approaches 0, so the iterates oscillate in their last bits and never converge. That's the same underlying issue as the extreme-tail oscillation in #2061, just triggered by the root's location rather than by p. Seeing two distinct symptoms of one root cause is what motivated converging on a single, options-free formulation.

1) What changes compared to master

This adds a dependency on Roots v3 (compat Roots = "3").

master solves quantiles with hand-rolled Newton/bisection loops and hand-tuned tolerances. This replaces them with Roots.find_zero, relying on Roots' default tolerances plus a bracketing fallback — no tol/max_iters/floor parameters:

• quantile_bisect → find_zero(x -> cdf(d, x) - p, (lx, rx), ITP())
• quantile_newton / cquantile_newton / invlogcdf_newton / invlogccdf_newton → find_zero((F, f), x, Newton(), ITP()). Passing a bracketing method after Newton() makes Roots switch to ITP the moment the iteration brackets the root (a sign change of F).
• Chernoff quantiles now route through quantile_newton (with the same Normal-based warm start) instead of three bespoke newton(...) branches. Those branches called a newton function that was never defined anywhere in the package, so quantile(Chernoff(), τ) threw UndefVarError for any non-tabulated τ (e.g. truncated(Chernoff(); lower=0.1)) — present since the original Chernoff PR (Chernoff distribution #840). This removes the dead code and fixes Error with Truncated Chernoff quantiles; missing function (newton) #1999.

The reason defaults suffice is that Roots' tolerances supply three distinct ingredients, each covering one of the failure modes:

Failure mode	Issue(s)	Ingredient that fixes it
quantile_newton stalls/oscillates for a root near zero	#2069 (SkewNormal)	Newton's absolute xatol = eps (replaces #2069's hand-added eps^(3/4) floor)
quantile_newton oscillates in the extreme tail (p ≈ 1)	#1571, #1898, #2061	the ITP bracketing fallback (switch on sign change)
quantile_bisect can't terminate for a bracket far from zero	#1611, #1807	ITP's relative xrtol = eps (the old absolute cbrt(eps)^2 fell below the float spacing)

Roots' full defaults for reference:

	Newton()	ITP() / bracketing
xatol	eps(T)	eps(T)^3
xrtol	eps(T)	eps(T)
atol / rtol	4·eps(S)	0 / 0
maxiters / maxevals	40	60
strict	false	true

I confirmed the mapping empirically: for SkewNormal's near-zero roots, Newton() converges without the ITP fallback (the absolute xatol does the work); the fallback is what rescues the extreme-tail cases. Because Roots' defaults already cover the near-zero regime, #2069 no longer needs its own quantile_newton change.

Bracketing algorithm: standardized on ITP — NonlinearSolve's documented recommendation, explicitly suited to smooth, well-behaved functions. CDF roots are simple (d(cdf)/dx = pdf > 0), so it's a natural fit. (For reference: Roots' robust default is Bisection, its performance picks are A42/AlefeldPotraShi; the newer ModAB was recently added to both packages but isn't yet a documented recommendation, and its strength — roots of higher multiplicity — doesn't apply here.)

Fixes: #1571, #1611, #1807, #1898, #1999, #2061 (and removes the need for #2069's separate quantile_newton fix).

2) What changes compared to #2066

#2066 also moves to Roots.find_zero, but keeps hand-tuned machinery on top. Relative to that branch, this PR:

• Drops the tol and max_iters parameters and the max(tol, cbrt(eps)^2) floors from all four _newton functions — pure Roots defaults instead.
• Drops the xatol = tol override in quantile_bisect — pure Roots defaults there too.
• Drops the approximate "collapsed bracket" handling (isapprox-style rx ≈ lx → middle), keeping only the exact rx == lx degenerate guard (Sampling from truncated mixture model results in an error #1501). A non-bracketing interval (Bug in quantile() for some special cases of a MixtureModel #1869) is a separate problem from the solver choice and is intentionally left out of scope.
• Restructures the tests: one regression per @testset (no cases loops), and isapprox defaults instead of fixed atol. Adds coverage for [InverseGaussian?] quantile() hangs for larger number of σ's (quantile_newton not converging) #1898 and a Float32 round-trip.

🤖 Generated with Claude Code

[codecov-commenter]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 86.74%. Comparing base (03ddda7) to head (2561f2e).
⚠️ Report is 2 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2070      +/-   ##
==========================================
+ Coverage   86.63%   86.74%   +0.11%     
==========================================
  Files         149      149              
  Lines        8921     8882      -39     
==========================================
- Hits         7729     7705      -24     
+ Misses       1192     1177      -15     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[devmotion]

Thank you @Cedriq1astaken!