1 | // Copyright 2008 Gautam Sewani
|
---|
2 | // Copyright 2008 John Maddock
|
---|
3 | //
|
---|
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)
|
---|
8 |
|
---|
9 | #ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_PDF_HPP
|
---|
10 | #define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_PDF_HPP
|
---|
11 |
|
---|
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>
|
---|
18 |
|
---|
19 | #ifdef BOOST_MATH_INSTRUMENT
|
---|
20 | #include <typeinfo>
|
---|
21 | #endif
|
---|
22 |
|
---|
23 | namespace boost{ namespace math{ namespace detail{
|
---|
24 |
|
---|
25 | template <class T, class Func>
|
---|
26 | void bubble_down_one(T* first, T* last, Func f)
|
---|
27 | {
|
---|
28 | using std::swap;
|
---|
29 | T* next = first;
|
---|
30 | ++next;
|
---|
31 | while((next != last) && (!f(*first, *next)))
|
---|
32 | {
|
---|
33 | swap(*first, *next);
|
---|
34 | ++first;
|
---|
35 | ++next;
|
---|
36 | }
|
---|
37 | }
|
---|
38 |
|
---|
39 | template <class T>
|
---|
40 | struct sort_functor
|
---|
41 | {
|
---|
42 | sort_functor(const T* exponents) : m_exponents(exponents){}
|
---|
43 | bool operator()(int i, int j)
|
---|
44 | {
|
---|
45 | return m_exponents[i] > m_exponents[j];
|
---|
46 | }
|
---|
47 | private:
|
---|
48 | const T* m_exponents;
|
---|
49 | };
|
---|
50 |
|
---|
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&)
|
---|
53 | {
|
---|
54 | BOOST_MATH_STD_USING
|
---|
55 |
|
---|
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());
|
---|
62 |
|
---|
63 | T bases[9] = {
|
---|
64 | T(n) + Lanczos::g() + 0.5f,
|
---|
65 | T(r) + Lanczos::g() + 0.5f,
|
---|
66 | T(N - n) + Lanczos::g() + 0.5f,
|
---|
67 | T(N - r) + Lanczos::g() + 0.5f,
|
---|
68 | 1 / (T(N) + Lanczos::g() + 0.5f),
|
---|
69 | 1 / (T(x) + Lanczos::g() + 0.5f),
|
---|
70 | 1 / (T(n - x) + Lanczos::g() + 0.5f),
|
---|
71 | 1 / (T(r - x) + Lanczos::g() + 0.5f),
|
---|
72 | 1 / (T(N - n - r + x) + Lanczos::g() + 0.5f)
|
---|
73 | };
|
---|
74 | T exponents[9] = {
|
---|
75 | n + 0.5f,
|
---|
76 | r + 0.5f,
|
---|
77 | N - n + 0.5f,
|
---|
78 | N - r + 0.5f,
|
---|
79 | N + 0.5f,
|
---|
80 | x + 0.5f,
|
---|
81 | n - x + 0.5f,
|
---|
82 | r - x + 0.5f,
|
---|
83 | N - n - r + x + 0.5f
|
---|
84 | };
|
---|
85 | int base_e_factors[9] = {
|
---|
86 | -1, -1, -1, -1, 1, 1, 1, 1, 1
|
---|
87 | };
|
---|
88 | int sorted_indexes[9] = {
|
---|
89 | 0, 1, 2, 3, 4, 5, 6, 7, 8
|
---|
90 | };
|
---|
91 | #ifdef BOOST_MATH_INSTRUMENT
|
---|
92 | BOOST_MATH_INSTRUMENT_FPU
|
---|
93 | for(unsigned i = 0; i < 9; ++i)
|
---|
94 | {
|
---|
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]);
|
---|
100 | }
|
---|
101 | #endif
|
---|
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)
|
---|
106 | {
|
---|
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]);
|
---|
112 | }
|
---|
113 | #endif
|
---|
114 |
|
---|
115 | do{
|
---|
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))
|
---|
119 | {
|
---|
120 | return 0;
|
---|
121 | }
|
---|
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));
|
---|
124 |
|
---|
125 | #ifdef BOOST_MATH_INSTRUMENT
|
---|
126 | for(unsigned i = 0; i < 9; ++i)
|
---|
127 | {
|
---|
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]);
|
---|
133 | }
|
---|
134 | #endif
|
---|
135 | }while(exponents[sorted_indexes[1]] > 1);
|
---|
136 |
|
---|
137 | //
|
---|
138 | // Combine equal powers:
|
---|
139 | //
|
---|
140 | int j = 8;
|
---|
141 | while(exponents[sorted_indexes[j]] == 0) --j;
|
---|
142 | while(j)
|
---|
143 | {
|
---|
144 | while(j && (exponents[sorted_indexes[j-1]] == exponents[sorted_indexes[j]]))
|
---|
145 | {
|
---|
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));
|
---|
150 | --j;
|
---|
151 | }
|
---|
152 | --j;
|
---|
153 |
|
---|
154 | #ifdef BOOST_MATH_INSTRUMENT
|
---|
155 | BOOST_MATH_INSTRUMENT_VARIABLE(j);
|
---|
156 | for(unsigned i = 0; i < 9; ++i)
|
---|
157 | {
|
---|
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]);
|
---|
163 | }
|
---|
164 | #endif
|
---|
165 | }
|
---|
166 |
|
---|
167 | #ifdef BOOST_MATH_INSTRUMENT
|
---|
168 | BOOST_MATH_INSTRUMENT_FPU
|
---|
169 | for(unsigned i = 0; i < 9; ++i)
|
---|
170 | {
|
---|
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]);
|
---|
176 | }
|
---|
177 | #endif
|
---|
178 |
|
---|
179 | T result;
|
---|
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]]);
|
---|
182 | {
|
---|
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]]);
|
---|
185 | }
|
---|
186 | BOOST_MATH_INSTRUMENT_VARIABLE(result);
|
---|
187 | for(unsigned i = 1; (i < 9) && (exponents[sorted_indexes[i]] > 0); ++i)
|
---|
188 | {
|
---|
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]])));
|
---|
196 | else
|
---|
197 | result *= pow(bases[sorted_indexes[i]] * exp(static_cast<T>(base_e_factors[sorted_indexes[i]])), exponents[sorted_indexes[i]]);
|
---|
198 |
|
---|
199 | BOOST_MATH_INSTRUMENT_VARIABLE(result);
|
---|
200 | }
|
---|
201 |
|
---|
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))
|
---|
206 | /
|
---|
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)));
|
---|
212 |
|
---|
213 | BOOST_MATH_INSTRUMENT_VARIABLE(result);
|
---|
214 | return result;
|
---|
215 | }
|
---|
216 |
|
---|
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)
|
---|
219 | {
|
---|
220 | BOOST_MATH_STD_USING
|
---|
221 | return exp(
|
---|
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));
|
---|
231 | }
|
---|
232 |
|
---|
233 | template <class T>
|
---|
234 | inline T integer_power(const T& x, int ex)
|
---|
235 | {
|
---|
236 | if(ex < 0)
|
---|
237 | return 1 / integer_power(x, -ex);
|
---|
238 | switch(ex)
|
---|
239 | {
|
---|
240 | case 0:
|
---|
241 | return 1;
|
---|
242 | case 1:
|
---|
243 | return x;
|
---|
244 | case 2:
|
---|
245 | return x * x;
|
---|
246 | case 3:
|
---|
247 | return x * x * x;
|
---|
248 | case 4:
|
---|
249 | return boost::math::pow<4>(x);
|
---|
250 | case 5:
|
---|
251 | return boost::math::pow<5>(x);
|
---|
252 | case 6:
|
---|
253 | return boost::math::pow<6>(x);
|
---|
254 | case 7:
|
---|
255 | return boost::math::pow<7>(x);
|
---|
256 | case 8:
|
---|
257 | return boost::math::pow<8>(x);
|
---|
258 | }
|
---|
259 | BOOST_MATH_STD_USING
|
---|
260 | #ifdef __SUNPRO_CC
|
---|
261 | return pow(x, T(ex));
|
---|
262 | #else
|
---|
263 | return pow(x, ex);
|
---|
264 | #endif
|
---|
265 | }
|
---|
266 | template <class T>
|
---|
267 | struct hypergeometric_pdf_prime_loop_result_entry
|
---|
268 | {
|
---|
269 | T value;
|
---|
270 | const hypergeometric_pdf_prime_loop_result_entry* next;
|
---|
271 | };
|
---|
272 |
|
---|
273 | #ifdef BOOST_MSVC
|
---|
274 | #pragma warning(push)
|
---|
275 | #pragma warning(disable:4510 4512 4610)
|
---|
276 | #endif
|
---|
277 |
|
---|
278 | struct hypergeometric_pdf_prime_loop_data
|
---|
279 | {
|
---|
280 | const unsigned x;
|
---|
281 | const unsigned r;
|
---|
282 | const unsigned n;
|
---|
283 | const unsigned N;
|
---|
284 | unsigned prime_index;
|
---|
285 | unsigned current_prime;
|
---|
286 | };
|
---|
287 |
|
---|
288 | #ifdef BOOST_MSVC
|
---|
289 | #pragma warning(pop)
|
---|
290 | #endif
|
---|
291 |
|
---|
292 | template <class T>
|
---|
293 | T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hypergeometric_pdf_prime_loop_result_entry<T>& result)
|
---|
294 | {
|
---|
295 | while(data.current_prime <= data.N)
|
---|
296 | {
|
---|
297 | unsigned base = data.current_prime;
|
---|
298 | int prime_powers = 0;
|
---|
299 | while(base <= data.N)
|
---|
300 | {
|
---|
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;
|
---|
311 | }
|
---|
312 | if(prime_powers)
|
---|
313 | {
|
---|
314 | T p = integer_power<T>(data.current_prime, prime_powers);
|
---|
315 | if((p > 1) && (tools::max_value<T>() / p < result.value))
|
---|
316 | {
|
---|
317 | //
|
---|
318 | // The next calculation would overflow, use recursion
|
---|
319 | // to sidestep the issue:
|
---|
320 | //
|
---|
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);
|
---|
324 | }
|
---|
325 | if((p < 1) && (tools::min_value<T>() / p > result.value))
|
---|
326 | {
|
---|
327 | //
|
---|
328 | // The next calculation would underflow, use recursion
|
---|
329 | // to sidestep the issue:
|
---|
330 | //
|
---|
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);
|
---|
334 | }
|
---|
335 | result.value *= p;
|
---|
336 | }
|
---|
337 | data.current_prime = prime(++data.prime_index);
|
---|
338 | }
|
---|
339 | //
|
---|
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
|
---|
344 | // with "result".
|
---|
345 | //
|
---|
346 | // All that remains is to multiply them together, taking
|
---|
347 | // care not to overflow or underflow.
|
---|
348 | //
|
---|
349 | // Enumerate partial results >= 1 in variable i
|
---|
350 | // and partial results < 1 in variable j:
|
---|
351 | //
|
---|
352 | hypergeometric_pdf_prime_loop_result_entry<T> const *i, *j;
|
---|
353 | i = &result;
|
---|
354 | while(i && i->value < 1)
|
---|
355 | i = i->next;
|
---|
356 | j = &result;
|
---|
357 | while(j && j->value >= 1)
|
---|
358 | j = j->next;
|
---|
359 |
|
---|
360 | T prod = 1;
|
---|
361 |
|
---|
362 | while(i || j)
|
---|
363 | {
|
---|
364 | while(i && ((prod <= 1) || (j == 0)))
|
---|
365 | {
|
---|
366 | prod *= i->value;
|
---|
367 | i = i->next;
|
---|
368 | while(i && i->value < 1)
|
---|
369 | i = i->next;
|
---|
370 | }
|
---|
371 | while(j && ((prod >= 1) || (i == 0)))
|
---|
372 | {
|
---|
373 | prod *= j->value;
|
---|
374 | j = j->next;
|
---|
375 | while(j && j->value >= 1)
|
---|
376 | j = j->next;
|
---|
377 | }
|
---|
378 | }
|
---|
379 |
|
---|
380 | return prod;
|
---|
381 | }
|
---|
382 |
|
---|
383 | template <class T, class Policy>
|
---|
384 | inline T hypergeometric_pdf_prime_imp(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&)
|
---|
385 | {
|
---|
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);
|
---|
389 | }
|
---|
390 |
|
---|
391 | template <class T, class Policy>
|
---|
392 | T hypergeometric_pdf_factorial_imp(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&)
|
---|
393 | {
|
---|
394 | BOOST_MATH_STD_USING
|
---|
395 | BOOST_ASSERT(N < boost::math::max_factorial<T>::value);
|
---|
396 | T result = boost::math::unchecked_factorial<T>(n);
|
---|
397 | T num[3] = {
|
---|
398 | boost::math::unchecked_factorial<T>(r),
|
---|
399 | boost::math::unchecked_factorial<T>(N - n),
|
---|
400 | boost::math::unchecked_factorial<T>(N - r)
|
---|
401 | };
|
---|
402 | T denom[5] = {
|
---|
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)
|
---|
408 | };
|
---|
409 | int i = 0;
|
---|
410 | int j = 0;
|
---|
411 | while((i < 3) || (j < 5))
|
---|
412 | {
|
---|
413 | while((j < 5) && ((result >= 1) || (i >= 3)))
|
---|
414 | {
|
---|
415 | result /= denom[j];
|
---|
416 | ++j;
|
---|
417 | }
|
---|
418 | while((i < 3) && ((result <= 1) || (j >= 5)))
|
---|
419 | {
|
---|
420 | result *= num[i];
|
---|
421 | ++i;
|
---|
422 | }
|
---|
423 | }
|
---|
424 | return result;
|
---|
425 | }
|
---|
426 |
|
---|
427 |
|
---|
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&)
|
---|
431 | {
|
---|
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<
|
---|
437 | Policy,
|
---|
438 | policies::promote_float<false>,
|
---|
439 | policies::promote_double<false>,
|
---|
440 | policies::discrete_quantile<>,
|
---|
441 | policies::assert_undefined<> >::type forwarding_policy;
|
---|
442 |
|
---|
443 | value_type result;
|
---|
444 | if(N <= boost::math::max_factorial<value_type>::value)
|
---|
445 | {
|
---|
446 | //
|
---|
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:
|
---|
450 | //
|
---|
451 | result = detail::hypergeometric_pdf_factorial_imp<value_type>(x, r, n, N, forwarding_policy());
|
---|
452 | }
|
---|
453 | else if(N <= boost::math::prime(boost::math::max_prime - 1))
|
---|
454 | {
|
---|
455 | //
|
---|
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:
|
---|
459 | //
|
---|
460 | result = detail::hypergeometric_pdf_prime_imp<value_type>(x, r, n, N, forwarding_policy());
|
---|
461 | }
|
---|
462 | else
|
---|
463 | {
|
---|
464 | //
|
---|
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:
|
---|
469 | //
|
---|
470 | result = detail::hypergeometric_pdf_lanczos_imp(value_type(), x, r, n, N, evaluation_type(), forwarding_policy());
|
---|
471 | }
|
---|
472 |
|
---|
473 | if(result > 1)
|
---|
474 | {
|
---|
475 | result = 1;
|
---|
476 | }
|
---|
477 | if(result < 0)
|
---|
478 | {
|
---|
479 | result = 0;
|
---|
480 | }
|
---|
481 |
|
---|
482 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(result, "boost::math::hypergeometric_pdf<%1%>(%1%,%1%,%1%,%1%)");
|
---|
483 | }
|
---|
484 |
|
---|
485 | }}} // namespaces
|
---|
486 |
|
---|
487 | #endif
|
---|
488 |
|
---|