1 // Copyright 2008 Gautam Sewani
2 // Copyright 2008 John Maddock
4 // Use, modification and distribution are subject to the
5 // Boost Software License, Version 1.0.
6 // (See accompanying file LICENSE_1_0.txt
7 // or copy at http://www.boost.org/LICENSE_1_0.txt)
9 #ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_PDF_HPP
10 #define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_PDF_HPP
12 #include <boost/math/constants/constants.hpp>
13 #include <boost/math/special_functions/lanczos.hpp>
14 #include <boost/math/special_functions/gamma.hpp>
15 #include <boost/math/special_functions/pow.hpp>
16 #include <boost/math/special_functions/prime.hpp>
17 #include <boost/math/policies/error_handling.hpp>
19 #ifdef BOOST_MATH_INSTRUMENT
23 namespace boost{ namespace math{ namespace detail{
25 template <class T, class Func>
26 void bubble_down_one(T* first, T* last, Func f)
31 while((next != last) && (!f(*first, *next)))
42 sort_functor(const T* exponents) : m_exponents(exponents){}
43 bool operator()(int i, int j)
45 return m_exponents[i] > m_exponents[j];
51 template <class T, class Lanczos, class Policy>
52 T hypergeometric_pdf_lanczos_imp(T /*dummy*/, unsigned x, unsigned r, unsigned n, unsigned N, const Lanczos&, const Policy&)
56 BOOST_MATH_INSTRUMENT_FPU
57 BOOST_MATH_INSTRUMENT_VARIABLE(x);
58 BOOST_MATH_INSTRUMENT_VARIABLE(r);
59 BOOST_MATH_INSTRUMENT_VARIABLE(n);
60 BOOST_MATH_INSTRUMENT_VARIABLE(N);
61 BOOST_MATH_INSTRUMENT_VARIABLE(typeid(Lanczos).name());
64 T(n) + static_cast<T>(Lanczos::g()) + 0.5f,
65 T(r) + static_cast<T>(Lanczos::g()) + 0.5f,
66 T(N - n) + static_cast<T>(Lanczos::g()) + 0.5f,
67 T(N - r) + static_cast<T>(Lanczos::g()) + 0.5f,
68 1 / (T(N) + static_cast<T>(Lanczos::g()) + 0.5f),
69 1 / (T(x) + static_cast<T>(Lanczos::g()) + 0.5f),
70 1 / (T(n - x) + static_cast<T>(Lanczos::g()) + 0.5f),
71 1 / (T(r - x) + static_cast<T>(Lanczos::g()) + 0.5f),
72 1 / (T(N - n - r + x) + static_cast<T>(Lanczos::g()) + 0.5f)
83 N - n - r + x + T(0.5f)
85 int base_e_factors[9] = {
86 -1, -1, -1, -1, 1, 1, 1, 1, 1
88 int sorted_indexes[9] = {
89 0, 1, 2, 3, 4, 5, 6, 7, 8
91 #ifdef BOOST_MATH_INSTRUMENT
92 BOOST_MATH_INSTRUMENT_FPU
93 for(unsigned i = 0; i < 9; ++i)
95 BOOST_MATH_INSTRUMENT_VARIABLE(i);
96 BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
97 BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
98 BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
99 BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
102 std::sort(sorted_indexes, sorted_indexes + 9, sort_functor<T>(exponents));
103 #ifdef BOOST_MATH_INSTRUMENT
104 BOOST_MATH_INSTRUMENT_FPU
105 for(unsigned i = 0; i < 9; ++i)
107 BOOST_MATH_INSTRUMENT_VARIABLE(i);
108 BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
109 BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
110 BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
111 BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
116 exponents[sorted_indexes[0]] -= exponents[sorted_indexes[1]];
117 bases[sorted_indexes[1]] *= bases[sorted_indexes[0]];
118 if((bases[sorted_indexes[1]] < tools::min_value<T>()) && (exponents[sorted_indexes[1]] != 0))
122 base_e_factors[sorted_indexes[1]] += base_e_factors[sorted_indexes[0]];
123 bubble_down_one(sorted_indexes, sorted_indexes + 9, sort_functor<T>(exponents));
125 #ifdef BOOST_MATH_INSTRUMENT
126 for(unsigned i = 0; i < 9; ++i)
128 BOOST_MATH_INSTRUMENT_VARIABLE(i);
129 BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
130 BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
131 BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
132 BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
135 }while(exponents[sorted_indexes[1]] > 1);
138 // Combine equal powers:
141 while(exponents[sorted_indexes[j]] == 0) --j;
144 while(j && (exponents[sorted_indexes[j-1]] == exponents[sorted_indexes[j]]))
146 bases[sorted_indexes[j-1]] *= bases[sorted_indexes[j]];
147 exponents[sorted_indexes[j]] = 0;
148 base_e_factors[sorted_indexes[j-1]] += base_e_factors[sorted_indexes[j]];
149 bubble_down_one(sorted_indexes + j, sorted_indexes + 9, sort_functor<T>(exponents));
154 #ifdef BOOST_MATH_INSTRUMENT
155 BOOST_MATH_INSTRUMENT_VARIABLE(j);
156 for(unsigned i = 0; i < 9; ++i)
158 BOOST_MATH_INSTRUMENT_VARIABLE(i);
159 BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
160 BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
161 BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
162 BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
167 #ifdef BOOST_MATH_INSTRUMENT
168 BOOST_MATH_INSTRUMENT_FPU
169 for(unsigned i = 0; i < 9; ++i)
171 BOOST_MATH_INSTRUMENT_VARIABLE(i);
172 BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
173 BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
174 BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
175 BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
180 BOOST_MATH_INSTRUMENT_VARIABLE(bases[sorted_indexes[0]] * exp(static_cast<T>(base_e_factors[sorted_indexes[0]])));
181 BOOST_MATH_INSTRUMENT_VARIABLE(exponents[sorted_indexes[0]]);
183 BOOST_FPU_EXCEPTION_GUARD
184 result = pow(bases[sorted_indexes[0]] * exp(static_cast<T>(base_e_factors[sorted_indexes[0]])), exponents[sorted_indexes[0]]);
186 BOOST_MATH_INSTRUMENT_VARIABLE(result);
187 for(unsigned i = 1; (i < 9) && (exponents[sorted_indexes[i]] > 0); ++i)
189 BOOST_FPU_EXCEPTION_GUARD
190 if(result < tools::min_value<T>())
191 return 0; // short circuit further evaluation
192 if(exponents[sorted_indexes[i]] == 1)
193 result *= bases[sorted_indexes[i]] * exp(static_cast<T>(base_e_factors[sorted_indexes[i]]));
194 else if(exponents[sorted_indexes[i]] == 0.5f)
195 result *= sqrt(bases[sorted_indexes[i]] * exp(static_cast<T>(base_e_factors[sorted_indexes[i]])));
197 result *= pow(bases[sorted_indexes[i]] * exp(static_cast<T>(base_e_factors[sorted_indexes[i]])), exponents[sorted_indexes[i]]);
199 BOOST_MATH_INSTRUMENT_VARIABLE(result);
202 result *= Lanczos::lanczos_sum_expG_scaled(static_cast<T>(n + 1))
203 * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(r + 1))
204 * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N - n + 1))
205 * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N - r + 1))
207 ( Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N + 1))
208 * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(x + 1))
209 * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(n - x + 1))
210 * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(r - x + 1))
211 * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N - n - r + x + 1)));
213 BOOST_MATH_INSTRUMENT_VARIABLE(result);
217 template <class T, class Policy>
218 T hypergeometric_pdf_lanczos_imp(T /*dummy*/, unsigned x, unsigned r, unsigned n, unsigned N, const boost::math::lanczos::undefined_lanczos&, const Policy& pol)
222 boost::math::lgamma(T(n + 1), pol)
223 + boost::math::lgamma(T(r + 1), pol)
224 + boost::math::lgamma(T(N - n + 1), pol)
225 + boost::math::lgamma(T(N - r + 1), pol)
226 - boost::math::lgamma(T(N + 1), pol)
227 - boost::math::lgamma(T(x + 1), pol)
228 - boost::math::lgamma(T(n - x + 1), pol)
229 - boost::math::lgamma(T(r - x + 1), pol)
230 - boost::math::lgamma(T(N - n - r + x + 1), pol));
234 inline T integer_power(const T& x, int ex)
237 return 1 / integer_power(x, -ex);
249 return boost::math::pow<4>(x);
251 return boost::math::pow<5>(x);
253 return boost::math::pow<6>(x);
255 return boost::math::pow<7>(x);
257 return boost::math::pow<8>(x);
261 return pow(x, T(ex));
267 struct hypergeometric_pdf_prime_loop_result_entry
270 const hypergeometric_pdf_prime_loop_result_entry* next;
274 #pragma warning(push)
275 #pragma warning(disable:4510 4512 4610)
278 struct hypergeometric_pdf_prime_loop_data
284 unsigned prime_index;
285 unsigned current_prime;
293 T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hypergeometric_pdf_prime_loop_result_entry<T>& result)
295 while(data.current_prime <= data.N)
297 unsigned base = data.current_prime;
298 int prime_powers = 0;
299 while(base <= data.N)
301 prime_powers += data.n / base;
302 prime_powers += data.r / base;
303 prime_powers += (data.N - data.n) / base;
304 prime_powers += (data.N - data.r) / base;
305 prime_powers -= data.N / base;
306 prime_powers -= data.x / base;
307 prime_powers -= (data.n - data.x) / base;
308 prime_powers -= (data.r - data.x) / base;
309 prime_powers -= (data.N - data.n - data.r + data.x) / base;
310 base *= data.current_prime;
314 T p = integer_power<T>(static_cast<T>(data.current_prime), prime_powers);
315 if((p > 1) && (tools::max_value<T>() / p < result.value))
318 // The next calculation would overflow, use recursion
319 // to sidestep the issue:
321 hypergeometric_pdf_prime_loop_result_entry<T> t = { p, &result };
322 data.current_prime = prime(++data.prime_index);
323 return hypergeometric_pdf_prime_loop_imp<T>(data, t);
325 if((p < 1) && (tools::min_value<T>() / p > result.value))
328 // The next calculation would underflow, use recursion
329 // to sidestep the issue:
331 hypergeometric_pdf_prime_loop_result_entry<T> t = { p, &result };
332 data.current_prime = prime(++data.prime_index);
333 return hypergeometric_pdf_prime_loop_imp<T>(data, t);
337 data.current_prime = prime(++data.prime_index);
340 // When we get to here we have run out of prime factors,
341 // the overall result is the product of all the partial
342 // results we have accumulated on the stack so far, these
343 // are in a linked list starting with "data.head" and ending
346 // All that remains is to multiply them together, taking
347 // care not to overflow or underflow.
349 // Enumerate partial results >= 1 in variable i
350 // and partial results < 1 in variable j:
352 hypergeometric_pdf_prime_loop_result_entry<T> const *i, *j;
354 while(i && i->value < 1)
357 while(j && j->value >= 1)
364 while(i && ((prod <= 1) || (j == 0)))
368 while(i && i->value < 1)
371 while(j && ((prod >= 1) || (i == 0)))
375 while(j && j->value >= 1)
383 template <class T, class Policy>
384 inline T hypergeometric_pdf_prime_imp(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&)
386 hypergeometric_pdf_prime_loop_result_entry<T> result = { 1, 0 };
387 hypergeometric_pdf_prime_loop_data data = { x, r, n, N, 0, prime(0) };
388 return hypergeometric_pdf_prime_loop_imp<T>(data, result);
391 template <class T, class Policy>
392 T hypergeometric_pdf_factorial_imp(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&)
395 BOOST_ASSERT(N <= boost::math::max_factorial<T>::value);
396 T result = boost::math::unchecked_factorial<T>(n);
398 boost::math::unchecked_factorial<T>(r),
399 boost::math::unchecked_factorial<T>(N - n),
400 boost::math::unchecked_factorial<T>(N - r)
403 boost::math::unchecked_factorial<T>(N),
404 boost::math::unchecked_factorial<T>(x),
405 boost::math::unchecked_factorial<T>(n - x),
406 boost::math::unchecked_factorial<T>(r - x),
407 boost::math::unchecked_factorial<T>(N - n - r + x)
411 while((i < 3) || (j < 5))
413 while((j < 5) && ((result >= 1) || (i >= 3)))
418 while((i < 3) && ((result <= 1) || (j >= 5)))
428 template <class T, class Policy>
429 inline typename tools::promote_args<T>::type
430 hypergeometric_pdf(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&)
432 BOOST_FPU_EXCEPTION_GUARD
433 typedef typename tools::promote_args<T>::type result_type;
434 typedef typename policies::evaluation<result_type, Policy>::type value_type;
435 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
436 typedef typename policies::normalise<
438 policies::promote_float<false>,
439 policies::promote_double<false>,
440 policies::discrete_quantile<>,
441 policies::assert_undefined<> >::type forwarding_policy;
444 if(N <= boost::math::max_factorial<value_type>::value)
447 // If N is small enough then we can evaluate the PDF via the factorials
448 // directly: table lookup of the factorials gives the best performance
449 // of the methods available:
451 result = detail::hypergeometric_pdf_factorial_imp<value_type>(x, r, n, N, forwarding_policy());
453 else if(N <= boost::math::prime(boost::math::max_prime - 1))
456 // If N is no larger than the largest prime number in our lookup table
457 // (104729) then we can use prime factorisation to evaluate the PDF,
458 // this is slow but accurate:
460 result = detail::hypergeometric_pdf_prime_imp<value_type>(x, r, n, N, forwarding_policy());
465 // Catch all case - use the lanczos approximation - where available -
466 // to evaluate the ratio of factorials. This is reasonably fast
467 // (almost as quick as using logarithmic evaluation in terms of lgamma)
468 // but only a few digits better in accuracy than using lgamma:
470 result = detail::hypergeometric_pdf_lanczos_imp(value_type(), x, r, n, N, evaluation_type(), forwarding_policy());
482 return policies::checked_narrowing_cast<result_type, forwarding_policy>(result, "boost::math::hypergeometric_pdf<%1%>(%1%,%1%,%1%,%1%)");