diff any/include/boost/math/distributions/detail/hypergeometric_pdf.hpp @ 160:cff480c41f97

Add some cross-platform Boost headers

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/any/include/boost/math/distributions/detail/hypergeometric_pdf.hpp	Sat Feb 16 16:31:25 2019 +0000
@@ -0,0 +1,488 @@
+// Copyright 2008 Gautam Sewani
+// Copyright 2008 John Maddock
+//
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_PDF_HPP
+#define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_PDF_HPP
+
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/special_functions/lanczos.hpp>
+#include <boost/math/special_functions/gamma.hpp>
+#include <boost/math/special_functions/pow.hpp>
+#include <boost/math/special_functions/prime.hpp>
+#include <boost/math/policies/error_handling.hpp>
+
+#ifdef BOOST_MATH_INSTRUMENT
+#include <typeinfo>
+#endif
+
+namespace boost{ namespace math{ namespace detail{
+
+template <class T, class Func>
+void bubble_down_one(T* first, T* last, Func f)
+{
+   using std::swap;
+   T* next = first;
+   ++next;
+   while((next != last) && (!f(*first, *next)))
+   {
+      swap(*first, *next);
+      ++first;
+      ++next;
+   }
+}
+
+template <class T>
+struct sort_functor
+{
+   sort_functor(const T* exponents) : m_exponents(exponents){}
+   bool operator()(int i, int j)
+   {
+      return m_exponents[i] > m_exponents[j];
+   }
+private:
+   const T* m_exponents;
+};
+
+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&)
+{
+   BOOST_MATH_STD_USING
+
+   BOOST_MATH_INSTRUMENT_FPU
+   BOOST_MATH_INSTRUMENT_VARIABLE(x);
+   BOOST_MATH_INSTRUMENT_VARIABLE(r);
+   BOOST_MATH_INSTRUMENT_VARIABLE(n);
+   BOOST_MATH_INSTRUMENT_VARIABLE(N);
+   BOOST_MATH_INSTRUMENT_VARIABLE(typeid(Lanczos).name());
+
+   T bases[9] = {
+      T(n) + static_cast<T>(Lanczos::g()) + 0.5f,
+      T(r) + static_cast<T>(Lanczos::g()) + 0.5f,
+      T(N - n) + static_cast<T>(Lanczos::g()) + 0.5f,
+      T(N - r) + static_cast<T>(Lanczos::g()) + 0.5f,
+      1 / (T(N) + static_cast<T>(Lanczos::g()) + 0.5f),
+      1 / (T(x) + static_cast<T>(Lanczos::g()) + 0.5f),
+      1 / (T(n - x) + static_cast<T>(Lanczos::g()) + 0.5f),
+      1 / (T(r - x) + static_cast<T>(Lanczos::g()) + 0.5f),
+      1 / (T(N - n - r + x) + static_cast<T>(Lanczos::g()) + 0.5f)
+   };
+   T exponents[9] = {
+      n + T(0.5f),
+      r + T(0.5f),
+      N - n + T(0.5f),
+      N - r + T(0.5f),
+      N + T(0.5f),
+      x + T(0.5f),
+      n - x + T(0.5f),
+      r - x + T(0.5f),
+      N - n - r + x + T(0.5f)
+   };
+   int base_e_factors[9] = {
+      -1, -1, -1, -1, 1, 1, 1, 1, 1
+   };
+   int sorted_indexes[9] = {
+      0, 1, 2, 3, 4, 5, 6, 7, 8
+   };
+#ifdef BOOST_MATH_INSTRUMENT
+   BOOST_MATH_INSTRUMENT_FPU
+   for(unsigned i = 0; i < 9; ++i)
+   {
+      BOOST_MATH_INSTRUMENT_VARIABLE(i);
+      BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
+      BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
+      BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
+      BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
+   }
+#endif
+   std::sort(sorted_indexes, sorted_indexes + 9, sort_functor<T>(exponents));
+#ifdef BOOST_MATH_INSTRUMENT
+   BOOST_MATH_INSTRUMENT_FPU
+   for(unsigned i = 0; i < 9; ++i)
+   {
+      BOOST_MATH_INSTRUMENT_VARIABLE(i);
+      BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
+      BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
+      BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
+      BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
+   }
+#endif
+
+   do{
+      exponents[sorted_indexes[0]] -= exponents[sorted_indexes[1]];
+      bases[sorted_indexes[1]] *= bases[sorted_indexes[0]];
+      if((bases[sorted_indexes[1]] < tools::min_value<T>()) && (exponents[sorted_indexes[1]] != 0))
+      {
+         return 0;
+      }
+      base_e_factors[sorted_indexes[1]] += base_e_factors[sorted_indexes[0]];
+      bubble_down_one(sorted_indexes, sorted_indexes + 9, sort_functor<T>(exponents));
+
+#ifdef BOOST_MATH_INSTRUMENT
+      for(unsigned i = 0; i < 9; ++i)
+      {
+         BOOST_MATH_INSTRUMENT_VARIABLE(i);
+         BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
+         BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
+         BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
+         BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
+      }
+#endif
+   }while(exponents[sorted_indexes[1]] > 1);
+
+   //
+   // Combine equal powers:
+   //
+   int j = 8;
+   while(exponents[sorted_indexes[j]] == 0) --j;
+   while(j)
+   {
+      while(j && (exponents[sorted_indexes[j-1]] == exponents[sorted_indexes[j]]))
+      {
+         bases[sorted_indexes[j-1]] *= bases[sorted_indexes[j]];
+         exponents[sorted_indexes[j]] = 0;
+         base_e_factors[sorted_indexes[j-1]] += base_e_factors[sorted_indexes[j]];
+         bubble_down_one(sorted_indexes + j, sorted_indexes + 9, sort_functor<T>(exponents));
+         --j;
+      }
+      --j;
+
+#ifdef BOOST_MATH_INSTRUMENT
+      BOOST_MATH_INSTRUMENT_VARIABLE(j);
+      for(unsigned i = 0; i < 9; ++i)
+      {
+         BOOST_MATH_INSTRUMENT_VARIABLE(i);
+         BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
+         BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
+         BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
+         BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
+      }
+#endif
+   }
+
+#ifdef BOOST_MATH_INSTRUMENT
+   BOOST_MATH_INSTRUMENT_FPU
+   for(unsigned i = 0; i < 9; ++i)
+   {
+      BOOST_MATH_INSTRUMENT_VARIABLE(i);
+      BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
+      BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
+      BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
+      BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
+   }
+#endif
+
+   T result;
+   BOOST_MATH_INSTRUMENT_VARIABLE(bases[sorted_indexes[0]] * exp(static_cast<T>(base_e_factors[sorted_indexes[0]])));
+   BOOST_MATH_INSTRUMENT_VARIABLE(exponents[sorted_indexes[0]]);
+   {
+      BOOST_FPU_EXCEPTION_GUARD
+      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)
+   {
+      BOOST_FPU_EXCEPTION_GUARD
+      if(result < tools::min_value<T>())
+         return 0; // short circuit further evaluation
+      if(exponents[sorted_indexes[i]] == 1)
+         result *= bases[sorted_indexes[i]] * exp(static_cast<T>(base_e_factors[sorted_indexes[i]]));
+      else if(exponents[sorted_indexes[i]] == 0.5f)
+         result *= sqrt(bases[sorted_indexes[i]] * exp(static_cast<T>(base_e_factors[sorted_indexes[i]])));
+      else
+         result *= pow(bases[sorted_indexes[i]] * exp(static_cast<T>(base_e_factors[sorted_indexes[i]])), exponents[sorted_indexes[i]]);
+   
+      BOOST_MATH_INSTRUMENT_VARIABLE(result);
+   }
+
+   result *= Lanczos::lanczos_sum_expG_scaled(static_cast<T>(n + 1))
+      * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(r + 1))
+      * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N - n + 1))
+      * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N - r + 1))
+      / 
+      ( Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N + 1))
+         * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(x + 1))
+         * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(n - x + 1))
+         * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(r - x + 1))
+         * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N - n - r + x + 1)));
+   
+   BOOST_MATH_INSTRUMENT_VARIABLE(result);
+   return result;
+}
+
+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)
+{
+   BOOST_MATH_STD_USING
+   return exp(
+      boost::math::lgamma(T(n + 1), pol)
+      + boost::math::lgamma(T(r + 1), pol)
+      + boost::math::lgamma(T(N - n + 1), pol)
+      + boost::math::lgamma(T(N - r + 1), pol)
+      - boost::math::lgamma(T(N + 1), pol)
+      - boost::math::lgamma(T(x + 1), pol)
+      - boost::math::lgamma(T(n - x + 1), pol)
+      - boost::math::lgamma(T(r - x + 1), pol)
+      - boost::math::lgamma(T(N - n - r + x + 1), pol));
+}
+
+template <class T>
+inline T integer_power(const T& x, int ex)
+{
+   if(ex < 0)
+      return 1 / integer_power(x, -ex);
+   switch(ex)
+   {
+   case 0:
+      return 1;
+   case 1:
+      return x;
+   case 2:
+      return x * x;
+   case 3:
+      return x * x * x;
+   case 4:
+      return boost::math::pow<4>(x);
+   case 5:
+      return boost::math::pow<5>(x);
+   case 6:
+      return boost::math::pow<6>(x);
+   case 7:
+      return boost::math::pow<7>(x);
+   case 8:
+      return boost::math::pow<8>(x);
+   }
+   BOOST_MATH_STD_USING
+#ifdef __SUNPRO_CC
+   return pow(x, T(ex));
+#else
+   return pow(x, ex);
+#endif
+}
+template <class T>
+struct hypergeometric_pdf_prime_loop_result_entry
+{
+   T value;
+   const hypergeometric_pdf_prime_loop_result_entry* next;
+};
+
+#ifdef BOOST_MSVC
+#pragma warning(push)
+#pragma warning(disable:4510 4512 4610)
+#endif
+
+struct hypergeometric_pdf_prime_loop_data
+{
+   const unsigned x;
+   const unsigned r;
+   const unsigned n;
+   const unsigned N;
+   unsigned prime_index;
+   unsigned current_prime;
+};
+
+#ifdef BOOST_MSVC
+#pragma warning(pop)
+#endif
+
+template <class T>
+T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hypergeometric_pdf_prime_loop_result_entry<T>& result)
+{
+   while(data.current_prime <= data.N)
+   {
+      unsigned base = data.current_prime;
+      int prime_powers = 0;
+      while(base <= data.N)
+      {
+         prime_powers += data.n / base;
+         prime_powers += data.r / base;
+         prime_powers += (data.N - data.n) / base;
+         prime_powers += (data.N - data.r) / base;
+         prime_powers -= data.N / base;
+         prime_powers -= data.x / base;
+         prime_powers -= (data.n - data.x) / base;
+         prime_powers -= (data.r - data.x) / base;
+         prime_powers -= (data.N - data.n - data.r + data.x) / base;
+         base *= data.current_prime;
+      }
+      if(prime_powers)
+      {
+         T p = integer_power<T>(static_cast<T>(data.current_prime), 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);
+            return hypergeometric_pdf_prime_loop_imp<T>(data, t);
+         }
+         if((p < 1) && (tools::min_value<T>() / p > result.value))
+         {
+            //
+            // The next calculation would underflow, use recursion
+            // to sidestep the issue:
+            //
+            hypergeometric_pdf_prime_loop_result_entry<T> t = { p, &result };
+            data.current_prime = prime(++data.prime_index);
+            return hypergeometric_pdf_prime_loop_imp<T>(data, t);
+         }
+         result.value *= p;
+      }
+      data.current_prime = prime(++data.prime_index);
+   }
+   //
+   // When we get to here we have run out of prime factors,
+   // the overall result is the product of all the partial
+   // results we have accumulated on the stack so far, these
+   // are in a linked list starting with "data.head" and ending
+   // with "result".
+   //
+   // All that remains is to multiply them together, taking
+   // care not to overflow or underflow.
+   //
+   // Enumerate partial results >= 1 in variable i
+   // and partial results < 1 in variable j:
+   //
+   hypergeometric_pdf_prime_loop_result_entry<T> const *i, *j;
+   i = &result;
+   while(i && i->value < 1)
+      i = i->next;
+   j = &result;
+   while(j && j->value >= 1)
+      j = j->next;
+
+   T prod = 1;
+
+   while(i || j)
+   {
+      while(i && ((prod <= 1) || (j == 0)))
+      {
+         prod *= i->value;
+         i = i->next;
+         while(i && i->value < 1)
+            i = i->next;
+      }
+      while(j && ((prod >= 1) || (i == 0)))
+      {
+         prod *= j->value;
+         j = j->next;
+         while(j && j->value >= 1)
+            j = j->next;
+      }
+   }
+
+   return prod;
+}
+
+template <class T, class Policy>
+inline T hypergeometric_pdf_prime_imp(unsigned x, unsigned r, unsigned n, unsigned 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&)
+{
+   BOOST_MATH_STD_USING
+   BOOST_ASSERT(N <= boost::math::max_factorial<T>::value);
+   T result = boost::math::unchecked_factorial<T>(n);
+   T num[3] = {
+      boost::math::unchecked_factorial<T>(r),
+      boost::math::unchecked_factorial<T>(N - n),
+      boost::math::unchecked_factorial<T>(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)
+   };
+   int i = 0;
+   int j = 0;
+   while((i < 3) || (j < 5))
+   {
+      while((j < 5) && ((result >= 1) || (i >= 3)))
+      {
+         result /= denom[j];
+         ++j;
+      }
+      while((i < 3) && ((result <= 1) || (j >= 5)))
+      {
+         result *= num[i];
+         ++i;
+      }
+   }
+   return result;
+}
+
+
+template <class T, class Policy>
+inline typename tools::promote_args<T>::type 
+   hypergeometric_pdf(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&)
+{
+   BOOST_FPU_EXCEPTION_GUARD
+   typedef typename tools::promote_args<T>::type result_type;
+   typedef typename policies::evaluation<result_type, Policy>::type value_type;
+   typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
+   typedef typename policies::normalise<
+      Policy, 
+      policies::promote_float<false>, 
+      policies::promote_double<false>, 
+      policies::discrete_quantile<>,
+      policies::assert_undefined<> >::type forwarding_policy;
+
+   value_type result;
+   if(N <= boost::math::max_factorial<value_type>::value)
+   {
+      //
+      // If N is small enough then we can evaluate the PDF via the factorials
+      // directly: table lookup of the factorials gives the best performance
+      // of the methods available:
+      //
+      result = detail::hypergeometric_pdf_factorial_imp<value_type>(x, r, n, N, forwarding_policy());
+   }
+   else if(N <= boost::math::prime(boost::math::max_prime - 1))
+   {
+      //
+      // If N is no larger than the largest prime number in our lookup table
+      // (104729) then we can use prime factorisation to evaluate the PDF,
+      // this is slow but accurate:
+      //
+      result = detail::hypergeometric_pdf_prime_imp<value_type>(x, r, n, N, forwarding_policy());
+   }
+   else
+   {
+      //
+      // Catch all case - use the lanczos approximation - where available - 
+      // to evaluate the ratio of factorials.  This is reasonably fast
+      // (almost as quick as using logarithmic evaluation in terms of lgamma)
+      // but only a few digits better in accuracy than using lgamma:
+      //
+      result = detail::hypergeometric_pdf_lanczos_imp(value_type(), x, r, n, N, evaluation_type(), forwarding_policy());
+   }
+
+   if(result > 1)
+   {
+      result = 1;
+   }
+   if(result < 0)
+   {
+      result = 0;
+   }
+
+   return policies::checked_narrowing_cast<result_type, forwarding_policy>(result, "boost::math::hypergeometric_pdf<%1%>(%1%,%1%,%1%,%1%)");
+}
+
+}}} // namespaces
+
+#endif
+