bug: `inv_beta_reg` panics, does not return, and is non-monotone for small `x`

[FBruzzesi]

Same context and premise as in #432

---

Three separate issues in the lower tail, so please treat them independently. All three are also reachable from Beta::inverse_cdf.

#	Issue	Trigger
1	panic	x below about 1.6e-162
2	"never" returns	x below about 1e-50 at (a, b) = (200, 2)
3	constant and non-monotone results	any root below about 1e-17

Environment: statrs 0.19.0 (default-features = false, features = ["std"]), rustc 1.97.1, aarch64-apple-darwin, --release, so debug_assert! is off.

Issue 1: panic for very small x

use statrs::function::beta::inv_beta_reg;

fn main() {
    println!("{}", inv_beta_reg(200.0, 2.0, 1e-165)); // panics immediately
}

thread 'main' panicked at .../statrs-0.19.0/src/function/beta.rs:132:31:
called `Result::unwrap()` on an `Err` value: XOutOfRange

Also reachable from the public API:

use statrs::distribution::{Beta, ContinuousCDF};
Beta::new(200.0, 2.0).unwrap().inverse_cdf(1e-170); // panics

Issue 2: no termination for a band of x

use statrs::function::beta::inv_beta_reg;

fn main() {
    println!("{}", inv_beta_reg(200.0, 2.0, 1e-40)); // returns in microseconds
    println!("{}", inv_beta_reg(200.0, 2.0, 1e-60)); // does not return
}

(I killed the second call after ~5 mins)

Issue 3: the result stops moving, and is not monotone

use statrs::function::beta::inv_beta_reg;

fn main() {
    for q in [1e-2, 1e-4, 1e-10, 1e-20, 1e-30] {
        println!("inv_beta_reg(0.1, 500, {q:e}) = {:e}", inv_beta_reg(0.1, 500.0, q));
    }
}

Reference from scipy.stats.beta(0.1, 500).ppf(q), cross-checked with mpmath:

q	statrs	true
1e-2	2.5310279925347366e-18	1.2157036049543745e-23
1e-4	5.297481902267595e-17	1.2157036049543762e-43
1e-10	5.164813474103027e-17	1.215703604954385e-103
1e-20	5.164813342904473e-17	1.2157036049543978e-203
1e-30	5.164813342904473e-17	1.2157036049544177e-303

From q = 1e-10 down, every call returns 5.1648133429e-17. The results are also not monotone in q: 1e-2 gives 2.53e-18, which is smaller than the 5.30e-17 returned for 1e-4, while the quantile function is increasing. Anything that brackets or bisects on the inverse CDF relies on that ordering.

The same floor shows up on better-behaved shapes, less dramatically:

q	inv_beta_reg(2, 200, q)	true	rel. err.
1e-20	7.053459200759641e-13	7.053456158916007e-13	4.3e-7
1e-30	7.180186771679575e-16	7.053456158585999e-18	1e2
1e-40	7.179955800324264e-16	7.053456158585986e-23	1e7

---

FYI I would be interested in trying opening a PR for one or more of these issues