beta_distribution quantile returning incorrect values when precision policy is used.

[OliverHallam]

We've run into an issue where the beta quantile function returns incorrect results for particular inputs, when using a precision policy.

Here is a simple repro of one particular case we hit:

#include <iostream>
#include <boost/math/distributions/beta.hpp>

using namespace boost::math;
using namespace boost::math::policies;

int main()
{
    double alpha = 1.78672964680176;
    double beta = 2583.71522530538;
    double p = 0.0273644014531048;

    beta_distribution<double, policy<digits10<8>>> dist{alpha, beta};

    auto x = quantile(dist, p);

    std::cout << x << "\n";
}

The expected result is 7.33064e-05, but the result returned is 0.253938.

If you change the alpha/beta values slightly the correct value is returned. Removing the digits policy also seems to eliminate the error.

In my brief investigation it looks like the initial approximation is correct, but it then jumps to the wrong value in the first iteration of halley_iterate. My uneducated guess is that it looks like its using the same precision to evaluate the beta function when its rootfinding to find the inverse - intuitively I'm not sure that's valid, although this judgement is well beyond my current mathematical prowess.

FWIW we were observing this bug on boost 1.69 but is still reproducible in 1.91 so its not something that has changed recently. I have reproduced this in Clang, but we were originally seeing this on production code built with ICC so I don't think its a compiler issue.

[OliverHallam]

This should have been reported against boostorg/math - can someone move it for me?