`InverseGamma::inverse_cdf` returns the wrong quantile at shape 151

[DRMacIver]

use statrs::distribution::{ContinuousCDF, InverseGamma};

fn main() {
    for shape in [140.0, 151.0] {
        let d = InverseGamma::new(shape, 4.5).unwrap();
        for p in [0.5, 0.9] {
            let q = d.inverse_cdf(p);
            println!("InverseGamma({shape}, 4.5).inverse_cdf({p}) = {q}, cdf of that = {}", d.cdf(q));
        }
    }
}

Output:

InverseGamma(140, 4.5).inverse_cdf(0.5) = 0.0322195377721101, cdf of that = 0.5000000000000567
InverseGamma(140, 4.5).inverse_cdf(0.9) = 0.03598080373442238, cdf of that = 0.8999999999999981
InverseGamma(151, 4.5).inverse_cdf(0.5) = 0.0234375, cdf of that = 0.0009628869544769414
InverseGamma(151, 4.5).inverse_cdf(0.9) = 0.0390625, cdf of that = 0.9991951730696086

cdf(inverse_cdf(p)) should be p, as it is at shape 140. At shape 151 the median comes back as the 0.001 quantile instead, and nothing reports a failure.

Possibly a duplicate of #453: pdf is inf across the mode at this shape, and fixing that would fix this. Filing separately because the fault is in newton_raphson_quantile, which accepts a Newton step of exactly zero as convergence, so it will keep returning wrong quantiles for whichever density overflows next. Close as a duplicate if you'd rather track it there.

Tested on statrs 0.19.1 and current main (50fcf4d).

(This was another bug found using hegel.)