Prevent passing denormals in calculation.

include/boost/math/concepts/distributions.hpp

@@ -343,7 +343,7 @@ struct DistributionConcept
    template <class R, class P>
    static void test_extra_members(const boost::math::hypergeometric_distribution<R, P>& d)
    {
-      unsigned u = d.defective();
+      std::uint64_t u = d.defective();
       u = d.sample_count();
       u = d.total();
       suppress_unused_variable_warning(u);

include/boost/math/concepts/std_real_concept.hpp

@@ -90,6 +90,9 @@ class std_real_concept
    std_real_concept& operator=(float c) { m_value = c; return *this; }
    std_real_concept& operator=(double c) { m_value = c; return *this; }
    std_real_concept& operator=(long double c) { m_value = c; return *this; }
+#ifdef BOOST_MATH_USE_FLOAT128
+   std_real_concept& operator=(BOOST_MATH_FLOAT128_TYPE c) { m_value = c; return *this; }
+#endif
 
    // Access:
    std_real_concept_base_type value()const{ return m_value; }

include/boost/math/distributions/detail/hypergeometric_pdf.hpp

@@ -297,7 +297,7 @@ T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hy
    while(data.current_prime <= data.N)
    {
       std::uint64_t base = data.current_prime;
-      std::int64_t prime_powers = 0;
+      std::uint64_t prime_powers = 0;
       while(base <= data.N)
       {
          prime_powers += data.n / base;
@@ -313,15 +313,15 @@ T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hy
       }
       if(prime_powers)
       {
-         T p = integer_power<T>(static_cast<T>(data.current_prime), prime_powers);
+         T p = integer_power<T>(static_cast<T>(data.current_prime), static_cast<int>(prime_powers));
          if((p > 1) && (tools::max_value<T>() / p < result.value))
          {
             //
             // The next calculation would overflow, use recursion
             // to sidestep the issue:
             //
             hypergeometric_pdf_prime_loop_result_entry<T> t = { p, &result };
-            data.current_prime = prime(++data.prime_index);
+            data.current_prime = prime(static_cast<unsigned>(++data.prime_index));
             return hypergeometric_pdf_prime_loop_imp<T>(data, t);
          }
          if((p < 1) && (tools::min_value<T>() / p > result.value))
@@ -331,12 +331,12 @@ T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hy
             // to sidestep the issue:
             //
             hypergeometric_pdf_prime_loop_result_entry<T> t = { p, &result };
-            data.current_prime = prime(++data.prime_index);
+            data.current_prime = prime(static_cast<unsigned>(++data.prime_index));
             return hypergeometric_pdf_prime_loop_imp<T>(data, t);
          }
          result.value *= p;
       }
-      data.current_prime = prime(++data.prime_index);
+      data.current_prime = prime(static_cast<unsigned>(++data.prime_index));
    }
    //
    // When we get to here we have run out of prime factors,
@@ -395,18 +395,18 @@ T hypergeometric_pdf_factorial_imp(std::uint64_t x, std::uint64_t r, std::uint64
 {
    BOOST_MATH_STD_USING
    BOOST_MATH_ASSERT(N <= boost::math::max_factorial<T>::value);
-   T result = boost::math::unchecked_factorial<T>(n);
+   T result = boost::math::unchecked_factorial<T>(static_cast<unsigned>(n));
    T num[3] = {
-      boost::math::unchecked_factorial<T>(r),
-      boost::math::unchecked_factorial<T>(N - n),
-      boost::math::unchecked_factorial<T>(N - r)
+      boost::math::unchecked_factorial<T>(static_cast<unsigned>(r)),
+      boost::math::unchecked_factorial<T>(static_cast<unsigned>(N - n)),
+      boost::math::unchecked_factorial<T>(static_cast<unsigned>(N - r))
    };
    T denom[5] = {
-      boost::math::unchecked_factorial<T>(N),
-      boost::math::unchecked_factorial<T>(x),
-      boost::math::unchecked_factorial<T>(n - x),
-      boost::math::unchecked_factorial<T>(r - x),
-      boost::math::unchecked_factorial<T>(N - n - r + x)
+      boost::math::unchecked_factorial<T>(static_cast<unsigned>(N)),
+      boost::math::unchecked_factorial<T>(static_cast<unsigned>(x)),
+      boost::math::unchecked_factorial<T>(static_cast<unsigned>(n - x)),
+      boost::math::unchecked_factorial<T>(static_cast<unsigned>(r - x)),
+      boost::math::unchecked_factorial<T>(static_cast<unsigned>(N - n - r + x))
    };
    std::size_t i = 0;
    std::size_t j = 0;

include/boost/math/distributions/detail/hypergeometric_quantile.hpp

@@ -14,7 +14,7 @@
 namespace boost{ namespace math{ namespace detail{
 
 template <class T>
-inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned lbound, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_down>&)
+inline std::uint64_t round_x_from_p(std::uint64_t x, T p, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t /*ubound*/, const policies::discrete_quantile<policies::integer_round_down>&)
 {
    if((p < cum * fudge_factor) && (x != lbound))
    {
@@ -25,7 +25,7 @@ inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned
 }
 
 template <class T>
-inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned /*lbound*/, unsigned ubound, const policies::discrete_quantile<policies::integer_round_up>&)
+inline std::uint64_t round_x_from_p(std::uint64_t x, T p, T cum, T fudge_factor, std::uint64_t /*lbound*/, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_up>&)
 {
    if((cum < p * fudge_factor) && (x != ubound))
    {
@@ -36,29 +36,29 @@ inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned
 }
 
 template <class T>
-inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_inwards>&)
+inline std::uint64_t round_x_from_p(std::uint64_t x, T p, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_inwards>&)
 {
    if(p >= 0.5)
       return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
    return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
 }
 
 template <class T>
-inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_outwards>&)
+inline std::uint64_t round_x_from_p(std::uint64_t x, T p, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_outwards>&)
 {
    if(p >= 0.5)
       return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
    return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
 }
 
 template <class T>
-inline unsigned round_x_from_p(unsigned x, T /*p*/, T /*cum*/, T /*fudge_factor*/, unsigned /*lbound*/, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&)
+inline std::uint64_t round_x_from_p(std::uint64_t x, T /*p*/, T /*cum*/, T /*fudge_factor*/, std::uint64_t /*lbound*/, std::uint64_t /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&)
 {
    return x;
 }
 
 template <class T>
-inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned lbound, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_down>&)
+inline std::uint64_t round_x_from_q(std::uint64_t x, T q, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t /*ubound*/, const policies::discrete_quantile<policies::integer_round_down>&)
 {
    if((q * fudge_factor > cum) && (x != lbound))
    {
@@ -69,7 +69,7 @@ inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned
 }
 
 template <class T>
-inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned /*lbound*/, unsigned ubound, const policies::discrete_quantile<policies::integer_round_up>&)
+inline std::uint64_t round_x_from_q(std::uint64_t x, T q, T cum, T fudge_factor, std::uint64_t /*lbound*/, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_up>&)
 {
    if((q < cum * fudge_factor) && (x != ubound))
    {
@@ -80,29 +80,29 @@ inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned
 }
 
 template <class T>
-inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_inwards>&)
+inline std::uint64_t round_x_from_q(std::uint64_t x, T q, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_inwards>&)
 {
    if(q < 0.5)
       return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
    return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
 }
 
 template <class T>
-inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_outwards>&)
+inline std::uint64_t round_x_from_q(std::uint64_t x, T q, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_outwards>&)
 {
    if(q >= 0.5)
       return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
    return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
 }
 
 template <class T>
-inline unsigned round_x_from_q(unsigned x, T /*q*/, T /*cum*/, T /*fudge_factor*/, unsigned /*lbound*/, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&)
+inline std::uint64_t round_x_from_q(std::uint64_t x, T /*q*/, T /*cum*/, T /*fudge_factor*/, std::uint64_t /*lbound*/, std::uint64_t /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&)
 {
    return x;
 }
 
 template <class T, class Policy>
-unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned N, const Policy& pol)
+std::uint64_t hypergeometric_quantile_imp(T p, T q, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Policy& pol)
 {
 #ifdef _MSC_VER
 #  pragma warning(push)
@@ -113,8 +113,8 @@ unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned
    BOOST_FPU_EXCEPTION_GUARD
    T result;
    T fudge_factor = 1 + tools::epsilon<T>() * ((N <= boost::math::prime(boost::math::max_prime - 1)) ? 50 : 2 * N);
-   unsigned base = static_cast<unsigned>((std::max)(0, static_cast<int>(n + r) - static_cast<int>(N)));
-   unsigned lim = (std::min)(r, n);
+   std::uint64_t base = static_cast<std::uint64_t>((std::max)(0, static_cast<int>(n + r) - static_cast<int>(N)));
+   std::uint64_t lim = (std::min)(r, n);
 
    BOOST_MATH_INSTRUMENT_VARIABLE(p);
    BOOST_MATH_INSTRUMENT_VARIABLE(q);
@@ -127,7 +127,7 @@ unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned
 
    if(p <= 0.5)
    {
-      unsigned x = base;
+      std::uint64_t x = base;
       result = hypergeometric_pdf<T>(x, r, n, N, pol);
       T diff = result;
       if (diff == 0)
@@ -175,7 +175,7 @@ unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned
    }
    else
    {
-      unsigned x = lim;
+      std::uint64_t x = lim;
       result = 0;
       T diff = hypergeometric_pdf<T>(x, r, n, N, pol);
       if (diff == 0)
@@ -225,7 +225,7 @@ unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned
 }
 
 template <class T, class Policy>
-inline unsigned hypergeometric_quantile(T p, T q, unsigned r, unsigned n, unsigned N, const Policy&)
+inline std::uint64_t hypergeometric_quantile(T p, T q, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Policy&)
 {
    BOOST_FPU_EXCEPTION_GUARD
    typedef typename tools::promote_args<T>::type result_type;

include/boost/math/distributions/non_central_beta.hpp

@@ -63,7 +63,7 @@ namespace boost
                ? detail::ibeta_imp(T(a + k), b, x, pol, false, true, &xterm)
                : detail::ibeta_imp(b, T(a + k), y, pol, true, true, &xterm);
 
-            while (beta * pois == 0)
+            while (fabs(beta * pois) < tools::min_value<T>())
             {
                if ((k == 0) || (pois == 0))
                   return init_val;
@@ -91,7 +91,7 @@ namespace boost
             {
                T term = beta * pois;
                sum += term;
-               if(((fabs(term/sum) < errtol) && (last_term >= term)) || (term == 0))
+               if(((fabs(term/sum) < errtol) && (fabs(last_term) >= fabs(term))) || (term == 0))
                {
                   count = k - i;
                   break;
@@ -106,6 +106,7 @@ namespace boost
 
                last_term = term;
             }
+            last_term = 0;
             for(auto i = k + 1; ; ++i)
             {
                poisf *= l2 / i;
@@ -114,10 +115,11 @@ namespace boost
 
                T term = poisf * betaf;
                sum += term;
-               if((fabs(term/sum) < errtol) || (term == 0))
+               if(((fabs(term/sum) < errtol) && (fabs(last_term) >= fabs(term))) || (term == 0))
                {
                   break;
                }
+               last_term = term;
                if(static_cast<std::uintmax_t>(count + i - k) > max_iter)
                {
                   return policies::raise_evaluation_error("cdf(non_central_beta_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE
@@ -554,7 +556,7 @@ namespace boost
                ibeta_derivative(a + k, b, x, pol)
                : ibeta_derivative(b, a + k, y, pol);
 
-            while (beta * pois == 0)
+            while (fabs(beta * pois) < tools::min_value<T>())
             {
                if ((k == 0) || (pois == 0))
                   return 0;  // Nothing else we can do!
@@ -579,33 +581,40 @@ namespace boost
             // Stable backwards recursion first:
             //
             std::uintmax_t count = k;
+            T ratio = 0;
+            T old_ratio = 0;
             for(auto i = k; i >= 0; --i)
             {
                T term = beta * pois;
                sum += term;
-               if((fabs(term/sum) < errtol) || (term == 0))
+               ratio = fabs(term / sum);
+               if(((ratio < errtol) && (ratio < old_ratio)) || (term == 0))
                {
                   count = k - i;
                   break;
                }
+               ratio = old_ratio;
                pois *= i / l2;
 
                if (a + b + i != 1)
                {
                   beta *= (a + i - 1) / (x * (a + i + b - 1));
                }
             }
+            old_ratio = 0;
             for(auto i = k + 1; ; ++i)
             {
                poisf *= l2 / i;
                betaf *= x * (a + b + i - 1) / (a + i - 1);
 
                T term = poisf * betaf;
                sum += term;
-               if((fabs(term/sum) < errtol) || (term == 0))
+               ratio = fabs(term / sum);
+               if(((ratio < errtol) && (ratio < old_ratio)) || (term == 0))
                {
                   break;
                }
+               old_ratio = ratio;
                if(static_cast<std::uintmax_t>(count + i - k) > max_iter)
                {
                   return policies::raise_evaluation_error("pdf(non_central_beta_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE