Failing gamma_p_inv() test case (when Briggs algos used)

[ckormanyos]

The following test (and other similar tests) fails on gamma_p_inv() for cpp_double_double when Briggs-style algorithms are used.

For cpp_double_double we get NaN, whereby the expected small-valued result should simply be zero for this type.

#include <iomanip>
#include <iostream>
#include <string>

#include <boost/math/special_functions/gamma.hpp>
#include <boost/multiprecision/cpp_double_fp.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>

#if defined(SC_)
#undef SC_
#endif

#define SC_(x) std::string( BOOST_STRINGIZE(x) )

int main()
{
  using boost::multiprecision::cpp_double_double;

  using float_type = typename cpp_double_double::backend_type::float_type;

  using dec_float_type = boost::multiprecision::number<boost::multiprecision::cpp_dec_float<31>, boost::multiprecision::et_off>;

  // Test case for gamma_p_inv.
  const std::string nums_as_str[] = { SC_(0.1730655412757187150418758392333984375e-5), SC_(0.12698681652545928955078125), SC_(0.2239623606222809074122747811596115646210220735131141509259977248899758059576948436798908594057794725e-517862), SC_(0.4348301951174619607003912855228982264838968134589390827069898370149065135278987288014463439625604227e-34079) };

  const auto a1 = cpp_double_double(nums_as_str[0U]);
  const auto p1 = cpp_double_double(nums_as_str[1U]);

  const auto a2 = dec_float_type(nums_as_str[0U]);
  const auto p2 = dec_float_type(nums_as_str[1U]);

  const auto gp1 = boost::math::gamma_p_inv(a1, p1);
  const auto gp2 = boost::math::gamma_p_inv(a2, p2);

  const auto ctrl1 = cpp_double_double(nums_as_str[2U]);
  const auto ctrl2 = dec_float_type   (nums_as_str[2U]);

  const auto delta1 = cpp_double_double(fabs(1 - fabs(gp1 / ctrl1)));
  const auto delta2 = dec_float_type   (fabs(1 - fabs(gp2 / ctrl2)));

  const auto flg = std::cout.flags();

  std::cout << std::setprecision(static_cast<std::streamsize>(std::numeric_limits<cpp_double_double>::digits10)) << gp1 << std::endl;
  std::cout << std::setprecision(static_cast<std::streamsize>(std::numeric_limits<dec_float_type   >::digits10)) << gp2 << std::endl;

  std::cout << std::setprecision(static_cast<std::streamsize>(std::numeric_limits<cpp_double_double>::digits10)) << ctrl1 << std::endl;
  std::cout << std::setprecision(static_cast<std::streamsize>(std::numeric_limits<dec_float_type   >::digits10)) << ctrl2 << std::endl;

  std::cout << std::setprecision(static_cast<std::streamsize>(std::numeric_limits<cpp_double_double>::digits10)) << delta1 << std::endl;
  std::cout << std::setprecision(static_cast<std::streamsize>(std::numeric_limits<dec_float_type   >::digits10)) << delta2 << std::endl;

  std::cout.flags(flg);
}