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

Add some cross-platform Boost headers
author Chris Cannam <cannam@all-day-breakfast.com>
date Sat, 16 Feb 2019 16:31:25 +0000
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/any/include/boost/math/distributions/detail/hypergeometric_cdf.hpp	Sat Feb 16 16:31:25 2019 +0000
@@ -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
+