Added boost header

test/external/boost/math/distributions/detail/common_error_handling.hpp

@@ -0,0 +1,178 @@
+// Copyright John Maddock 2006, 2007.
+// Copyright Paul A. Bristow 2006, 2007.
+
+// 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_COMMON_ERROR_HANDLING_HPP
+#define BOOST_MATH_DISTRIBUTIONS_COMMON_ERROR_HANDLING_HPP
+
+#include <boost/math/policies/error_handling.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
+// using boost::math::isfinite;
+
+namespace boost{ namespace math{ namespace detail
+{
+
+template <class RealType, class Policy>
+inline bool check_probability(const char* function, RealType const& prob, RealType* result, const Policy& pol)
+{
+   if((prob < 0) || (prob > 1) || !(boost::math::isfinite)(prob))
+   {
+      *result = policies::raise_domain_error<RealType>(
+         function,
+         "Probability argument is %1%, but must be >= 0 and <= 1 !", prob, pol);
+      return false;
+   }
+   return true;
+}
+
+template <class RealType, class Policy>
+inline bool check_df(const char* function, RealType const& df, RealType* result, const Policy& pol)
+{
+   if((df <= 0) || !(boost::math::isfinite)(df))
+   {
+      *result = policies::raise_domain_error<RealType>(
+         function,
+         "Degrees of freedom argument is %1%, but must be > 0 !", df, pol);
+      return false;
+   }
+   return true;
+}
+
+template <class RealType, class Policy>
+inline bool check_scale(
+      const char* function,
+      RealType scale,
+      RealType* result,
+      const Policy& pol)
+{
+   if((scale <= 0) || !(boost::math::isfinite)(scale))
+   { // Assume scale == 0 is NOT valid for any distribution.
+      *result = policies::raise_domain_error<RealType>(
+         function,
+         "Scale parameter is %1%, but must be > 0 !", scale, pol);
+      return false;
+   }
+   return true;
+}
+
+template <class RealType, class Policy>
+inline bool check_location(
+      const char* function,
+      RealType location,
+      RealType* result,
+      const Policy& pol)
+{
+   if(!(boost::math::isfinite)(location))
+   {
+      *result = policies::raise_domain_error<RealType>(
+         function,
+         "Location parameter is %1%, but must be finite!", location, pol);
+      return false;
+   }
+   return true;
+}
+
+template <class RealType, class Policy>
+inline bool check_x(
+      const char* function,
+      RealType x,
+      RealType* result,
+      const Policy& pol)
+{
+   if(!(boost::math::isfinite)(x))
+   {
+      *result = policies::raise_domain_error<RealType>(
+         function,
+         "Random variate x is %1%, but must be finite!", x, pol);
+      return false;
+   }
+   return true;
+   // Note that this test catches both infinity and NaN.
+   // Some special cases permit x to be infinite, so these must be tested 1st,
+   // leaving this test to catch any NaNs.  see Normal and cauchy for example.
+} // bool check_x
+
+template <class RealType, class Policy>
+inline bool check_x_gt0(
+      const char* function,
+      RealType x,
+      RealType* result,
+      const Policy& pol)
+{
+   if(x <= 0)
+   {
+      *result = policies::raise_domain_error<RealType>(
+         function,
+         "Random variate x is %1%, but must be > 0!", x, pol);
+      return false;
+   }
+
+   return true;
+   // Note that this test catches both infinity and NaN.
+   // Some special cases permit x to be infinite, so these must be tested 1st,
+   // leaving this test to catch any NaNs.  See Normal and cauchy for example.
+} // bool check_x_gt0
+
+template <class RealType, class Policy>
+inline bool check_positive_x(
+      const char* function,
+      RealType x,
+      RealType* result,
+      const Policy& pol)
+{
+   if(!(boost::math::isfinite)(x) || (x < 0))
+   {
+      *result = policies::raise_domain_error<RealType>(
+         function,
+         "Random variate x is %1%, but must be finite and >= 0!", x, pol);
+      return false;
+   }
+   return true;
+   // Note that this test catches both infinity and NaN.
+   // Some special cases permit x to be infinite, so these must be tested 1st,
+   // leaving this test to catch any NaNs.  see Normal and cauchy for example.
+}
+
+template <class RealType, class Policy>
+inline bool check_non_centrality(
+      const char* function,
+      RealType ncp,
+      RealType* result,
+      const Policy& pol)
+{
+   if((ncp < 0) || !(boost::math::isfinite)(ncp))
+   { // Assume scale == 0 is NOT valid for any distribution.
+      *result = policies::raise_domain_error<RealType>(
+         function,
+         "Non centrality parameter is %1%, but must be > 0 !", ncp, pol);
+      return false;
+   }
+   return true;
+}
+
+template <class RealType, class Policy>
+inline bool check_finite(
+      const char* function,
+      RealType x,
+      RealType* result,
+      const Policy& pol)
+{
+   if(!(boost::math::isfinite)(x))
+   { // Assume scale == 0 is NOT valid for any distribution.
+      *result = policies::raise_domain_error<RealType>(
+         function,
+         "Parameter is %1%, but must be finite !", x, pol);
+      return false;
+   }
+   return true;
+}
+
+} // namespace detail
+} // namespace math
+} // namespace boost
+
+#endif // BOOST_MATH_DISTRIBUTIONS_COMMON_ERROR_HANDLING_HPP

test/external/boost/math/distributions/detail/derived_accessors.hpp

@@ -0,0 +1,163 @@
+//  Copyright John Maddock 2006.
+//  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_STATS_DERIVED_HPP
+#define BOOST_STATS_DERIVED_HPP
+
+// This file implements various common properties of distributions
+// that can be implemented in terms of other properties:
+// variance OR standard deviation (see note below),
+// hazard, cumulative hazard (chf), coefficient_of_variation.
+//
+// Note that while both variance and standard_deviation are provided
+// here, each distribution MUST SPECIALIZE AT LEAST ONE OF THESE
+// otherwise these two versions will just call each other over and over
+// until stack space runs out ...
+
+// Of course there may be more efficient means of implementing these
+// that are specific to a particular distribution, but these generic
+// versions give these properties "for free" with most distributions.
+//
+// In order to make use of this header, it must be included AT THE END
+// of the distribution header, AFTER the distribution and its core
+// property accessors have been defined: this is so that compilers
+// that implement 2-phase lookup and early-type-checking of templates
+// can find the definitions refered to herein.
+//
+
+#include <boost/type_traits/is_same.hpp>
+#include <boost/static_assert.hpp>
+
+#ifdef BOOST_MSVC
+# pragma warning(push)
+# pragma warning(disable: 4723) // potential divide by 0
+// Suppressing spurious warning in coefficient_of_variation
+#endif
+
+namespace boost{ namespace math{
+
+template <class Distribution>
+typename Distribution::value_type variance(const Distribution& dist);
+
+template <class Distribution>
+inline typename Distribution::value_type standard_deviation(const Distribution& dist)
+{
+   BOOST_MATH_STD_USING  // ADL of sqrt.
+   return sqrt(variance(dist));
+}
+
+template <class Distribution>
+inline typename Distribution::value_type variance(const Distribution& dist)
+{
+   typename Distribution::value_type result = standard_deviation(dist);
+   return result * result;
+}
+
+template <class Distribution, class RealType>
+inline typename Distribution::value_type hazard(const Distribution& dist, const RealType& x)
+{ // hazard function
+  // http://www.itl.nist.gov/div898/handbook/eda/section3/eda362.htm#HAZ
+   typedef typename Distribution::value_type value_type;
+   typedef typename Distribution::policy_type policy_type;
+   value_type p = cdf(complement(dist, x));
+   value_type d = pdf(dist, x);
+   if(d > p * tools::max_value<value_type>())
+      return policies::raise_overflow_error<value_type>(
+      "boost::math::hazard(const Distribution&, %1%)", 0, policy_type());
+   if(d == 0)
+   {
+      // This protects against 0/0, but is it the right thing to do?
+      return 0;
+   }
+   return d / p;
+}
+
+template <class Distribution, class RealType>
+inline typename Distribution::value_type chf(const Distribution& dist, const RealType& x)
+{ // cumulative hazard function.
+  // http://www.itl.nist.gov/div898/handbook/eda/section3/eda362.htm#HAZ
+   BOOST_MATH_STD_USING
+   return -log(cdf(complement(dist, x)));
+}
+
+template <class Distribution>
+inline typename Distribution::value_type coefficient_of_variation(const Distribution& dist)
+{
+   typedef typename Distribution::value_type value_type;
+   typedef typename Distribution::policy_type policy_type;
+
+   using std::abs;
+
+   value_type m = mean(dist);
+   value_type d = standard_deviation(dist);
+   if((abs(m) < 1) && (d > abs(m) * tools::max_value<value_type>()))
+   { // Checks too that m is not zero,
+      return policies::raise_overflow_error<value_type>("boost::math::coefficient_of_variation(const Distribution&, %1%)", 0, policy_type());
+   }
+   return d / m; // so MSVC warning on zerodivide is spurious, and suppressed.
+}
+//
+// Next follow overloads of some of the standard accessors with mixed
+// argument types. We just use a typecast to forward on to the "real"
+// implementation with all arguments of the same type:
+//
+template <class Distribution, class RealType>
+inline typename Distribution::value_type pdf(const Distribution& dist, const RealType& x)
+{
+   typedef typename Distribution::value_type value_type;
+   return pdf(dist, static_cast<value_type>(x));
+}
+template <class Distribution, class RealType>
+inline typename Distribution::value_type cdf(const Distribution& dist, const RealType& x)
+{
+   typedef typename Distribution::value_type value_type;
+   return cdf(dist, static_cast<value_type>(x));
+}
+template <class Distribution, class RealType>
+inline typename Distribution::value_type quantile(const Distribution& dist, const RealType& x)
+{
+   typedef typename Distribution::value_type value_type;
+   return quantile(dist, static_cast<value_type>(x));
+}
+/*
+template <class Distribution, class RealType>
+inline typename Distribution::value_type chf(const Distribution& dist, const RealType& x)
+{
+   typedef typename Distribution::value_type value_type;
+   return chf(dist, static_cast<value_type>(x));
+}
+*/
+template <class Distribution, class RealType>
+inline typename Distribution::value_type cdf(const complemented2_type<Distribution, RealType>& c)
+{
+   typedef typename Distribution::value_type value_type;
+   return cdf(complement(c.dist, static_cast<value_type>(c.param)));
+}
+
+template <class Distribution, class RealType>
+inline typename Distribution::value_type quantile(const complemented2_type<Distribution, RealType>& c)
+{
+   typedef typename Distribution::value_type value_type;
+   return quantile(complement(c.dist, static_cast<value_type>(c.param)));
+}
+
+template <class Dist>
+inline typename Dist::value_type median(const Dist& d)
+{ // median - default definition for those distributions for which a
+  // simple closed form is not known,
+  // and for which a domain_error and/or NaN generating function is NOT defined.
+  typedef typename Dist::value_type value_type;
+  return quantile(d, static_cast<value_type>(0.5f));
+}
+
+} // namespace math
+} // namespace boost
+
+
+#ifdef BOOST_MSVC
+# pragma warning(pop)
+#endif
+
+#endif // BOOST_STATS_DERIVED_HPP

test/external/boost/math/distributions/detail/generic_mode.hpp

@@ -0,0 +1,149 @@
+// Copyright John Maddock 2008.
+
+// 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_MODE_HPP
+#define BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP
+
+#include <boost/math/tools/minima.hpp> // function minimization for mode
+#include <boost/math/policies/error_handling.hpp>
+#include <boost/math/distributions/fwd.hpp>
+
+namespace boost{ namespace math{ namespace detail{
+
+template <class Dist>
+struct pdf_minimizer
+{
+   pdf_minimizer(const Dist& d)
+      : dist(d) {}
+
+   typename Dist::value_type operator()(const typename Dist::value_type& x)
+   {
+      return -pdf(dist, x);
+   }
+private:
+   Dist dist;
+};
+
+template <class Dist>
+typename Dist::value_type generic_find_mode(const Dist& dist, typename Dist::value_type guess, const char* function, typename Dist::value_type step = 0)
+{
+   BOOST_MATH_STD_USING
+   typedef typename Dist::value_type value_type;
+   typedef typename Dist::policy_type policy_type;
+   //
+   // Need to begin by bracketing the maxima of the PDF:
+   //
+   value_type maxval;
+   value_type upper_bound = guess;
+   value_type lower_bound;
+   value_type v = pdf(dist, guess);
+   if(v == 0)
+   {
+      //
+      // Oops we don't know how to handle this, or even in which
+      // direction we should move in, treat as an evaluation error:
+      //
+      policies::raise_evaluation_error(
+         function, 
+         "Could not locate a starting location for the search for the mode, original guess was %1%", guess, policy_type());
+   }
+   do
+   {
+      maxval = v;
+      if(step != 0)
+         upper_bound += step;
+      else
+         upper_bound *= 2;
+      v = pdf(dist, upper_bound);
+   }while(maxval < v);
+
+   lower_bound = upper_bound;
+   do
+   {
+      maxval = v;
+      if(step != 0)
+         lower_bound -= step;
+      else
+         lower_bound /= 2;
+      v = pdf(dist, lower_bound);
+   }while(maxval < v);
+
+   boost::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>();
+
+   value_type result = tools::brent_find_minima(
+      pdf_minimizer<Dist>(dist), 
+      lower_bound, 
+      upper_bound, 
+      policies::digits<value_type, policy_type>(), 
+      max_iter).first;
+   if(max_iter >= policies::get_max_root_iterations<policy_type>())
+   {
+      return policies::raise_evaluation_error<value_type>(
+         function, 
+         "Unable to locate solution in a reasonable time:"
+         " either there is no answer to the mode of the distribution"
+         " or the answer is infinite.  Current best guess is %1%", result, policy_type());
+   }
+   return result;
+}
+//
+// As above,but confined to the interval [0,1]:
+//
+template <class Dist>
+typename Dist::value_type generic_find_mode_01(const Dist& dist, typename Dist::value_type guess, const char* function)
+{
+   BOOST_MATH_STD_USING
+   typedef typename Dist::value_type value_type;
+   typedef typename Dist::policy_type policy_type;
+   //
+   // Need to begin by bracketing the maxima of the PDF:
+   //
+   value_type maxval;
+   value_type upper_bound = guess;
+   value_type lower_bound;
+   value_type v = pdf(dist, guess);
+   do
+   {
+      maxval = v;
+      upper_bound = 1 - (1 - upper_bound) / 2;
+      if(upper_bound == 1)
+         return 1;
+      v = pdf(dist, upper_bound);
+   }while(maxval < v);
+
+   lower_bound = upper_bound;
+   do
+   {
+      maxval = v;
+      lower_bound /= 2;
+      if(lower_bound < tools::min_value<value_type>())
+         return 0;
+      v = pdf(dist, lower_bound);
+   }while(maxval < v);
+
+   boost::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>();
+
+   value_type result = tools::brent_find_minima(
+      pdf_minimizer<Dist>(dist), 
+      lower_bound, 
+      upper_bound, 
+      policies::digits<value_type, policy_type>(), 
+      max_iter).first;
+   if(max_iter >= policies::get_max_root_iterations<policy_type>())
+   {
+      return policies::raise_evaluation_error<value_type>(
+         function, 
+         "Unable to locate solution in a reasonable time:"
+         " either there is no answer to the mode of the distribution"
+         " or the answer is infinite.  Current best guess is %1%", result, policy_type());
+   }
+   return result;
+}
+
+}}} // namespaces
+
+#endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP

test/external/boost/math/distributions/detail/generic_quantile.hpp

@@ -0,0 +1,91 @@
+//  Copyright John Maddock 2008.
+//  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_DISTIBUTIONS_DETAIL_GENERIC_QUANTILE_HPP
+#define BOOST_MATH_DISTIBUTIONS_DETAIL_GENERIC_QUANTILE_HPP
+
+namespace boost{ namespace math{ namespace detail{
+
+template <class Dist>
+struct generic_quantile_finder
+{
+   typedef typename Dist::value_type value_type;
+   typedef typename Dist::policy_type policy_type;
+
+   generic_quantile_finder(const Dist& d, value_type t, bool c)
+      : dist(d), target(t), comp(c) {}
+
+   value_type operator()(const value_type& x)
+   {
+      return comp ?
+         value_type(target - cdf(complement(dist, x)))
+         : value_type(cdf(dist, x) - target);
+   }
+
+private:
+   Dist dist;
+   value_type target;
+   bool comp;
+};
+
+template <class T, class Policy>
+inline T check_range_result(const T& x, const Policy& pol, const char* function)
+{
+   if((x >= 0) && (x < tools::min_value<T>()))
+      return policies::raise_underflow_error<T>(function, 0, pol);
+   if(x <= -tools::max_value<T>())
+      return -policies::raise_overflow_error<T>(function, 0, pol);
+   if(x >= tools::max_value<T>())
+      return policies::raise_overflow_error<T>(function, 0, pol);
+   return x;
+}
+
+template <class Dist>
+typename Dist::value_type generic_quantile(const Dist& dist, const typename Dist::value_type& p, const typename Dist::value_type& guess, bool comp, const char* function)
+{
+   typedef typename Dist::value_type value_type;
+   typedef typename Dist::policy_type policy_type;
+   typedef typename policies::normalise<
+      policy_type, 
+      policies::promote_float<false>, 
+      policies::promote_double<false>, 
+      policies::discrete_quantile<>,
+      policies::assert_undefined<> >::type forwarding_policy;
+
+   //
+   // Special cases first:
+   //
+   if(p == 0)
+   {
+      return comp
+      ? check_range_result(range(dist).second, forwarding_policy(), function)
+      : check_range_result(range(dist).first, forwarding_policy(), function);
+   }
+   if(p == 1)
+   {
+      return !comp
+      ? check_range_result(range(dist).second, forwarding_policy(), function)
+      : check_range_result(range(dist).first, forwarding_policy(), function);
+   }
+
+   generic_quantile_finder<Dist> f(dist, p, comp);
+   tools::eps_tolerance<value_type> tol(policies::digits<value_type, forwarding_policy>() - 3);
+   boost::uintmax_t max_iter = policies::get_max_root_iterations<forwarding_policy>();
+   std::pair<value_type, value_type> ir = tools::bracket_and_solve_root(
+      f, guess, value_type(2), true, tol, max_iter, forwarding_policy());
+   value_type result = ir.first + (ir.second - ir.first) / 2;
+   if(max_iter >= policies::get_max_root_iterations<forwarding_policy>())
+   {
+      policies::raise_evaluation_error<value_type>(function, "Unable to locate solution in a reasonable time:"
+         " either there is no answer to quantile"
+         " or the answer is infinite.  Current best guess is %1%", result, forwarding_policy());
+   }
+   return result;
+}
+
+}}} // namespaces
+
+#endif // BOOST_MATH_DISTIBUTIONS_DETAIL_GENERIC_QUANTILE_HPP
+

test/external/boost/math/distributions/detail/hypergeometric_cdf.hpp

@@ -0,0 +1,100 @@
+// 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_CDF_HPP
+#define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_CDF_HPP
+
+#include <boost/math/policies/error_handling.hpp>
+#include <boost/math/distributions/detail/hypergeometric_pdf.hpp>
+
+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)
+   {
+#ifdef BOOST_MSVC
+#  pragma warning(push)
+#  pragma warning(disable:4267)
+#endif
+      BOOST_MATH_STD_USING
+      T result = 0;
+      T mode = floor(T(r + 1) * T(n + 1) / (N + 2));
+      if(x < mode)
+      {
+         result = hypergeometric_pdf<T>(x, r, n, N, pol);
+         T diff = result;
+         unsigned lower_limit = static_cast<unsigned>((std::max)(0, (int)(n + r) - (int)(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));
+            result += diff;
+            BOOST_MATH_INSTRUMENT_VARIABLE(x);
+            BOOST_MATH_INSTRUMENT_VARIABLE(diff);
+            BOOST_MATH_INSTRUMENT_VARIABLE(result);
+            if(x == lower_limit)
+               break;
+            --x;
+         }
+      }
+      else
+      {
+         invert = !invert;
+         unsigned upper_limit = (std::min)(r, n);
+         if(x != upper_limit)
+         {
+            ++x;
+            result = hypergeometric_pdf<T>(x, r, n, N, pol);
+            T diff = result;
+            while((x <= upper_limit) && (diff > (invert ? T(1) : result) * tools::epsilon<T>()))
+            {
+               diff = T(n - x) * T(r - x) * diff / (T(x + 1) * T((N + x + 1) - n - r));
+               result += diff;
+               ++x;
+               BOOST_MATH_INSTRUMENT_VARIABLE(x);
+               BOOST_MATH_INSTRUMENT_VARIABLE(diff);
+               BOOST_MATH_INSTRUMENT_VARIABLE(result);
+            }
+         }
+      }
+      if(invert)
+         result = 1 - result;
+      return result;
+#ifdef BOOST_MSVC
+#  pragma warning(pop)
+#endif
+   }
+
+   template <class T, class Policy>
+   inline T hypergeometric_cdf(unsigned x, unsigned r, unsigned n, unsigned N, bool invert, 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 policies::normalise<
+         Policy, 
+         policies::promote_float<false>, 
+         policies::promote_double<false>, 
+         policies::discrete_quantile<>,
+         policies::assert_undefined<> >::type forwarding_policy;
+
+      value_type result;
+      result = detail::hypergeometric_cdf_imp<value_type>(x, r, n, N, invert, 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_cdf<%1%>(%1%,%1%,%1%,%1%)");
+   }
+
+}}} // namespaces
+
+#endif
+

test/external/boost/math/distributions/detail/hypergeometric_pdf.hpp

@@ -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) + Lanczos::g() + 0.5f,
+      T(r) + Lanczos::g() + 0.5f,
+      T(N - n) + Lanczos::g() + 0.5f,
+      T(N - r) + Lanczos::g() + 0.5f,
+      1 / (T(N) + Lanczos::g() + 0.5f),
+      1 / (T(x) + Lanczos::g() + 0.5f),
+      1 / (T(n - x) + Lanczos::g() + 0.5f),
+      1 / (T(r - x) + Lanczos::g() + 0.5f),
+      1 / (T(N - n - r + x) + 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>(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
+

test/external/boost/math/distributions/detail/hypergeometric_quantile.hpp

@@ -0,0 +1,199 @@
+// 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_QUANTILE_HPP
+#define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_QUANTILE_HPP
+
+#include <boost/math/policies/error_handling.hpp>
+#include <boost/math/distributions/detail/hypergeometric_pdf.hpp>
+
+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>&)
+{
+   if((p < cum * fudge_factor) && (x != lbound))
+   {
+      BOOST_MATH_INSTRUMENT_VARIABLE(x-1);
+      return --x;
+   }
+   return x;
+}
+
+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>&)
+{
+   if((cum < p * fudge_factor) && (x != ubound))
+   {
+      BOOST_MATH_INSTRUMENT_VARIABLE(x+1);
+      return ++x;
+   }
+   return x;
+}
+
+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>&)
+{
+   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>&)
+{
+   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>&)
+{
+   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>&)
+{
+   if((q * fudge_factor > cum) && (x != lbound))
+   {
+      BOOST_MATH_INSTRUMENT_VARIABLE(x-1);
+      return --x;
+   }
+   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_up>&)
+{
+   if((q < cum * fudge_factor) && (x != ubound))
+   {
+      BOOST_MATH_INSTRUMENT_VARIABLE(x+1);
+      return ++x;
+   }
+   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_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>&)
+{
+   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>&)
+{
+   return x;
+}
+
+template <class T, class Policy>
+unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned N, const Policy& pol)
+{
+#ifdef BOOST_MSVC
+#  pragma warning(push)
+#  pragma warning(disable:4267)
+#endif
+   typedef typename Policy::discrete_quantile_type discrete_quantile_type;
+   BOOST_MATH_STD_USING
+   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, (int)(n + r) - (int)(N)));
+   unsigned lim = (std::min)(r, n);
+
+   BOOST_MATH_INSTRUMENT_VARIABLE(p);
+   BOOST_MATH_INSTRUMENT_VARIABLE(q);
+   BOOST_MATH_INSTRUMENT_VARIABLE(r);
+   BOOST_MATH_INSTRUMENT_VARIABLE(n);
+   BOOST_MATH_INSTRUMENT_VARIABLE(N);
+   BOOST_MATH_INSTRUMENT_VARIABLE(fudge_factor);
+   BOOST_MATH_INSTRUMENT_VARIABLE(base);
+   BOOST_MATH_INSTRUMENT_VARIABLE(lim);
+
+   if(p <= 0.5)
+   {
+      unsigned x = base;
+      result = hypergeometric_pdf<T>(x, r, n, N, pol);
+      T diff = result;
+      while(result < p)
+      {
+         diff = (diff > tools::min_value<T>() * 8) 
+            ? T(n - x) * T(r - x) * diff / (T(x + 1) * T(N + x + 1 - n - r))
+            : hypergeometric_pdf<T>(x + 1, r, n, N, pol);
+         if(result + diff / 2 > p)
+            break;
+         ++x;
+         result += diff;
+#ifdef BOOST_MATH_INSTRUMENT
+         if(diff != 0)
+         {
+            BOOST_MATH_INSTRUMENT_VARIABLE(x);
+            BOOST_MATH_INSTRUMENT_VARIABLE(diff);
+            BOOST_MATH_INSTRUMENT_VARIABLE(result);
+         }
+#endif
+      }
+      return round_x_from_p(x, p, result, fudge_factor, base, lim, discrete_quantile_type());
+   }
+   else
+   {
+      unsigned x = lim;
+      result = 0;
+      T diff = hypergeometric_pdf<T>(x, r, n, N, pol);
+      while(result + diff / 2 < q)
+      {
+         result += diff;
+         diff = (diff > tools::min_value<T>() * 8)
+            ? x * T(N + x - n - r) * diff / (T(1 + n - x) * T(1 + r - x))
+            : hypergeometric_pdf<T>(x - 1, r, n, N, pol);
+         --x;
+#ifdef BOOST_MATH_INSTRUMENT
+         if(diff != 0)
+         {
+            BOOST_MATH_INSTRUMENT_VARIABLE(x);
+            BOOST_MATH_INSTRUMENT_VARIABLE(diff);
+            BOOST_MATH_INSTRUMENT_VARIABLE(result);
+         }
+#endif
+      }
+      return round_x_from_q(x, q, result, fudge_factor, base, lim, discrete_quantile_type());
+   }
+#ifdef BOOST_MSVC
+#  pragma warning(pop)
+#endif
+}
+
+template <class T, class Policy>
+inline unsigned hypergeometric_quantile(T p, T q, 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 policies::normalise<
+      Policy, 
+      policies::promote_float<false>, 
+      policies::promote_double<false>, 
+      policies::assert_undefined<> >::type forwarding_policy;
+
+   return detail::hypergeometric_quantile_imp<value_type>(p, q, r, n, N, forwarding_policy());
+}
+
+}}} // namespaces
+
+#endif
+

test/external/boost/math/distributions/detail/inv_discrete_quantile.hpp

@@ -0,0 +1,481 @@
+//  Copyright John Maddock 2007.
+//  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_INV_DISCRETE_QUANTILE
+#define BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE
+
+#include <algorithm>
+
+namespace boost{ namespace math{ namespace detail{
+
+//
+// Functor for root finding algorithm:
+//
+template <class Dist>
+struct distribution_quantile_finder
+{
+   typedef typename Dist::value_type value_type;
+   typedef typename Dist::policy_type policy_type;
+
+   distribution_quantile_finder(const Dist d, value_type p, value_type q)
+      : dist(d), target(p < q ? p : q), comp(p < q ? false : true) {}
+
+   value_type operator()(value_type const& x)
+   {
+      return comp ? value_type(target - cdf(complement(dist, x))) : value_type(cdf(dist, x) - target);
+   }
+
+private:
+   Dist dist;
+   value_type target;
+   bool comp;
+};
+//
+// The purpose of adjust_bounds, is to toggle the last bit of the
+// range so that both ends round to the same integer, if possible.
+// If they do both round the same then we terminate the search
+// for the root *very* quickly when finding an integer result.
+// At the point that this function is called we know that "a" is
+// below the root and "b" above it, so this change can not result
+// in the root no longer being bracketed.
+//
+template <class Real, class Tol>
+void adjust_bounds(Real& /* a */, Real& /* b */, Tol const& /* tol */){}
+
+template <class Real>
+void adjust_bounds(Real& /* a */, Real& b, tools::equal_floor const& /* tol */)
+{
+   BOOST_MATH_STD_USING
+   b -= tools::epsilon<Real>() * b;
+}
+
+template <class Real>
+void adjust_bounds(Real& a, Real& /* b */, tools::equal_ceil const& /* tol */)
+{
+   BOOST_MATH_STD_USING
+   a += tools::epsilon<Real>() * a;
+}
+
+template <class Real>
+void adjust_bounds(Real& a, Real& b, tools::equal_nearest_integer const& /* tol */)
+{
+   BOOST_MATH_STD_USING
+   a += tools::epsilon<Real>() * a;
+   b -= tools::epsilon<Real>() * b;
+}
+//
+// This is where all the work is done:
+//
+template <class Dist, class Tolerance>
+typename Dist::value_type 
+   do_inverse_discrete_quantile(
+      const Dist& dist,
+      const typename Dist::value_type& p,
+      const typename Dist::value_type& q,
+      typename Dist::value_type guess,
+      const typename Dist::value_type& multiplier,
+      typename Dist::value_type adder,
+      const Tolerance& tol,
+      boost::uintmax_t& max_iter)
+{
+   typedef typename Dist::value_type value_type;
+   typedef typename Dist::policy_type policy_type;
+
+   static const char* function = "boost::math::do_inverse_discrete_quantile<%1%>";
+
+   BOOST_MATH_STD_USING
+
+   distribution_quantile_finder<Dist> f(dist, p, q);
+   //
+   // Max bounds of the distribution:
+   //
+   value_type min_bound, max_bound;
+   boost::math::tie(min_bound, max_bound) = support(dist);
+
+   if(guess > max_bound)
+      guess = max_bound;
+   if(guess < min_bound)
+      guess = min_bound;
+
+   value_type fa = f(guess);
+   boost::uintmax_t count = max_iter - 1;
+   value_type fb(fa), a(guess), b =0; // Compiler warning C4701: potentially uninitialized local variable 'b' used
+
+   if(fa == 0)
+      return guess;
+
+   //
+   // For small expected results, just use a linear search:
+   //
+   if(guess < 10)
+   {
+      b = a;
+      while((a < 10) && (fa * fb >= 0))
+      {
+         if(fb <= 0)
+         {
+            a = b;
+            b = a + 1;
+            if(b > max_bound)
+               b = max_bound;
+            fb = f(b);
+            --count;
+            if(fb == 0)
+               return b;
+         }
+         else
+         {
+            b = a;
+            a = (std::max)(value_type(b - 1), value_type(0));
+            if(a < min_bound)
+               a = min_bound;
+            fa = f(a);
+            --count;
+            if(fa == 0)
+               return a;
+         }
+      }
+   }
+   //
+   // Try and bracket using a couple of additions first, 
+   // we're assuming that "guess" is likely to be accurate
+   // to the nearest int or so:
+   //
+   else if(adder != 0)
+   {
+      //
+      // If we're looking for a large result, then bump "adder" up
+      // by a bit to increase our chances of bracketing the root:
+      //
+      //adder = (std::max)(adder, 0.001f * guess);
+      if(fa < 0)
+      {
+         b = a + adder;
+         if(b > max_bound)
+            b = max_bound;
+      }
+      else
+      {
+         b = (std::max)(value_type(a - adder), value_type(0));
+         if(b < min_bound)
+            b = min_bound;
+      }
+      fb = f(b);
+      --count;
+      if(fb == 0)
+         return b;
+      if(count && (fa * fb >= 0))
+      {
+         //
+         // We didn't bracket the root, try 
+         // once more:
+         //
+         a = b;
+         fa = fb;
+         if(fa < 0)
+         {
+            b = a + adder;
+            if(b > max_bound)
+               b = max_bound;
+         }
+         else
+         {
+            b = (std::max)(value_type(a - adder), value_type(0));
+            if(b < min_bound)
+               b = min_bound;
+         }
+         fb = f(b);
+         --count;
+      }
+      if(a > b)
+      {
+         using std::swap;
+         swap(a, b);
+         swap(fa, fb);
+      }
+   }
+   //
+   // If the root hasn't been bracketed yet, try again
+   // using the multiplier this time:
+   //
+   if((boost::math::sign)(fb) == (boost::math::sign)(fa))
+   {
+      if(fa < 0)
+      {
+         //
+         // Zero is to the right of x2, so walk upwards
+         // until we find it:
+         //
+         while((boost::math::sign)(fb) == (boost::math::sign)(fa))
+         {
+            if(count == 0)
+               policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, policy_type());
+            a = b;
+            fa = fb;
+            b *= multiplier;
+            if(b > max_bound)
+               b = max_bound;
+            fb = f(b);
+            --count;
+            BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);
+         }
+      }
+      else
+      {
+         //
+         // Zero is to the left of a, so walk downwards
+         // until we find it:
+         //
+         while((boost::math::sign)(fb) == (boost::math::sign)(fa))
+         {
+            if(fabs(a) < tools::min_value<value_type>())
+            {
+               // Escape route just in case the answer is zero!
+               max_iter -= count;
+               max_iter += 1;
+               return 0;
+            }
+            if(count == 0)
+               policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, policy_type());
+            b = a;
+            fb = fa;
+            a /= multiplier;
+            if(a < min_bound)
+               a = min_bound;
+            fa = f(a);
+            --count;
+            BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);
+         }
+      }
+   }
+   max_iter -= count;
+   if(fa == 0)
+      return a;
+   if(fb == 0)
+      return b;
+   //
+   // Adjust bounds so that if we're looking for an integer
+   // result, then both ends round the same way:
+   //
+   adjust_bounds(a, b, tol);
+   //
+   // We don't want zero or denorm lower bounds:
+   //
+   if(a < tools::min_value<value_type>())
+      a = tools::min_value<value_type>();
+   //
+   // Go ahead and find the root:
+   //
+   std::pair<value_type, value_type> r = toms748_solve(f, a, b, fa, fb, tol, count, policy_type());
+   max_iter += count;
+   BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count);
+   return (r.first + r.second) / 2;
+}
+//
+// Now finally are the public API functions.
+// There is one overload for each policy,
+// each one is responsible for selecting the correct
+// termination condition, and rounding the result
+// to an int where required.
+//
+template <class Dist>
+inline typename Dist::value_type 
+   inverse_discrete_quantile(
+      const Dist& dist,
+      const typename Dist::value_type& p,
+      const typename Dist::value_type& q,
+      const typename Dist::value_type& guess,
+      const typename Dist::value_type& multiplier,
+      const typename Dist::value_type& adder,
+      const policies::discrete_quantile<policies::real>&,
+      boost::uintmax_t& max_iter)
+{
+   if(p <= pdf(dist, 0))
+      return 0;
+   return do_inverse_discrete_quantile(
+      dist, 
+      p, 
+      q,
+      guess, 
+      multiplier, 
+      adder, 
+      tools::eps_tolerance<typename Dist::value_type>(policies::digits<typename Dist::value_type, typename Dist::policy_type>()),
+      max_iter);
+}
+
+template <class Dist>
+inline typename Dist::value_type 
+   inverse_discrete_quantile(
+      const Dist& dist,
+      const typename Dist::value_type& p,
+      const typename Dist::value_type& q,
+      const typename Dist::value_type& guess,
+      const typename Dist::value_type& multiplier,
+      const typename Dist::value_type& adder,
+      const policies::discrete_quantile<policies::integer_round_outwards>&,
+      boost::uintmax_t& max_iter)
+{
+   typedef typename Dist::value_type value_type;
+   BOOST_MATH_STD_USING
+   if(p <= pdf(dist, 0))
+      return 0;
+   //
+   // What happens next depends on whether we're looking for an 
+   // upper or lower quantile:
+   //
+   if(p < 0.5f)
+      return floor(do_inverse_discrete_quantile(
+         dist, 
+         p, 
+         q,
+         (guess < 1 ? value_type(1) : (value_type)floor(guess)), 
+         multiplier, 
+         adder, 
+         tools::equal_floor(),
+         max_iter));
+   // else:
+   return ceil(do_inverse_discrete_quantile(
+      dist, 
+      p, 
+      q,
+      (value_type)ceil(guess), 
+      multiplier, 
+      adder, 
+      tools::equal_ceil(),
+      max_iter));
+}
+
+template <class Dist>
+inline typename Dist::value_type 
+   inverse_discrete_quantile(
+      const Dist& dist,
+      const typename Dist::value_type& p,
+      const typename Dist::value_type& q,
+      const typename Dist::value_type& guess,
+      const typename Dist::value_type& multiplier,
+      const typename Dist::value_type& adder,
+      const policies::discrete_quantile<policies::integer_round_inwards>&,
+      boost::uintmax_t& max_iter)
+{
+   typedef typename Dist::value_type value_type;
+   BOOST_MATH_STD_USING
+   if(p <= pdf(dist, 0))
+      return 0;
+   //
+   // What happens next depends on whether we're looking for an 
+   // upper or lower quantile:
+   //
+   if(p < 0.5f)
+      return ceil(do_inverse_discrete_quantile(
+         dist, 
+         p, 
+         q,
+         ceil(guess), 
+         multiplier, 
+         adder, 
+         tools::equal_ceil(),
+         max_iter));
+   // else:
+   return floor(do_inverse_discrete_quantile(
+      dist, 
+      p, 
+      q,
+      (guess < 1 ? value_type(1) : floor(guess)), 
+      multiplier, 
+      adder, 
+      tools::equal_floor(),
+      max_iter));
+}
+
+template <class Dist>
+inline typename Dist::value_type 
+   inverse_discrete_quantile(
+      const Dist& dist,
+      const typename Dist::value_type& p,
+      const typename Dist::value_type& q,
+      const typename Dist::value_type& guess,
+      const typename Dist::value_type& multiplier,
+      const typename Dist::value_type& adder,
+      const policies::discrete_quantile<policies::integer_round_down>&,
+      boost::uintmax_t& max_iter)
+{
+   typedef typename Dist::value_type value_type;
+   BOOST_MATH_STD_USING
+   if(p <= pdf(dist, 0))
+      return 0;
+   return floor(do_inverse_discrete_quantile(
+      dist, 
+      p, 
+      q,
+      (guess < 1 ? value_type(1) : floor(guess)), 
+      multiplier, 
+      adder, 
+      tools::equal_floor(),
+      max_iter));
+}
+
+template <class Dist>
+inline typename Dist::value_type 
+   inverse_discrete_quantile(
+      const Dist& dist,
+      const typename Dist::value_type& p,
+      const typename Dist::value_type& q,
+      const typename Dist::value_type& guess,
+      const typename Dist::value_type& multiplier,
+      const typename Dist::value_type& adder,
+      const policies::discrete_quantile<policies::integer_round_up>&,
+      boost::uintmax_t& max_iter)
+{
+   BOOST_MATH_STD_USING
+   if(p <= pdf(dist, 0))
+      return 0;
+   return ceil(do_inverse_discrete_quantile(
+      dist, 
+      p, 
+      q,
+      ceil(guess), 
+      multiplier, 
+      adder, 
+      tools::equal_ceil(),
+      max_iter));
+}
+
+template <class Dist>
+inline typename Dist::value_type 
+   inverse_discrete_quantile(
+      const Dist& dist,
+      const typename Dist::value_type& p,
+      const typename Dist::value_type& q,
+      const typename Dist::value_type& guess,
+      const typename Dist::value_type& multiplier,
+      const typename Dist::value_type& adder,
+      const policies::discrete_quantile<policies::integer_round_nearest>&,
+      boost::uintmax_t& max_iter)
+{
+   typedef typename Dist::value_type value_type;
+   BOOST_MATH_STD_USING
+   if(p <= pdf(dist, 0))
+      return 0;
+   //
+   // Note that we adjust the guess to the nearest half-integer:
+   // this increase the chances that we will bracket the root
+   // with two results that both round to the same integer quickly.
+   //
+   return floor(do_inverse_discrete_quantile(
+      dist, 
+      p, 
+      q,
+      (guess < 0.5f ? value_type(1.5f) : floor(guess + 0.5f) + 0.5f), 
+      multiplier, 
+      adder, 
+      tools::equal_nearest_integer(),
+      max_iter) + 0.5f);
+}
+
+}}} // namespaces
+
+#endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE
+
+