Binomial distribution quantile edge case

[dschmitz89]

SciPy users once more found an edge case in one of the foundational distributions in scipy/scipy#22622. For extremely low probabilities, the quantile of the Binomial complementary CDF does not round trip. The reason I think this is fixable is that scipy's new root finder is able to find the right solution here so I guess a little tuning of the root finder could help here.

Reproducer:

#include <iostream>
#include <iomanip>
#include <boost/math/distributions/binomial.hpp>

int main() {
    using boost::math::binomial_distribution;
    using boost::math::complement;
    using boost::math::cdf;
    using boost::math::quantile;

    double n = 10000;
    double p = 4e-5;
    
    binomial_distribution<double> dist(n, p);

    // scipy: binom.sf(16, 10000, 4e-5) -> 3.268e-22
    double sf_value = cdf(complement(dist, 16.0));
    std::cout << "binom.sf(16, 10000, 4e-5) = " << std::scientific << std::setprecision(4) << sf_value << std::endl;

    // scipy: binom.isf(1e-20, 10000, 4e-5) -> 10000, expects <=16
    double q = 1e-20;
    double isf_value = quantile(complement(dist, q));
    std::cout << "binom.isf(1e-20, 10000, 4e-5) = " << std::fixed << std::setprecision(1) << isf_value << std::endl;
    # binom.sf(16, 10000, 4e-5) = 3.2685e-22
    #binom.isf(3.2685e-22, 10000, 4e-5) = 10000.0
    return 0;
}