Merge pull request #987 from mborland/18511

Fix for scipy issues 18506 and 18511

Author: mborland

doc/distributions/hypergeometric.qbk

@@ -15,11 +15,11 @@
       typedef RealType value_type;
       typedef Policy   policy_type;
       // Construct:
-      hypergeometric_distribution(unsigned r, unsigned n, unsigned N); // r=defective/failures/success, n=trials/draws, N=total population.
+      hypergeometric_distribution(uint64_t r, uint64_t n, uint64_t N); // r=defective/failures/success, n=trials/draws, N=total population.
       // Accessors:
-      unsigned total()const;
-      unsigned defective()const;
-      unsigned sample_count()const;
+      uint64_t total()const;
+      uint64_t defective()const;
+      uint64_t sample_count()const;
    };
 
    typedef hypergeometric_distribution<> hypergeometric;
@@ -56,20 +56,20 @@ then we obtain basically the same graphs:
 
 [h4 Member Functions]
 
-   hypergeometric_distribution(unsigned r, unsigned n, unsigned N);
+   hypergeometric_distribution(uint64_t r, uint64_t n, uint64_t N);
 
 Constructs a hypergeometric distribution with a population of /N/ objects,
 of which /r/ are defective, and from which /n/ are sampled.
 
-   unsigned total()const;
+   uint64_t total()const;
 
 Returns the total number of objects /N/.
 
-   unsigned defective()const;
+   uint64_t defective()const;
 
 Returns the number of objects /r/ in population /N/ which are defective.
 
-   unsigned sample_count()const;
+   uint64_t sample_count()const;
 
 Returns the number of objects /n/ which are sampled from the population /N/.
 
@@ -87,7 +87,7 @@ and Python
 All the [link math_toolkit.dist_ref.nmp usual non-member accessor functions]
 that are generic to all distributions are supported: __usual_accessors.
 
-The domain of the random variable is the unsigned integers in the range
+The domain of the random variable are the 64-bit unsigned integers in the range
 \[max(0, n + r - N), min(n, r)\].  A __domain_error is raised if the
 random variable is outside this range, or is not an integral value.
 

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

@@ -10,11 +10,12 @@
 
 #include <boost/math/policies/error_handling.hpp>
 #include <boost/math/distributions/detail/hypergeometric_pdf.hpp>
+#include <cstdint>
 
 namespace boost{ namespace math{ namespace detail{
 
    template <class T, class Policy>
-   T hypergeometric_cdf_imp(unsigned x, unsigned r, unsigned n, unsigned N, bool invert, const Policy& pol)
+   T hypergeometric_cdf_imp(std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, bool invert, const Policy& pol)
    {
 #ifdef _MSC_VER
 #  pragma warning(push)
@@ -27,7 +28,7 @@ namespace boost{ namespace math{ namespace detail{
       {
          result = hypergeometric_pdf<T>(x, r, n, N, pol);
          T diff = result;
-         unsigned lower_limit = static_cast<unsigned>((std::max)(0, static_cast<int>(n + r) - static_cast<int>(N)));
+         const auto lower_limit = static_cast<std::uint64_t>((std::max)(INT64_C(0), static_cast<std::int64_t>(n + r) - static_cast<std::int64_t>(N)));
          while(diff > (invert ? T(1) : result) * tools::epsilon<T>())
          {
             diff = T(x) * T((N + x) - n - r) * diff / (T(1 + n - x) * T(1 + r - x));
@@ -43,7 +44,7 @@ namespace boost{ namespace math{ namespace detail{
       else
       {
          invert = !invert;
-         unsigned upper_limit = (std::min)(r, n);
+         const auto upper_limit = (std::min)(r, n);
          if(x != upper_limit)
          {
             ++x;
@@ -69,7 +70,7 @@ namespace boost{ namespace math{ namespace detail{
    }
 
    template <class T, class Policy>
-   inline T hypergeometric_cdf(unsigned x, unsigned r, unsigned n, unsigned N, bool invert, const Policy&)
+   inline T hypergeometric_cdf(std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, bool invert, const Policy&)
    {
       BOOST_FPU_EXCEPTION_GUARD
       typedef typename tools::promote_args<T>::type result_type;

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

@@ -15,6 +15,8 @@
 #include <boost/math/special_functions/pow.hpp>
 #include <boost/math/special_functions/prime.hpp>
 #include <boost/math/policies/error_handling.hpp>
+#include <algorithm>
+#include <cstdint>
 
 #ifdef BOOST_MATH_INSTRUMENT
 #include <typeinfo>
@@ -40,7 +42,7 @@ template <class T>
 struct sort_functor
 {
    sort_functor(const T* exponents) : m_exponents(exponents){}
-   bool operator()(int i, int j)
+   bool operator()(std::size_t i, std::size_t j)
    {
       return m_exponents[i] > m_exponents[j];
    }
@@ -49,7 +51,7 @@ struct sort_functor
 };
 
 template <class T, class Lanczos, class Policy>
-T hypergeometric_pdf_lanczos_imp(T /*dummy*/, unsigned x, unsigned r, unsigned n, unsigned N, const Lanczos&, const Policy&)
+T hypergeometric_pdf_lanczos_imp(T /*dummy*/, std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Lanczos&, const Policy&)
 {
    BOOST_MATH_STD_USING
 
@@ -137,7 +139,7 @@ T hypergeometric_pdf_lanczos_imp(T /*dummy*/, unsigned x, unsigned r, unsigned n
    //
    // Combine equal powers:
    //
-   int j = 8;
+   std::size_t j = 8;
    while(exponents[sorted_indexes[j]] == 0) --j;
    while(j)
    {
@@ -184,7 +186,7 @@ T hypergeometric_pdf_lanczos_imp(T /*dummy*/, unsigned x, unsigned r, unsigned n
       result = pow(bases[sorted_indexes[0]] * exp(static_cast<T>(base_e_factors[sorted_indexes[0]])), exponents[sorted_indexes[0]]);
    }
    BOOST_MATH_INSTRUMENT_VARIABLE(result);
-   for(unsigned i = 1; (i < 9) && (exponents[sorted_indexes[i]] > 0); ++i)
+   for(std::size_t i = 1; (i < 9) && (exponents[sorted_indexes[i]] > 0); ++i)
    {
       BOOST_FPU_EXCEPTION_GUARD
       if(result < tools::min_value<T>())
@@ -215,7 +217,7 @@ T hypergeometric_pdf_lanczos_imp(T /*dummy*/, unsigned x, unsigned r, unsigned n
 }
 
 template <class T, class Policy>
-T hypergeometric_pdf_lanczos_imp(T /*dummy*/, unsigned x, unsigned r, unsigned n, unsigned N, const boost::math::lanczos::undefined_lanczos&, const Policy& pol)
+T hypergeometric_pdf_lanczos_imp(T /*dummy*/, std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, const boost::math::lanczos::undefined_lanczos&, const Policy& pol)
 {
    BOOST_MATH_STD_USING
    return exp(
@@ -277,12 +279,12 @@ struct hypergeometric_pdf_prime_loop_result_entry
 
 struct hypergeometric_pdf_prime_loop_data
 {
-   const unsigned x;
-   const unsigned r;
-   const unsigned n;
-   const unsigned N;
-   unsigned prime_index;
-   unsigned current_prime;
+   const std::uint64_t x;
+   const std::uint64_t r;
+   const std::uint64_t n;
+   const std::uint64_t N;
+   std::size_t prime_index;
+   std::uint64_t current_prime;
 };
 
 #ifdef _MSC_VER
@@ -294,8 +296,8 @@ T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hy
 {
    while(data.current_prime <= data.N)
    {
-      unsigned base = data.current_prime;
-      int prime_powers = 0;
+      std::uint64_t base = data.current_prime;
+      std::int64_t prime_powers = 0;
       while(base <= data.N)
       {
          prime_powers += data.n / base;
@@ -381,15 +383,15 @@ T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hy
 }
 
 template <class T, class Policy>
-inline T hypergeometric_pdf_prime_imp(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&)
+inline T hypergeometric_pdf_prime_imp(std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Policy&)
 {
    hypergeometric_pdf_prime_loop_result_entry<T> result = { 1, 0 };
    hypergeometric_pdf_prime_loop_data data = { x, r, n, N, 0, prime(0) };
    return hypergeometric_pdf_prime_loop_imp<T>(data, result);
 }
 
 template <class T, class Policy>
-T hypergeometric_pdf_factorial_imp(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&)
+T hypergeometric_pdf_factorial_imp(std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Policy&)
 {
    BOOST_MATH_STD_USING
    BOOST_MATH_ASSERT(N <= boost::math::max_factorial<T>::value);
@@ -406,8 +408,8 @@ T hypergeometric_pdf_factorial_imp(unsigned x, unsigned r, unsigned n, unsigned
       boost::math::unchecked_factorial<T>(r - x),
       boost::math::unchecked_factorial<T>(N - n - r + x)
    };
-   int i = 0;
-   int j = 0;
+   std::size_t i = 0;
+   std::size_t j = 0;
    while((i < 3) || (j < 5))
    {
       while((j < 5) && ((result >= 1) || (i >= 3)))
@@ -427,7 +429,7 @@ T hypergeometric_pdf_factorial_imp(unsigned x, unsigned r, unsigned n, unsigned
 
 template <class T, class Policy>
 inline typename tools::promote_args<T>::type 
-   hypergeometric_pdf(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&)
+   hypergeometric_pdf(std::uint64_t x, 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/hypergeometric.hpp

@@ -16,6 +16,7 @@
 #include <boost/math/distributions/detail/hypergeometric_cdf.hpp>
 #include <boost/math/distributions/detail/hypergeometric_quantile.hpp>
 #include <boost/math/special_functions/fpclassify.hpp>
+#include <cstdint>
 
 namespace boost { namespace math {
 
@@ -26,25 +27,25 @@ namespace boost { namespace math {
       typedef RealType value_type;
       typedef Policy policy_type;
 
-      hypergeometric_distribution(unsigned r, unsigned n, unsigned N) // Constructor. r=defective/failures/success, n=trials/draws, N=total population.
+      hypergeometric_distribution(std::uint64_t r, std::uint64_t n, std::uint64_t N) // Constructor. r=defective/failures/success, n=trials/draws, N=total population.
          : m_n(n), m_N(N), m_r(r)
       {
          static const char* function = "boost::math::hypergeometric_distribution<%1%>::hypergeometric_distribution";
          RealType ret;
          check_params(function, &ret);
       }
       // Accessor functions.
-      unsigned total()const
+      std::uint64_t total() const
       {
          return m_N;
       }
 
-      unsigned defective()const // successes/failures/events
+      std::uint64_t defective() const // successes/failures/events
       {
          return m_r;
       }
 
-      unsigned sample_count()const
+      std::uint64_t sample_count()const
       {
          return m_n;
       }
@@ -65,9 +66,9 @@ namespace boost { namespace math {
          }
          return true;
       }
-      bool check_x(unsigned x, const char* function, RealType* result)const
+      bool check_x(std::uint64_t x, const char* function, RealType* result)const
       {
-         if(x < static_cast<unsigned>((std::max)(0, static_cast<int>(m_n + m_r) - static_cast<int>(m_N))))
+         if(x < static_cast<std::uint64_t>((std::max)(INT64_C(0), static_cast<std::int64_t>(m_n + m_r) - static_cast<std::int64_t>(m_N))))
          {
             *result = boost::math::policies::raise_domain_error<RealType>(
                function, "Random variable out of range: must be > 0 and > m + r - N but got %1%", static_cast<RealType>(x), Policy());
@@ -84,40 +85,40 @@ namespace boost { namespace math {
 
    private:
       // Data members:
-      unsigned m_n;  // number of items picked or drawn.
-      unsigned m_N; // number of "total" items.
-      unsigned m_r; // number of "defective/successes/failures/events items.
+      std::uint64_t m_n;  // number of items picked or drawn.
+      std::uint64_t m_N; // number of "total" items.
+      std::uint64_t m_r; // number of "defective/successes/failures/events items.
 
    }; // class hypergeometric_distribution
 
    typedef hypergeometric_distribution<double> hypergeometric;
 
    template <class RealType, class Policy>
-   inline const std::pair<unsigned, unsigned> range(const hypergeometric_distribution<RealType, Policy>& dist)
+   inline const std::pair<std::uint64_t, std::uint64_t> range(const hypergeometric_distribution<RealType, Policy>& dist)
    { // Range of permissible values for random variable x.
 #ifdef _MSC_VER
 #  pragma warning(push)
 #  pragma warning(disable:4267)
 #endif
-      unsigned r = dist.defective();
-      unsigned n = dist.sample_count();
-      unsigned N = dist.total();
-      unsigned l = static_cast<unsigned>((std::max)(0, static_cast<int>(n + r) - static_cast<int>(N)));
-      unsigned u = (std::min)(r, n);
-      return std::pair<unsigned, unsigned>(l, u);
+      const auto r = dist.defective();
+      const auto n = dist.sample_count();
+      const auto N = dist.total();
+      const auto l = static_cast<std::uint64_t>((std::max)(INT64_C(0), static_cast<std::int64_t>(n + r) - static_cast<std::int64_t>(N)));
+      const auto u = (std::min)(r, n);
+      return std::make_pair(l, u);
 #ifdef _MSC_VER
 #  pragma warning(pop)
 #endif
    }
 
    template <class RealType, class Policy>
-   inline const std::pair<unsigned, unsigned> support(const hypergeometric_distribution<RealType, Policy>& d)
+   inline const std::pair<std::uint64_t, std::uint64_t> support(const hypergeometric_distribution<RealType, Policy>& d)
    {
       return range(d);
    }
 
    template <class RealType, class Policy>
-   inline RealType pdf(const hypergeometric_distribution<RealType, Policy>& dist, const unsigned& x)
+   inline RealType pdf(const hypergeometric_distribution<RealType, Policy>& dist, const std::uint64_t& x)
    {
       static const char* function = "boost::math::pdf(const hypergeometric_distribution<%1%>&, const %1%&)";
       RealType result = 0;
@@ -136,7 +137,7 @@ namespace boost { namespace math {
       BOOST_MATH_STD_USING
       static const char* function = "boost::math::pdf(const hypergeometric_distribution<%1%>&, const %1%&)";
       RealType r = static_cast<RealType>(x);
-      unsigned u = itrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type());
+      auto u = static_cast<std::uint64_t>(lltrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type()));
       if(u != r)
       {
          return boost::math::policies::raise_domain_error<RealType>(
@@ -146,7 +147,7 @@ namespace boost { namespace math {
    }
 
    template <class RealType, class Policy>
-   inline RealType cdf(const hypergeometric_distribution<RealType, Policy>& dist, const unsigned& x)
+   inline RealType cdf(const hypergeometric_distribution<RealType, Policy>& dist, const std::uint64_t& x)
    {
       static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";
       RealType result = 0;
@@ -165,7 +166,7 @@ namespace boost { namespace math {
       BOOST_MATH_STD_USING
       static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";
       RealType r = static_cast<RealType>(x);
-      unsigned u = itrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type());
+      auto u = static_cast<std::uint64_t>(lltrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type()));
       if(u != r)
       {
          return boost::math::policies::raise_domain_error<RealType>(
@@ -175,7 +176,7 @@ namespace boost { namespace math {
    }
 
    template <class RealType, class Policy>
-   inline RealType cdf(const complemented2_type<hypergeometric_distribution<RealType, Policy>, unsigned>& c)
+   inline RealType cdf(const complemented2_type<hypergeometric_distribution<RealType, Policy>, std::uint64_t>& c)
    {
       static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";
       RealType result = 0;
@@ -194,7 +195,7 @@ namespace boost { namespace math {
       BOOST_MATH_STD_USING
       static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";
       RealType r = static_cast<RealType>(c.param);
-      unsigned u = itrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type());
+      auto u = static_cast<std::uint64_t>(lltrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type()));
       if(u != r)
       {
          return boost::math::policies::raise_domain_error<RealType>(
@@ -208,11 +209,14 @@ namespace boost { namespace math {
    {
       BOOST_MATH_STD_USING // for ADL of std functions
 
-         // Checking function argument
-         RealType result = 0;
+      // Checking function argument
+      RealType result = 0;
       const char* function = "boost::math::quantile(const hypergeometric_distribution<%1%>&, %1%)";
-      if (false == dist.check_params(function, &result)) return result;
-      if(false == detail::check_probability(function, p, &result, Policy())) return result;
+      if (false == dist.check_params(function, &result)) 
+         return result;
+
+      if(false == detail::check_probability(function, p, &result, Policy())) 
+         return result;
 
       return static_cast<RealType>(detail::hypergeometric_quantile(p, RealType(1 - p), dist.defective(), dist.sample_count(), dist.total(), Policy()));
    } // quantile
@@ -225,8 +229,10 @@ namespace boost { namespace math {
       // Checking function argument
       RealType result = 0;
       const char* function = "quantile(const complemented2_type<hypergeometric_distribution<%1%>, %1%>&)";
-      if (false == c.dist.check_params(function, &result)) return result;
-      if(false == detail::check_probability(function, c.param, &result, Policy())) return result;
+      if (false == c.dist.check_params(function, &result)) 
+         return result;
+      if (false == detail::check_probability(function, c.param, &result, Policy()))
+         return result;
 
       return static_cast<RealType>(detail::hypergeometric_quantile(RealType(1 - c.param), c.param, c.dist.defective(), c.dist.sample_count(), c.dist.total(), Policy()));
    } // quantile