1
2 // Copyright John Maddock 2006-7, 2013-14.
3 // Copyright Paul A. Bristow 2007, 2013-14.
4 // Copyright Nikhar Agrawal 2013-14
5 // Copyright Christopher Kormanyos 2013-14
6
7 // Use, modification and distribution are subject to the
8 // Boost Software License, Version 1.0. (See accompanying file
9 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
10
11 #ifndef BOOST_MATH_SF_GAMMA_HPP
12 #define BOOST_MATH_SF_GAMMA_HPP
13
14 #ifdef _MSC_VER
15 #pragma once
16 #endif
17
18 #include <boost/config.hpp>
19 #include <boost/math/tools/series.hpp>
20 #include <boost/math/tools/fraction.hpp>
21 #include <boost/math/tools/precision.hpp>
22 #include <boost/math/tools/promotion.hpp>
23 #include <boost/math/policies/error_handling.hpp>
24 #include <boost/math/constants/constants.hpp>
25 #include <boost/math/special_functions/math_fwd.hpp>
26 #include <boost/math/special_functions/log1p.hpp>
27 #include <boost/math/special_functions/trunc.hpp>
28 #include <boost/math/special_functions/powm1.hpp>
29 #include <boost/math/special_functions/sqrt1pm1.hpp>
30 #include <boost/math/special_functions/lanczos.hpp>
31 #include <boost/math/special_functions/fpclassify.hpp>
32 #include <boost/math/special_functions/detail/igamma_large.hpp>
33 #include <boost/math/special_functions/detail/unchecked_factorial.hpp>
34 #include <boost/math/special_functions/detail/lgamma_small.hpp>
35 #include <boost/math/special_functions/bernoulli.hpp>
36 #include <boost/type_traits/is_convertible.hpp>
37 #include <boost/assert.hpp>
38 #include <boost/mpl/greater.hpp>
39 #include <boost/mpl/equal_to.hpp>
40 #include <boost/mpl/greater.hpp>
41
42 #include <boost/config/no_tr1/cmath.hpp>
43 #include <algorithm>
44
45 #ifdef BOOST_MSVC
46 # pragma warning(push)
47 # pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
48 # pragma warning(disable: 4127) // conditional expression is constant.
49 # pragma warning(disable: 4100) // unreferenced formal parameter.
50 // Several variables made comments,
51 // but some difficulty as whether referenced on not may depend on macro values.
52 // So to be safe, 4100 warnings suppressed.
53 // TODO - revisit this?
54 #endif
55
56 namespace boost{ namespace math{
57
58 namespace detail{
59
60 template <class T>
is_odd(T v,const boost::true_type &)61 inline bool is_odd(T v, const boost::true_type&)
62 {
63 int i = static_cast<int>(v);
64 return i&1;
65 }
66 template <class T>
is_odd(T v,const boost::false_type &)67 inline bool is_odd(T v, const boost::false_type&)
68 {
69 // Oh dear can't cast T to int!
70 BOOST_MATH_STD_USING
71 T modulus = v - 2 * floor(v/2);
72 return static_cast<bool>(modulus != 0);
73 }
74 template <class T>
is_odd(T v)75 inline bool is_odd(T v)
76 {
77 return is_odd(v, ::boost::is_convertible<T, int>());
78 }
79
80 template <class T>
sinpx(T z)81 T sinpx(T z)
82 {
83 // Ad hoc function calculates x * sin(pi * x),
84 // taking extra care near when x is near a whole number.
85 BOOST_MATH_STD_USING
86 int sign = 1;
87 if(z < 0)
88 {
89 z = -z;
90 }
91 T fl = floor(z);
92 T dist;
93 if(is_odd(fl))
94 {
95 fl += 1;
96 dist = fl - z;
97 sign = -sign;
98 }
99 else
100 {
101 dist = z - fl;
102 }
103 BOOST_ASSERT(fl >= 0);
104 if(dist > 0.5)
105 dist = 1 - dist;
106 T result = sin(dist*boost::math::constants::pi<T>());
107 return sign*z*result;
108 } // template <class T> T sinpx(T z)
109 //
110 // tgamma(z), with Lanczos support:
111 //
112 template <class T, class Policy, class Lanczos>
113 T gamma_imp(T z, const Policy& pol, const Lanczos& l)
114 {
115 BOOST_MATH_STD_USING
116
117 T result = 1;
118
119 #ifdef BOOST_MATH_INSTRUMENT
120 static bool b = false;
121 if(!b)
122 {
123 std::cout << "tgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
124 b = true;
125 }
126 #endif
127 static const char* function = "boost::math::tgamma<%1%>(%1%)";
128
129 if(z <= 0)
130 {
131 if(floor(z) == z)
132 return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
133 if(z <= -20)
134 {
135 result = gamma_imp(T(-z), pol, l) * sinpx(z);
136 BOOST_MATH_INSTRUMENT_VARIABLE(result);
137 if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
138 return -boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
139 result = -boost::math::constants::pi<T>() / result;
140 if(result == 0)
141 return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
142 if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL)
143 return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);
144 BOOST_MATH_INSTRUMENT_VARIABLE(result);
145 return result;
146 }
147
148 // shift z to > 1:
149 while(z < 0)
150 {
151 result /= z;
152 z += 1;
153 }
154 }
155 BOOST_MATH_INSTRUMENT_VARIABLE(result);
156 if((floor(z) == z) && (z < max_factorial<T>::value))
157 {
158 result *= unchecked_factorial<T>(itrunc(z, pol) - 1);
159 BOOST_MATH_INSTRUMENT_VARIABLE(result);
160 }
161 else if (z < tools::root_epsilon<T>())
162 {
163 if (z < 1 / tools::max_value<T>())
164 result = policies::raise_overflow_error<T>(function, 0, pol);
165 result *= 1 / z - constants::euler<T>();
166 }
167 else
168 {
169 result *= Lanczos::lanczos_sum(z);
170 T zgh = (z + static_cast<T>(Lanczos::g()) - boost::math::constants::half<T>());
171 T lzgh = log(zgh);
172 BOOST_MATH_INSTRUMENT_VARIABLE(result);
173 BOOST_MATH_INSTRUMENT_VARIABLE(tools::log_max_value<T>());
174 if(z * lzgh > tools::log_max_value<T>())
175 {
176 // we're going to overflow unless this is done with care:
177 BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
178 if(lzgh * z / 2 > tools::log_max_value<T>())
179 return boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
180 T hp = pow(zgh, (z / 2) - T(0.25));
181 BOOST_MATH_INSTRUMENT_VARIABLE(hp);
182 result *= hp / exp(zgh);
183 BOOST_MATH_INSTRUMENT_VARIABLE(result);
184 if(tools::max_value<T>() / hp < result)
185 return boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
186 result *= hp;
187 BOOST_MATH_INSTRUMENT_VARIABLE(result);
188 }
189 else
190 {
191 BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
192 BOOST_MATH_INSTRUMENT_VARIABLE(pow(zgh, z - boost::math::constants::half<T>()));
193 BOOST_MATH_INSTRUMENT_VARIABLE(exp(zgh));
194 result *= pow(zgh, z - boost::math::constants::half<T>()) / exp(zgh);
195 BOOST_MATH_INSTRUMENT_VARIABLE(result);
196 }
197 }
198 return result;
199 }
200 //
201 // lgamma(z) with Lanczos support:
202 //
203 template <class T, class Policy, class Lanczos>
204 T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0)
205 {
206 #ifdef BOOST_MATH_INSTRUMENT
207 static bool b = false;
208 if(!b)
209 {
210 std::cout << "lgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
211 b = true;
212 }
213 #endif
214
215 BOOST_MATH_STD_USING
216
217 static const char* function = "boost::math::lgamma<%1%>(%1%)";
218
219 T result = 0;
220 int sresult = 1;
221 if(z <= -tools::root_epsilon<T>())
222 {
223 // reflection formula:
224 if(floor(z) == z)
225 return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol);
226
227 T t = sinpx(z);
228 z = -z;
229 if(t < 0)
230 {
231 t = -t;
232 }
233 else
234 {
235 sresult = -sresult;
236 }
237 result = log(boost::math::constants::pi<T>()) - lgamma_imp(z, pol, l) - log(t);
238 }
239 else if (z < tools::root_epsilon<T>())
240 {
241 if (0 == z)
242 return policies::raise_pole_error<T>(function, "Evaluation of lgamma at %1%.", z, pol);
243 if (fabs(z) < 1 / tools::max_value<T>())
244 result = -log(fabs(z));
245 else
246 result = log(fabs(1 / z - constants::euler<T>()));
247 if (z < 0)
248 sresult = -1;
249 }
250 else if(z < 15)
251 {
252 typedef typename policies::precision<T, Policy>::type precision_type;
253 typedef typename mpl::if_<
254 mpl::and_<
255 mpl::less_equal<precision_type, mpl::int_<64> >,
256 mpl::greater<precision_type, mpl::int_<0> >
257 >,
258 mpl::int_<64>,
259 typename mpl::if_<
260 mpl::and_<
261 mpl::less_equal<precision_type, mpl::int_<113> >,
262 mpl::greater<precision_type, mpl::int_<0> >
263 >,
264 mpl::int_<113>, mpl::int_<0> >::type
265 >::type tag_type;
266 result = lgamma_small_imp<T>(z, T(z - 1), T(z - 2), tag_type(), pol, l);
267 }
268 else if((z >= 3) && (z < 100) && (std::numeric_limits<T>::max_exponent >= 1024))
269 {
270 // taking the log of tgamma reduces the error, no danger of overflow here:
271 result = log(gamma_imp(z, pol, l));
272 }
273 else
274 {
275 // regular evaluation:
276 T zgh = static_cast<T>(z + Lanczos::g() - boost::math::constants::half<T>());
277 result = log(zgh) - 1;
278 result *= z - 0.5f;
279 result += log(Lanczos::lanczos_sum_expG_scaled(z));
280 }
281
282 if(sign)
283 *sign = sresult;
284 return result;
285 }
286
287 //
288 // Incomplete gamma functions follow:
289 //
290 template <class T>
291 struct upper_incomplete_gamma_fract
292 {
293 private:
294 T z, a;
295 int k;
296 public:
297 typedef std::pair<T,T> result_type;
298
upper_incomplete_gamma_fractboost::math::detail::upper_incomplete_gamma_fract299 upper_incomplete_gamma_fract(T a1, T z1)
300 : z(z1-a1+1), a(a1), k(0)
301 {
302 }
303
operator ()boost::math::detail::upper_incomplete_gamma_fract304 result_type operator()()
305 {
306 ++k;
307 z += 2;
308 return result_type(k * (a - k), z);
309 }
310 };
311
312 template <class T>
upper_gamma_fraction(T a,T z,T eps)313 inline T upper_gamma_fraction(T a, T z, T eps)
314 {
315 // Multiply result by z^a * e^-z to get the full
316 // upper incomplete integral. Divide by tgamma(z)
317 // to normalise.
318 upper_incomplete_gamma_fract<T> f(a, z);
319 return 1 / (z - a + 1 + boost::math::tools::continued_fraction_a(f, eps));
320 }
321
322 template <class T>
323 struct lower_incomplete_gamma_series
324 {
325 private:
326 T a, z, result;
327 public:
328 typedef T result_type;
lower_incomplete_gamma_seriesboost::math::detail::lower_incomplete_gamma_series329 lower_incomplete_gamma_series(T a1, T z1) : a(a1), z(z1), result(1){}
330
operator ()boost::math::detail::lower_incomplete_gamma_series331 T operator()()
332 {
333 T r = result;
334 a += 1;
335 result *= z/a;
336 return r;
337 }
338 };
339
340 template <class T, class Policy>
lower_gamma_series(T a,T z,const Policy & pol,T init_value=0)341 inline T lower_gamma_series(T a, T z, const Policy& pol, T init_value = 0)
342 {
343 // Multiply result by ((z^a) * (e^-z) / a) to get the full
344 // lower incomplete integral. Then divide by tgamma(a)
345 // to get the normalised value.
346 lower_incomplete_gamma_series<T> s(a, z);
347 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
348 T factor = policies::get_epsilon<T, Policy>();
349 T result = boost::math::tools::sum_series(s, factor, max_iter, init_value);
350 policies::check_series_iterations<T>("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol);
351 return result;
352 }
353
354 //
355 // Fully generic tgamma and lgamma use Stirling's approximation
356 // with Bernoulli numbers.
357 //
358 template<class T>
highest_bernoulli_index()359 std::size_t highest_bernoulli_index()
360 {
361 const float digits10_of_type = (std::numeric_limits<T>::is_specialized
362 ? static_cast<float>(std::numeric_limits<T>::digits10)
363 : static_cast<float>(boost::math::tools::digits<T>() * 0.301F));
364
365 // Find the high index n for Bn to produce the desired precision in Stirling's calculation.
366 return static_cast<std::size_t>(18.0F + (0.6F * digits10_of_type));
367 }
368
369 template<class T>
minimum_argument_for_bernoulli_recursion()370 T minimum_argument_for_bernoulli_recursion()
371 {
372 const float digits10_of_type = (std::numeric_limits<T>::is_specialized
373 ? static_cast<float>(std::numeric_limits<T>::digits10)
374 : static_cast<float>(boost::math::tools::digits<T>() * 0.301F));
375
376 return T(digits10_of_type * 1.7F);
377 }
378
379 // Forward declaration of the lgamma_imp template specialization.
380 template <class T, class Policy>
381 T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign = 0);
382
383 template <class T, class Policy>
384 T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&)
385 {
386 BOOST_MATH_STD_USING
387
388 static const char* function = "boost::math::tgamma<%1%>(%1%)";
389
390 // Check if the argument of tgamma is identically zero.
391 const bool is_at_zero = (z == 0);
392
393 if(is_at_zero)
394 return policies::raise_domain_error<T>(function, "Evaluation of tgamma at zero %1%.", z, pol);
395
396 const bool b_neg = (z < 0);
397
398 const bool floor_of_z_is_equal_to_z = (floor(z) == z);
399
400 // Special case handling of small factorials:
401 if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial<T>::value))
402 {
403 return boost::math::unchecked_factorial<T>(itrunc(z) - 1);
404 }
405
406 // Make a local, unsigned copy of the input argument.
407 T zz((!b_neg) ? z : -z);
408
409 // Special case for ultra-small z:
410 if(zz < tools::cbrt_epsilon<T>())
411 {
412 const T a0(1);
413 const T a1(boost::math::constants::euler<T>());
414 const T six_euler_squared((boost::math::constants::euler<T>() * boost::math::constants::euler<T>()) * 6);
415 const T a2((six_euler_squared - boost::math::constants::pi_sqr<T>()) / 12);
416
417 const T inverse_tgamma_series = z * ((a2 * z + a1) * z + a0);
418
419 return 1 / inverse_tgamma_series;
420 }
421
422 // Scale the argument up for the calculation of lgamma,
423 // and use downward recursion later for the final result.
424 const T min_arg_for_recursion = minimum_argument_for_bernoulli_recursion<T>();
425
426 int n_recur;
427
428 if(zz < min_arg_for_recursion)
429 {
430 n_recur = boost::math::itrunc(min_arg_for_recursion - zz) + 1;
431
432 zz += n_recur;
433 }
434 else
435 {
436 n_recur = 0;
437 }
438
439 const T log_gamma_value = lgamma_imp(zz, pol, lanczos::undefined_lanczos());
440
441 if(log_gamma_value > tools::log_max_value<T>())
442 return policies::raise_overflow_error<T>(function, 0, pol);
443
444 T gamma_value = exp(log_gamma_value);
445
446 // Rescale the result using downward recursion if necessary.
447 if(n_recur)
448 {
449 // The order of divides is important, if we keep subtracting 1 from zz
450 // we DO NOT get back to z (cancellation error). Further if z < epsilon
451 // we would end up dividing by zero. Also in order to prevent spurious
452 // overflow with the first division, we must save dividing by |z| till last,
453 // so the optimal order of divides is z+1, z+2, z+3...z+n_recur-1,z.
454 zz = fabs(z) + 1;
455 for(int k = 1; k < n_recur; ++k)
456 {
457 gamma_value /= zz;
458 zz += 1;
459 }
460 gamma_value /= fabs(z);
461 }
462
463 // Return the result, accounting for possible negative arguments.
464 if(b_neg)
465 {
466 // Provide special error analysis for:
467 // * arguments in the neighborhood of a negative integer
468 // * arguments exactly equal to a negative integer.
469
470 // Check if the argument of tgamma is exactly equal to a negative integer.
471 if(floor_of_z_is_equal_to_z)
472 return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
473
474 gamma_value *= sinpx(z);
475
476 BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value);
477
478 const bool result_is_too_large_to_represent = ( (abs(gamma_value) < 1)
479 && ((tools::max_value<T>() * abs(gamma_value)) < boost::math::constants::pi<T>()));
480
481 if(result_is_too_large_to_represent)
482 return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
483
484 gamma_value = -boost::math::constants::pi<T>() / gamma_value;
485 BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value);
486
487 if(gamma_value == 0)
488 return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
489
490 if((boost::math::fpclassify)(gamma_value) == static_cast<int>(FP_SUBNORMAL))
491 return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", gamma_value, pol);
492 }
493
494 return gamma_value;
495 }
496
497 template <class T, class Policy>
498 T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign)
499 {
500 BOOST_MATH_STD_USING
501
502 static const char* function = "boost::math::lgamma<%1%>(%1%)";
503
504 // Check if the argument of lgamma is identically zero.
505 const bool is_at_zero = (z == 0);
506
507 if(is_at_zero)
508 return policies::raise_domain_error<T>(function, "Evaluation of lgamma at zero %1%.", z, pol);
509
510 const bool b_neg = (z < 0);
511
512 const bool floor_of_z_is_equal_to_z = (floor(z) == z);
513
514 // Special case handling of small factorials:
515 if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial<T>::value))
516 {
517 return log(boost::math::unchecked_factorial<T>(itrunc(z) - 1));
518 }
519
520 // Make a local, unsigned copy of the input argument.
521 T zz((!b_neg) ? z : -z);
522
523 const T min_arg_for_recursion = minimum_argument_for_bernoulli_recursion<T>();
524
525 T log_gamma_value;
526
527 if (zz < min_arg_for_recursion)
528 {
529 // Here we simply take the logarithm of tgamma(). This is somewhat
530 // inefficient, but simple. The rationale is that the argument here
531 // is relatively small and overflow is not expected to be likely.
532 if (z > -tools::root_epsilon<T>())
533 {
534 // Reflection formula may fail if z is very close to zero, let the series
535 // expansion for tgamma close to zero do the work:
536 log_gamma_value = log(abs(gamma_imp(z, pol, lanczos::undefined_lanczos())));
537 if (sign)
538 {
539 *sign = z < 0 ? -1 : 1;
540 }
541 return log_gamma_value;
542 }
543 else
544 {
545 // No issue with spurious overflow in reflection formula,
546 // just fall through to regular code:
547 log_gamma_value = log(abs(gamma_imp(zz, pol, lanczos::undefined_lanczos())));
548 }
549 }
550 else
551 {
552 // Perform the Bernoulli series expansion of Stirling's approximation.
553
554 const std::size_t number_of_bernoullis_b2n = highest_bernoulli_index<T>();
555
556 T one_over_x_pow_two_n_minus_one = 1 / zz;
557 const T one_over_x2 = one_over_x_pow_two_n_minus_one * one_over_x_pow_two_n_minus_one;
558 T sum = (boost::math::bernoulli_b2n<T>(1) / 2) * one_over_x_pow_two_n_minus_one;
559 const T target_epsilon_to_break_loop = (sum * boost::math::tools::epsilon<T>()) * T(1.0E-10F);
560
561 for(std::size_t n = 2U; n < number_of_bernoullis_b2n; ++n)
562 {
563 one_over_x_pow_two_n_minus_one *= one_over_x2;
564
565 const std::size_t n2 = static_cast<std::size_t>(n * 2U);
566
567 const T term = (boost::math::bernoulli_b2n<T>(static_cast<int>(n)) * one_over_x_pow_two_n_minus_one) / (n2 * (n2 - 1U));
568
569 if((n >= 8U) && (abs(term) < target_epsilon_to_break_loop))
570 {
571 // We have reached the desired precision in Stirling's expansion.
572 // Adding additional terms to the sum of this divergent asymptotic
573 // expansion will not improve the result.
574
575 // Break from the loop.
576 break;
577 }
578
579 sum += term;
580 }
581
582 // Complete Stirling's approximation.
583 const T half_ln_two_pi = log(boost::math::constants::two_pi<T>()) / 2;
584
585 log_gamma_value = ((((zz - boost::math::constants::half<T>()) * log(zz)) - zz) + half_ln_two_pi) + sum;
586 }
587
588 int sign_of_result = 1;
589
590 if(b_neg)
591 {
592 // Provide special error analysis if the argument is exactly
593 // equal to a negative integer.
594
595 // Check if the argument of lgamma is exactly equal to a negative integer.
596 if(floor_of_z_is_equal_to_z)
597 return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol);
598
599 T t = sinpx(z);
600
601 if(t < 0)
602 {
603 t = -t;
604 }
605 else
606 {
607 sign_of_result = -sign_of_result;
608 }
609
610 log_gamma_value = - log_gamma_value
611 + log(boost::math::constants::pi<T>())
612 - log(t);
613 }
614
615 if(sign != static_cast<int*>(0U)) { *sign = sign_of_result; }
616
617 return log_gamma_value;
618 }
619
620 //
621 // This helper calculates tgamma(dz+1)-1 without cancellation errors,
622 // used by the upper incomplete gamma with z < 1:
623 //
624 template <class T, class Policy, class Lanczos>
625 T tgammap1m1_imp(T dz, Policy const& pol, const Lanczos& l)
626 {
627 BOOST_MATH_STD_USING
628
629 typedef typename policies::precision<T,Policy>::type precision_type;
630
631 typedef typename mpl::if_<
632 mpl::or_<
633 mpl::less_equal<precision_type, mpl::int_<0> >,
634 mpl::greater<precision_type, mpl::int_<113> >
635 >,
636 typename mpl::if_<
637 is_same<Lanczos, lanczos::lanczos24m113>,
638 mpl::int_<113>,
639 mpl::int_<0>
640 >::type,
641 typename mpl::if_<
642 mpl::less_equal<precision_type, mpl::int_<64> >,
643 mpl::int_<64>, mpl::int_<113> >::type
644 >::type tag_type;
645
646 T result;
647 if(dz < 0)
648 {
649 if(dz < -0.5)
650 {
651 // Best method is simply to subtract 1 from tgamma:
652 result = boost::math::tgamma(1+dz, pol) - 1;
653 BOOST_MATH_INSTRUMENT_CODE(result);
654 }
655 else
656 {
657 // Use expm1 on lgamma:
658 result = boost::math::expm1(-boost::math::log1p(dz, pol)
659 + lgamma_small_imp<T>(dz+2, dz + 1, dz, tag_type(), pol, l));
660 BOOST_MATH_INSTRUMENT_CODE(result);
661 }
662 }
663 else
664 {
665 if(dz < 2)
666 {
667 // Use expm1 on lgamma:
668 result = boost::math::expm1(lgamma_small_imp<T>(dz+1, dz, dz-1, tag_type(), pol, l), pol);
669 BOOST_MATH_INSTRUMENT_CODE(result);
670 }
671 else
672 {
673 // Best method is simply to subtract 1 from tgamma:
674 result = boost::math::tgamma(1+dz, pol) - 1;
675 BOOST_MATH_INSTRUMENT_CODE(result);
676 }
677 }
678
679 return result;
680 }
681
682 template <class T, class Policy>
tgammap1m1_imp(T dz,Policy const & pol,const::boost::math::lanczos::undefined_lanczos & l)683 inline T tgammap1m1_imp(T dz, Policy const& pol,
684 const ::boost::math::lanczos::undefined_lanczos& l)
685 {
686 BOOST_MATH_STD_USING // ADL of std names
687 //
688 // There should be a better solution than this, but the
689 // algebra isn't easy for the general case....
690 // Start by subracting 1 from tgamma:
691 //
692 T result = gamma_imp(T(1 + dz), pol, l) - 1;
693 BOOST_MATH_INSTRUMENT_CODE(result);
694 //
695 // Test the level of cancellation error observed: we loose one bit
696 // for each power of 2 the result is less than 1. If we would get
697 // more bits from our most precise lgamma rational approximation,
698 // then use that instead:
699 //
700 BOOST_MATH_INSTRUMENT_CODE((dz > -0.5));
701 BOOST_MATH_INSTRUMENT_CODE((dz < 2));
702 BOOST_MATH_INSTRUMENT_CODE((ldexp(1.0, boost::math::policies::digits<T, Policy>()) * fabs(result) < 1e34));
703 if((dz > -0.5) && (dz < 2) && (ldexp(1.0, boost::math::policies::digits<T, Policy>()) * fabs(result) < 1e34))
704 {
705 result = tgammap1m1_imp(dz, pol, boost::math::lanczos::lanczos24m113());
706 BOOST_MATH_INSTRUMENT_CODE(result);
707 }
708 return result;
709 }
710
711 //
712 // Series representation for upper fraction when z is small:
713 //
714 template <class T>
715 struct small_gamma2_series
716 {
717 typedef T result_type;
718
small_gamma2_seriesboost::math::detail::small_gamma2_series719 small_gamma2_series(T a_, T x_) : result(-x_), x(-x_), apn(a_+1), n(1){}
720
operator ()boost::math::detail::small_gamma2_series721 T operator()()
722 {
723 T r = result / (apn);
724 result *= x;
725 result /= ++n;
726 apn += 1;
727 return r;
728 }
729
730 private:
731 T result, x, apn;
732 int n;
733 };
734 //
735 // calculate power term prefix (z^a)(e^-z) used in the non-normalised
736 // incomplete gammas:
737 //
738 template <class T, class Policy>
739 T full_igamma_prefix(T a, T z, const Policy& pol)
740 {
741 BOOST_MATH_STD_USING
742
743 T prefix;
744 T alz = a * log(z);
745
746 if(z >= 1)
747 {
748 if((alz < tools::log_max_value<T>()) && (-z > tools::log_min_value<T>()))
749 {
750 prefix = pow(z, a) * exp(-z);
751 }
752 else if(a >= 1)
753 {
754 prefix = pow(z / exp(z/a), a);
755 }
756 else
757 {
758 prefix = exp(alz - z);
759 }
760 }
761 else
762 {
763 if(alz > tools::log_min_value<T>())
764 {
765 prefix = pow(z, a) * exp(-z);
766 }
767 else if(z/a < tools::log_max_value<T>())
768 {
769 prefix = pow(z / exp(z/a), a);
770 }
771 else
772 {
773 prefix = exp(alz - z);
774 }
775 }
776 //
777 // This error handling isn't very good: it happens after the fact
778 // rather than before it...
779 //
780 if((boost::math::fpclassify)(prefix) == (int)FP_INFINITE)
781 return policies::raise_overflow_error<T>("boost::math::detail::full_igamma_prefix<%1%>(%1%, %1%)", "Result of incomplete gamma function is too large to represent.", pol);
782
783 return prefix;
784 }
785 //
786 // Compute (z^a)(e^-z)/tgamma(a)
787 // most if the error occurs in this function:
788 //
789 template <class T, class Policy, class Lanczos>
790 T regularised_gamma_prefix(T a, T z, const Policy& pol, const Lanczos& l)
791 {
792 BOOST_MATH_STD_USING
793 T agh = a + static_cast<T>(Lanczos::g()) - T(0.5);
794 T prefix;
795 T d = ((z - a) - static_cast<T>(Lanczos::g()) + T(0.5)) / agh;
796
797 if(a < 1)
798 {
799 //
800 // We have to treat a < 1 as a special case because our Lanczos
801 // approximations are optimised against the factorials with a > 1,
802 // and for high precision types especially (128-bit reals for example)
803 // very small values of a can give rather eroneous results for gamma
804 // unless we do this:
805 //
806 // TODO: is this still required? Lanczos approx should be better now?
807 //
808 if(z <= tools::log_min_value<T>())
809 {
810 // Oh dear, have to use logs, should be free of cancellation errors though:
811 return exp(a * log(z) - z - lgamma_imp(a, pol, l));
812 }
813 else
814 {
815 // direct calculation, no danger of overflow as gamma(a) < 1/a
816 // for small a.
817 return pow(z, a) * exp(-z) / gamma_imp(a, pol, l);
818 }
819 }
820 else if((fabs(d*d*a) <= 100) && (a > 150))
821 {
822 // special case for large a and a ~ z.
823 prefix = a * boost::math::log1pmx(d, pol) + z * static_cast<T>(0.5 - Lanczos::g()) / agh;
824 prefix = exp(prefix);
825 }
826 else
827 {
828 //
829 // general case.
830 // direct computation is most accurate, but use various fallbacks
831 // for different parts of the problem domain:
832 //
833 T alz = a * log(z / agh);
834 T amz = a - z;
835 if(((std::min)(alz, amz) <= tools::log_min_value<T>()) || ((std::max)(alz, amz) >= tools::log_max_value<T>()))
836 {
837 T amza = amz / a;
838 if(((std::min)(alz, amz)/2 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/2 < tools::log_max_value<T>()))
839 {
840 // compute square root of the result and then square it:
841 T sq = pow(z / agh, a / 2) * exp(amz / 2);
842 prefix = sq * sq;
843 }
844 else if(((std::min)(alz, amz)/4 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/4 < tools::log_max_value<T>()) && (z > a))
845 {
846 // compute the 4th root of the result then square it twice:
847 T sq = pow(z / agh, a / 4) * exp(amz / 4);
848 prefix = sq * sq;
849 prefix *= prefix;
850 }
851 else if((amza > tools::log_min_value<T>()) && (amza < tools::log_max_value<T>()))
852 {
853 prefix = pow((z * exp(amza)) / agh, a);
854 }
855 else
856 {
857 prefix = exp(alz + amz);
858 }
859 }
860 else
861 {
862 prefix = pow(z / agh, a) * exp(amz);
863 }
864 }
865 prefix *= sqrt(agh / boost::math::constants::e<T>()) / Lanczos::lanczos_sum_expG_scaled(a);
866 return prefix;
867 }
868 //
869 // And again, without Lanczos support:
870 //
871 template <class T, class Policy>
872 T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined_lanczos&)
873 {
874 BOOST_MATH_STD_USING
875
876 T limit = (std::max)(T(10), a);
877 T sum = detail::lower_gamma_series(a, limit, pol) / a;
878 sum += detail::upper_gamma_fraction(a, limit, ::boost::math::policies::get_epsilon<T, Policy>());
879
880 if(a < 10)
881 {
882 // special case for small a:
883 T prefix = pow(z / 10, a);
884 prefix *= exp(10-z);
885 if(0 == prefix)
886 {
887 prefix = pow((z * exp((10-z)/a)) / 10, a);
888 }
889 prefix /= sum;
890 return prefix;
891 }
892
893 T zoa = z / a;
894 T amz = a - z;
895 T alzoa = a * log(zoa);
896 T prefix;
897 if(((std::min)(alzoa, amz) <= tools::log_min_value<T>()) || ((std::max)(alzoa, amz) >= tools::log_max_value<T>()))
898 {
899 T amza = amz / a;
900 if((amza <= tools::log_min_value<T>()) || (amza >= tools::log_max_value<T>()))
901 {
902 prefix = exp(alzoa + amz);
903 }
904 else
905 {
906 prefix = pow(zoa * exp(amza), a);
907 }
908 }
909 else
910 {
911 prefix = pow(zoa, a) * exp(amz);
912 }
913 prefix /= sum;
914 return prefix;
915 }
916 //
917 // Upper gamma fraction for very small a:
918 //
919 template <class T, class Policy>
tgamma_small_upper_part(T a,T x,const Policy & pol,T * pgam=0,bool invert=false,T * pderivative=0)920 inline T tgamma_small_upper_part(T a, T x, const Policy& pol, T* pgam = 0, bool invert = false, T* pderivative = 0)
921 {
922 BOOST_MATH_STD_USING // ADL of std functions.
923 //
924 // Compute the full upper fraction (Q) when a is very small:
925 //
926 T result;
927 result = boost::math::tgamma1pm1(a, pol);
928 if(pgam)
929 *pgam = (result + 1) / a;
930 T p = boost::math::powm1(x, a, pol);
931 result -= p;
932 result /= a;
933 detail::small_gamma2_series<T> s(a, x);
934 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>() - 10;
935 p += 1;
936 if(pderivative)
937 *pderivative = p / (*pgam * exp(x));
938 T init_value = invert ? *pgam : 0;
939 result = -p * tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, (init_value - result) / p);
940 policies::check_series_iterations<T>("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol);
941 if(invert)
942 result = -result;
943 return result;
944 }
945 //
946 // Upper gamma fraction for integer a:
947 //
948 template <class T, class Policy>
finite_gamma_q(T a,T x,Policy const & pol,T * pderivative=0)949 inline T finite_gamma_q(T a, T x, Policy const& pol, T* pderivative = 0)
950 {
951 //
952 // Calculates normalised Q when a is an integer:
953 //
954 BOOST_MATH_STD_USING
955 T e = exp(-x);
956 T sum = e;
957 if(sum != 0)
958 {
959 T term = sum;
960 for(unsigned n = 1; n < a; ++n)
961 {
962 term /= n;
963 term *= x;
964 sum += term;
965 }
966 }
967 if(pderivative)
968 {
969 *pderivative = e * pow(x, a) / boost::math::unchecked_factorial<T>(itrunc(T(a - 1), pol));
970 }
971 return sum;
972 }
973 //
974 // Upper gamma fraction for half integer a:
975 //
976 template <class T, class Policy>
977 T finite_half_gamma_q(T a, T x, T* p_derivative, const Policy& pol)
978 {
979 //
980 // Calculates normalised Q when a is a half-integer:
981 //
982 BOOST_MATH_STD_USING
983 T e = boost::math::erfc(sqrt(x), pol);
984 if((e != 0) && (a > 1))
985 {
986 T term = exp(-x) / sqrt(constants::pi<T>() * x);
987 term *= x;
988 static const T half = T(1) / 2;
989 term /= half;
990 T sum = term;
991 for(unsigned n = 2; n < a; ++n)
992 {
993 term /= n - half;
994 term *= x;
995 sum += term;
996 }
997 e += sum;
998 if(p_derivative)
999 {
1000 *p_derivative = 0;
1001 }
1002 }
1003 else if(p_derivative)
1004 {
1005 // We'll be dividing by x later, so calculate derivative * x:
1006 *p_derivative = sqrt(x) * exp(-x) / constants::root_pi<T>();
1007 }
1008 return e;
1009 }
1010 //
1011 // Main incomplete gamma entry point, handles all four incomplete gamma's:
1012 //
1013 template <class T, class Policy>
1014 T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
1015 const Policy& pol, T* p_derivative)
1016 {
1017 static const char* function = "boost::math::gamma_p<%1%>(%1%, %1%)";
1018 if(a <= 0)
1019 return policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
1020 if(x < 0)
1021 return policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
1022
1023 BOOST_MATH_STD_USING
1024
1025 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1026
1027 T result = 0; // Just to avoid warning C4701: potentially uninitialized local variable 'result' used
1028
1029 if(a >= max_factorial<T>::value && !normalised)
1030 {
1031 //
1032 // When we're computing the non-normalized incomplete gamma
1033 // and a is large the result is rather hard to compute unless
1034 // we use logs. There are really two options - if x is a long
1035 // way from a in value then we can reliably use methods 2 and 4
1036 // below in logarithmic form and go straight to the result.
1037 // Otherwise we let the regularized gamma take the strain
1038 // (the result is unlikely to unerflow in the central region anyway)
1039 // and combine with lgamma in the hopes that we get a finite result.
1040 //
1041 if(invert && (a * 4 < x))
1042 {
1043 // This is method 4 below, done in logs:
1044 result = a * log(x) - x;
1045 if(p_derivative)
1046 *p_derivative = exp(result);
1047 result += log(upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>()));
1048 }
1049 else if(!invert && (a > 4 * x))
1050 {
1051 // This is method 2 below, done in logs:
1052 result = a * log(x) - x;
1053 if(p_derivative)
1054 *p_derivative = exp(result);
1055 T init_value = 0;
1056 result += log(detail::lower_gamma_series(a, x, pol, init_value) / a);
1057 }
1058 else
1059 {
1060 result = gamma_incomplete_imp(a, x, true, invert, pol, p_derivative);
1061 if(result == 0)
1062 {
1063 if(invert)
1064 {
1065 // Try http://functions.wolfram.com/06.06.06.0039.01
1066 result = 1 + 1 / (12 * a) + 1 / (288 * a * a);
1067 result = log(result) - a + (a - 0.5f) * log(a) + log(boost::math::constants::root_two_pi<T>());
1068 if(p_derivative)
1069 *p_derivative = exp(a * log(x) - x);
1070 }
1071 else
1072 {
1073 // This is method 2 below, done in logs, we're really outside the
1074 // range of this method, but since the result is almost certainly
1075 // infinite, we should probably be OK:
1076 result = a * log(x) - x;
1077 if(p_derivative)
1078 *p_derivative = exp(result);
1079 T init_value = 0;
1080 result += log(detail::lower_gamma_series(a, x, pol, init_value) / a);
1081 }
1082 }
1083 else
1084 {
1085 result = log(result) + boost::math::lgamma(a, pol);
1086 }
1087 }
1088 if(result > tools::log_max_value<T>())
1089 return policies::raise_overflow_error<T>(function, 0, pol);
1090 return exp(result);
1091 }
1092
1093 BOOST_ASSERT((p_derivative == 0) || (normalised == true));
1094
1095 bool is_int, is_half_int;
1096 bool is_small_a = (a < 30) && (a <= x + 1) && (x < tools::log_max_value<T>());
1097 if(is_small_a)
1098 {
1099 T fa = floor(a);
1100 is_int = (fa == a);
1101 is_half_int = is_int ? false : (fabs(fa - a) == 0.5f);
1102 }
1103 else
1104 {
1105 is_int = is_half_int = false;
1106 }
1107
1108 int eval_method;
1109
1110 if(is_int && (x > 0.6))
1111 {
1112 // calculate Q via finite sum:
1113 invert = !invert;
1114 eval_method = 0;
1115 }
1116 else if(is_half_int && (x > 0.2))
1117 {
1118 // calculate Q via finite sum for half integer a:
1119 invert = !invert;
1120 eval_method = 1;
1121 }
1122 else if((x < tools::root_epsilon<T>()) && (a > 1))
1123 {
1124 eval_method = 6;
1125 }
1126 else if(x < 0.5)
1127 {
1128 //
1129 // Changeover criterion chosen to give a changeover at Q ~ 0.33
1130 //
1131 if(-0.4 / log(x) < a)
1132 {
1133 eval_method = 2;
1134 }
1135 else
1136 {
1137 eval_method = 3;
1138 }
1139 }
1140 else if(x < 1.1)
1141 {
1142 //
1143 // Changover here occurs when P ~ 0.75 or Q ~ 0.25:
1144 //
1145 if(x * 0.75f < a)
1146 {
1147 eval_method = 2;
1148 }
1149 else
1150 {
1151 eval_method = 3;
1152 }
1153 }
1154 else
1155 {
1156 //
1157 // Begin by testing whether we're in the "bad" zone
1158 // where the result will be near 0.5 and the usual
1159 // series and continued fractions are slow to converge:
1160 //
1161 bool use_temme = false;
1162 if(normalised && std::numeric_limits<T>::is_specialized && (a > 20))
1163 {
1164 T sigma = fabs((x-a)/a);
1165 if((a > 200) && (policies::digits<T, Policy>() <= 113))
1166 {
1167 //
1168 // This limit is chosen so that we use Temme's expansion
1169 // only if the result would be larger than about 10^-6.
1170 // Below that the regular series and continued fractions
1171 // converge OK, and if we use Temme's method we get increasing
1172 // errors from the dominant erfc term as it's (inexact) argument
1173 // increases in magnitude.
1174 //
1175 if(20 / a > sigma * sigma)
1176 use_temme = true;
1177 }
1178 else if(policies::digits<T, Policy>() <= 64)
1179 {
1180 // Note in this zone we can't use Temme's expansion for
1181 // types longer than an 80-bit real:
1182 // it would require too many terms in the polynomials.
1183 if(sigma < 0.4)
1184 use_temme = true;
1185 }
1186 }
1187 if(use_temme)
1188 {
1189 eval_method = 5;
1190 }
1191 else
1192 {
1193 //
1194 // Regular case where the result will not be too close to 0.5.
1195 //
1196 // Changeover here occurs at P ~ Q ~ 0.5
1197 // Note that series computation of P is about x2 faster than continued fraction
1198 // calculation of Q, so try and use the CF only when really necessary, especially
1199 // for small x.
1200 //
1201 if(x - (1 / (3 * x)) < a)
1202 {
1203 eval_method = 2;
1204 }
1205 else
1206 {
1207 eval_method = 4;
1208 invert = !invert;
1209 }
1210 }
1211 }
1212
1213 switch(eval_method)
1214 {
1215 case 0:
1216 {
1217 result = finite_gamma_q(a, x, pol, p_derivative);
1218 if(normalised == false)
1219 result *= boost::math::tgamma(a, pol);
1220 break;
1221 }
1222 case 1:
1223 {
1224 result = finite_half_gamma_q(a, x, p_derivative, pol);
1225 if(normalised == false)
1226 result *= boost::math::tgamma(a, pol);
1227 if(p_derivative && (*p_derivative == 0))
1228 *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
1229 break;
1230 }
1231 case 2:
1232 {
1233 // Compute P:
1234 result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
1235 if(p_derivative)
1236 *p_derivative = result;
1237 if(result != 0)
1238 {
1239 //
1240 // If we're going to be inverting the result then we can
1241 // reduce the number of series evaluations by quite
1242 // a few iterations if we set an initial value for the
1243 // series sum based on what we'll end up subtracting it from
1244 // at the end.
1245 // Have to be careful though that this optimization doesn't
1246 // lead to spurious numberic overflow. Note that the
1247 // scary/expensive overflow checks below are more often
1248 // than not bypassed in practice for "sensible" input
1249 // values:
1250 //
1251 T init_value = 0;
1252 bool optimised_invert = false;
1253 if(invert)
1254 {
1255 init_value = (normalised ? 1 : boost::math::tgamma(a, pol));
1256 if(normalised || (result >= 1) || (tools::max_value<T>() * result > init_value))
1257 {
1258 init_value /= result;
1259 if(normalised || (a < 1) || (tools::max_value<T>() / a > init_value))
1260 {
1261 init_value *= -a;
1262 optimised_invert = true;
1263 }
1264 else
1265 init_value = 0;
1266 }
1267 else
1268 init_value = 0;
1269 }
1270 result *= detail::lower_gamma_series(a, x, pol, init_value) / a;
1271 if(optimised_invert)
1272 {
1273 invert = false;
1274 result = -result;
1275 }
1276 }
1277 break;
1278 }
1279 case 3:
1280 {
1281 // Compute Q:
1282 invert = !invert;
1283 T g;
1284 result = tgamma_small_upper_part(a, x, pol, &g, invert, p_derivative);
1285 invert = false;
1286 if(normalised)
1287 result /= g;
1288 break;
1289 }
1290 case 4:
1291 {
1292 // Compute Q:
1293 result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
1294 if(p_derivative)
1295 *p_derivative = result;
1296 if(result != 0)
1297 result *= upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>());
1298 break;
1299 }
1300 case 5:
1301 {
1302 //
1303 // Use compile time dispatch to the appropriate
1304 // Temme asymptotic expansion. This may be dead code
1305 // if T does not have numeric limits support, or has
1306 // too many digits for the most precise version of
1307 // these expansions, in that case we'll be calling
1308 // an empty function.
1309 //
1310 typedef typename policies::precision<T, Policy>::type precision_type;
1311
1312 typedef typename mpl::if_<
1313 mpl::or_<mpl::equal_to<precision_type, mpl::int_<0> >,
1314 mpl::greater<precision_type, mpl::int_<113> > >,
1315 mpl::int_<0>,
1316 typename mpl::if_<
1317 mpl::less_equal<precision_type, mpl::int_<53> >,
1318 mpl::int_<53>,
1319 typename mpl::if_<
1320 mpl::less_equal<precision_type, mpl::int_<64> >,
1321 mpl::int_<64>,
1322 mpl::int_<113>
1323 >::type
1324 >::type
1325 >::type tag_type;
1326
1327 result = igamma_temme_large(a, x, pol, static_cast<tag_type const*>(0));
1328 if(x >= a)
1329 invert = !invert;
1330 if(p_derivative)
1331 *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
1332 break;
1333 }
1334 case 6:
1335 {
1336 // x is so small that P is necessarily very small too,
1337 // use http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/
1338 result = !normalised ? pow(x, a) / (a) : pow(x, a) / boost::math::tgamma(a + 1, pol);
1339 result *= 1 - a * x / (a + 1);
1340 }
1341 }
1342
1343 if(normalised && (result > 1))
1344 result = 1;
1345 if(invert)
1346 {
1347 T gam = normalised ? 1 : boost::math::tgamma(a, pol);
1348 result = gam - result;
1349 }
1350 if(p_derivative)
1351 {
1352 //
1353 // Need to convert prefix term to derivative:
1354 //
1355 if((x < 1) && (tools::max_value<T>() * x < *p_derivative))
1356 {
1357 // overflow, just return an arbitrarily large value:
1358 *p_derivative = tools::max_value<T>() / 2;
1359 }
1360
1361 *p_derivative /= x;
1362 }
1363
1364 return result;
1365 }
1366
1367 //
1368 // Ratios of two gamma functions:
1369 //
1370 template <class T, class Policy, class Lanczos>
1371 T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos& l)
1372 {
1373 BOOST_MATH_STD_USING
1374 if(z < tools::epsilon<T>())
1375 {
1376 //
1377 // We get spurious numeric overflow unless we're very careful, this
1378 // can occur either inside Lanczos::lanczos_sum(z) or in the
1379 // final combination of terms, to avoid this, split the product up
1380 // into 2 (or 3) parts:
1381 //
1382 // G(z) / G(L) = 1 / (z * G(L)) ; z < eps, L = z + delta = delta
1383 // z * G(L) = z * G(lim) * (G(L)/G(lim)) ; lim = largest factorial
1384 //
1385 if(boost::math::max_factorial<T>::value < delta)
1386 {
1387 T ratio = tgamma_delta_ratio_imp_lanczos(delta, T(boost::math::max_factorial<T>::value - delta), pol, l);
1388 ratio *= z;
1389 ratio *= boost::math::unchecked_factorial<T>(boost::math::max_factorial<T>::value - 1);
1390 return 1 / ratio;
1391 }
1392 else
1393 {
1394 return 1 / (z * boost::math::tgamma(z + delta, pol));
1395 }
1396 }
1397 T zgh = z + Lanczos::g() - constants::half<T>();
1398 T result;
1399 if(fabs(delta) < 10)
1400 {
1401 result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));
1402 }
1403 else
1404 {
1405 result = pow(zgh / (zgh + delta), z - constants::half<T>());
1406 }
1407 // Split the calculation up to avoid spurious overflow:
1408 result *= Lanczos::lanczos_sum(z) / Lanczos::lanczos_sum(T(z + delta));
1409 result *= pow(constants::e<T>() / (zgh + delta), delta);
1410 return result;
1411 }
1412 //
1413 // And again without Lanczos support this time:
1414 //
1415 template <class T, class Policy>
1416 T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const lanczos::undefined_lanczos&)
1417 {
1418 BOOST_MATH_STD_USING
1419 //
1420 // The upper gamma fraction is *very* slow for z < 6, actually it's very
1421 // slow to converge everywhere but recursing until z > 6 gets rid of the
1422 // worst of it's behaviour.
1423 //
1424 T prefix = 1;
1425 T zd = z + delta;
1426 while((zd < 6) && (z < 6))
1427 {
1428 prefix /= z;
1429 prefix *= zd;
1430 z += 1;
1431 zd += 1;
1432 }
1433 if(delta < 10)
1434 {
1435 prefix *= exp(-z * boost::math::log1p(delta / z, pol));
1436 }
1437 else
1438 {
1439 prefix *= pow(z / zd, z);
1440 }
1441 prefix *= pow(constants::e<T>() / zd, delta);
1442 T sum = detail::lower_gamma_series(z, z, pol) / z;
1443 sum += detail::upper_gamma_fraction(z, z, ::boost::math::policies::get_epsilon<T, Policy>());
1444 T sumd = detail::lower_gamma_series(zd, zd, pol) / zd;
1445 sumd += detail::upper_gamma_fraction(zd, zd, ::boost::math::policies::get_epsilon<T, Policy>());
1446 sum /= sumd;
1447 if(fabs(tools::max_value<T>() / prefix) < fabs(sum))
1448 return policies::raise_overflow_error<T>("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Result of tgamma is too large to represent.", pol);
1449 return sum * prefix;
1450 }
1451
1452 template <class T, class Policy>
1453 T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol)
1454 {
1455 BOOST_MATH_STD_USING
1456
1457 if((z <= 0) || (z + delta <= 0))
1458 {
1459 // This isn't very sofisticated, or accurate, but it does work:
1460 return boost::math::tgamma(z, pol) / boost::math::tgamma(z + delta, pol);
1461 }
1462
1463 if(floor(delta) == delta)
1464 {
1465 if(floor(z) == z)
1466 {
1467 //
1468 // Both z and delta are integers, see if we can just use table lookup
1469 // of the factorials to get the result:
1470 //
1471 if((z <= max_factorial<T>::value) && (z + delta <= max_factorial<T>::value))
1472 {
1473 return unchecked_factorial<T>((unsigned)itrunc(z, pol) - 1) / unchecked_factorial<T>((unsigned)itrunc(T(z + delta), pol) - 1);
1474 }
1475 }
1476 if(fabs(delta) < 20)
1477 {
1478 //
1479 // delta is a small integer, we can use a finite product:
1480 //
1481 if(delta == 0)
1482 return 1;
1483 if(delta < 0)
1484 {
1485 z -= 1;
1486 T result = z;
1487 while(0 != (delta += 1))
1488 {
1489 z -= 1;
1490 result *= z;
1491 }
1492 return result;
1493 }
1494 else
1495 {
1496 T result = 1 / z;
1497 while(0 != (delta -= 1))
1498 {
1499 z += 1;
1500 result /= z;
1501 }
1502 return result;
1503 }
1504 }
1505 }
1506 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1507 return tgamma_delta_ratio_imp_lanczos(z, delta, pol, lanczos_type());
1508 }
1509
1510 template <class T, class Policy>
1511 T tgamma_ratio_imp(T x, T y, const Policy& pol)
1512 {
1513 BOOST_MATH_STD_USING
1514
1515 if((x <= 0) || (boost::math::isinf)(x))
1516 return policies::raise_domain_error<T>("boost::math::tgamma_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got a=%1%).", x, pol);
1517 if((y <= 0) || (boost::math::isinf)(y))
1518 return policies::raise_domain_error<T>("boost::math::tgamma_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got b=%1%).", y, pol);
1519
1520 if(x <= tools::min_value<T>())
1521 {
1522 // Special case for denorms...Ugh.
1523 T shift = ldexp(T(1), tools::digits<T>());
1524 return shift * tgamma_ratio_imp(T(x * shift), y, pol);
1525 }
1526
1527 if((x < max_factorial<T>::value) && (y < max_factorial<T>::value))
1528 {
1529 // Rather than subtracting values, lets just call the gamma functions directly:
1530 return boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
1531 }
1532 T prefix = 1;
1533 if(x < 1)
1534 {
1535 if(y < 2 * max_factorial<T>::value)
1536 {
1537 // We need to sidestep on x as well, otherwise we'll underflow
1538 // before we get to factor in the prefix term:
1539 prefix /= x;
1540 x += 1;
1541 while(y >= max_factorial<T>::value)
1542 {
1543 y -= 1;
1544 prefix /= y;
1545 }
1546 return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
1547 }
1548 //
1549 // result is almost certainly going to underflow to zero, try logs just in case:
1550 //
1551 return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
1552 }
1553 if(y < 1)
1554 {
1555 if(x < 2 * max_factorial<T>::value)
1556 {
1557 // We need to sidestep on y as well, otherwise we'll overflow
1558 // before we get to factor in the prefix term:
1559 prefix *= y;
1560 y += 1;
1561 while(x >= max_factorial<T>::value)
1562 {
1563 x -= 1;
1564 prefix *= x;
1565 }
1566 return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
1567 }
1568 //
1569 // Result will almost certainly overflow, try logs just in case:
1570 //
1571 return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
1572 }
1573 //
1574 // Regular case, x and y both large and similar in magnitude:
1575 //
1576 return boost::math::tgamma_delta_ratio(x, y - x, pol);
1577 }
1578
1579 template <class T, class Policy>
1580 T gamma_p_derivative_imp(T a, T x, const Policy& pol)
1581 {
1582 BOOST_MATH_STD_USING
1583 //
1584 // Usual error checks first:
1585 //
1586 if(a <= 0)
1587 return policies::raise_domain_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
1588 if(x < 0)
1589 return policies::raise_domain_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
1590 //
1591 // Now special cases:
1592 //
1593 if(x == 0)
1594 {
1595 return (a > 1) ? 0 :
1596 (a == 1) ? 1 : policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
1597 }
1598 //
1599 // Normal case:
1600 //
1601 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1602 T f1 = detail::regularised_gamma_prefix(a, x, pol, lanczos_type());
1603 if((x < 1) && (tools::max_value<T>() * x < f1))
1604 {
1605 // overflow:
1606 return policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
1607 }
1608 if(f1 == 0)
1609 {
1610 // Underflow in calculation, use logs instead:
1611 f1 = a * log(x) - x - lgamma(a, pol) - log(x);
1612 f1 = exp(f1);
1613 }
1614 else
1615 f1 /= x;
1616
1617 return f1;
1618 }
1619
1620 template <class T, class Policy>
1621 inline typename tools::promote_args<T>::type
tgamma(T z,const Policy &,const mpl::true_)1622 tgamma(T z, const Policy& /* pol */, const mpl::true_)
1623 {
1624 BOOST_FPU_EXCEPTION_GUARD
1625 typedef typename tools::promote_args<T>::type result_type;
1626 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1627 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1628 typedef typename policies::normalise<
1629 Policy,
1630 policies::promote_float<false>,
1631 policies::promote_double<false>,
1632 policies::discrete_quantile<>,
1633 policies::assert_undefined<> >::type forwarding_policy;
1634 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::gamma_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type()), "boost::math::tgamma<%1%>(%1%)");
1635 }
1636
1637 template <class T, class Policy>
1638 struct igamma_initializer
1639 {
1640 struct init
1641 {
initboost::math::detail::igamma_initializer::init1642 init()
1643 {
1644 typedef typename policies::precision<T, Policy>::type precision_type;
1645
1646 typedef typename mpl::if_<
1647 mpl::or_<mpl::equal_to<precision_type, mpl::int_<0> >,
1648 mpl::greater<precision_type, mpl::int_<113> > >,
1649 mpl::int_<0>,
1650 typename mpl::if_<
1651 mpl::less_equal<precision_type, mpl::int_<53> >,
1652 mpl::int_<53>,
1653 typename mpl::if_<
1654 mpl::less_equal<precision_type, mpl::int_<64> >,
1655 mpl::int_<64>,
1656 mpl::int_<113>
1657 >::type
1658 >::type
1659 >::type tag_type;
1660
1661 do_init(tag_type());
1662 }
1663 template <int N>
do_initboost::math::detail::igamma_initializer::init1664 static void do_init(const mpl::int_<N>&)
1665 {
1666 boost::math::gamma_p(static_cast<T>(400), static_cast<T>(400), Policy());
1667 }
do_initboost::math::detail::igamma_initializer::init1668 static void do_init(const mpl::int_<53>&){}
force_instantiateboost::math::detail::igamma_initializer::init1669 void force_instantiate()const{}
1670 };
1671 static const init initializer;
force_instantiateboost::math::detail::igamma_initializer1672 static void force_instantiate()
1673 {
1674 initializer.force_instantiate();
1675 }
1676 };
1677
1678 template <class T, class Policy>
1679 const typename igamma_initializer<T, Policy>::init igamma_initializer<T, Policy>::initializer;
1680
1681 template <class T, class Policy>
1682 struct lgamma_initializer
1683 {
1684 struct init
1685 {
initboost::math::detail::lgamma_initializer::init1686 init()
1687 {
1688 typedef typename policies::precision<T, Policy>::type precision_type;
1689 typedef typename mpl::if_<
1690 mpl::and_<
1691 mpl::less_equal<precision_type, mpl::int_<64> >,
1692 mpl::greater<precision_type, mpl::int_<0> >
1693 >,
1694 mpl::int_<64>,
1695 typename mpl::if_<
1696 mpl::and_<
1697 mpl::less_equal<precision_type, mpl::int_<113> >,
1698 mpl::greater<precision_type, mpl::int_<0> >
1699 >,
1700 mpl::int_<113>, mpl::int_<0> >::type
1701 >::type tag_type;
1702 do_init(tag_type());
1703 }
do_initboost::math::detail::lgamma_initializer::init1704 static void do_init(const mpl::int_<64>&)
1705 {
1706 boost::math::lgamma(static_cast<T>(2.5), Policy());
1707 boost::math::lgamma(static_cast<T>(1.25), Policy());
1708 boost::math::lgamma(static_cast<T>(1.75), Policy());
1709 }
do_initboost::math::detail::lgamma_initializer::init1710 static void do_init(const mpl::int_<113>&)
1711 {
1712 boost::math::lgamma(static_cast<T>(2.5), Policy());
1713 boost::math::lgamma(static_cast<T>(1.25), Policy());
1714 boost::math::lgamma(static_cast<T>(1.5), Policy());
1715 boost::math::lgamma(static_cast<T>(1.75), Policy());
1716 }
do_initboost::math::detail::lgamma_initializer::init1717 static void do_init(const mpl::int_<0>&)
1718 {
1719 }
force_instantiateboost::math::detail::lgamma_initializer::init1720 void force_instantiate()const{}
1721 };
1722 static const init initializer;
force_instantiateboost::math::detail::lgamma_initializer1723 static void force_instantiate()
1724 {
1725 initializer.force_instantiate();
1726 }
1727 };
1728
1729 template <class T, class Policy>
1730 const typename lgamma_initializer<T, Policy>::init lgamma_initializer<T, Policy>::initializer;
1731
1732 template <class T1, class T2, class Policy>
1733 inline typename tools::promote_args<T1, T2>::type
tgamma(T1 a,T2 z,const Policy &,const mpl::false_)1734 tgamma(T1 a, T2 z, const Policy&, const mpl::false_)
1735 {
1736 BOOST_FPU_EXCEPTION_GUARD
1737 typedef typename tools::promote_args<T1, T2>::type result_type;
1738 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1739 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1740 typedef typename policies::normalise<
1741 Policy,
1742 policies::promote_float<false>,
1743 policies::promote_double<false>,
1744 policies::discrete_quantile<>,
1745 policies::assert_undefined<> >::type forwarding_policy;
1746
1747 igamma_initializer<value_type, forwarding_policy>::force_instantiate();
1748
1749 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
1750 detail::gamma_incomplete_imp(static_cast<value_type>(a),
1751 static_cast<value_type>(z), false, true,
1752 forwarding_policy(), static_cast<value_type*>(0)), "boost::math::tgamma<%1%>(%1%, %1%)");
1753 }
1754
1755 template <class T1, class T2>
1756 inline typename tools::promote_args<T1, T2>::type
tgamma(T1 a,T2 z,const mpl::false_ tag)1757 tgamma(T1 a, T2 z, const mpl::false_ tag)
1758 {
1759 return tgamma(a, z, policies::policy<>(), tag);
1760 }
1761
1762
1763 } // namespace detail
1764
1765 template <class T>
1766 inline typename tools::promote_args<T>::type
tgamma(T z)1767 tgamma(T z)
1768 {
1769 return tgamma(z, policies::policy<>());
1770 }
1771
1772 template <class T, class Policy>
1773 inline typename tools::promote_args<T>::type
lgamma(T z,int * sign,const Policy &)1774 lgamma(T z, int* sign, const Policy&)
1775 {
1776 BOOST_FPU_EXCEPTION_GUARD
1777 typedef typename tools::promote_args<T>::type result_type;
1778 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1779 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1780 typedef typename policies::normalise<
1781 Policy,
1782 policies::promote_float<false>,
1783 policies::promote_double<false>,
1784 policies::discrete_quantile<>,
1785 policies::assert_undefined<> >::type forwarding_policy;
1786
1787 detail::lgamma_initializer<value_type, forwarding_policy>::force_instantiate();
1788
1789 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::lgamma_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type(), sign), "boost::math::lgamma<%1%>(%1%)");
1790 }
1791
1792 template <class T>
1793 inline typename tools::promote_args<T>::type
lgamma(T z,int * sign)1794 lgamma(T z, int* sign)
1795 {
1796 return lgamma(z, sign, policies::policy<>());
1797 }
1798
1799 template <class T, class Policy>
1800 inline typename tools::promote_args<T>::type
lgamma(T x,const Policy & pol)1801 lgamma(T x, const Policy& pol)
1802 {
1803 return ::boost::math::lgamma(x, 0, pol);
1804 }
1805
1806 template <class T>
1807 inline typename tools::promote_args<T>::type
lgamma(T x)1808 lgamma(T x)
1809 {
1810 return ::boost::math::lgamma(x, 0, policies::policy<>());
1811 }
1812
1813 template <class T, class Policy>
1814 inline typename tools::promote_args<T>::type
tgamma1pm1(T z,const Policy &)1815 tgamma1pm1(T z, const Policy& /* pol */)
1816 {
1817 BOOST_FPU_EXCEPTION_GUARD
1818 typedef typename tools::promote_args<T>::type result_type;
1819 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1820 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1821 typedef typename policies::normalise<
1822 Policy,
1823 policies::promote_float<false>,
1824 policies::promote_double<false>,
1825 policies::discrete_quantile<>,
1826 policies::assert_undefined<> >::type forwarding_policy;
1827
1828 return policies::checked_narrowing_cast<typename remove_cv<result_type>::type, forwarding_policy>(detail::tgammap1m1_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type()), "boost::math::tgamma1pm1<%!%>(%1%)");
1829 }
1830
1831 template <class T>
1832 inline typename tools::promote_args<T>::type
tgamma1pm1(T z)1833 tgamma1pm1(T z)
1834 {
1835 return tgamma1pm1(z, policies::policy<>());
1836 }
1837
1838 //
1839 // Full upper incomplete gamma:
1840 //
1841 template <class T1, class T2>
1842 inline typename tools::promote_args<T1, T2>::type
tgamma(T1 a,T2 z)1843 tgamma(T1 a, T2 z)
1844 {
1845 //
1846 // Type T2 could be a policy object, or a value, select the
1847 // right overload based on T2:
1848 //
1849 typedef typename policies::is_policy<T2>::type maybe_policy;
1850 return detail::tgamma(a, z, maybe_policy());
1851 }
1852 template <class T1, class T2, class Policy>
1853 inline typename tools::promote_args<T1, T2>::type
tgamma(T1 a,T2 z,const Policy & pol)1854 tgamma(T1 a, T2 z, const Policy& pol)
1855 {
1856 return detail::tgamma(a, z, pol, mpl::false_());
1857 }
1858 //
1859 // Full lower incomplete gamma:
1860 //
1861 template <class T1, class T2, class Policy>
1862 inline typename tools::promote_args<T1, T2>::type
tgamma_lower(T1 a,T2 z,const Policy &)1863 tgamma_lower(T1 a, T2 z, const Policy&)
1864 {
1865 BOOST_FPU_EXCEPTION_GUARD
1866 typedef typename tools::promote_args<T1, T2>::type result_type;
1867 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1868 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1869 typedef typename policies::normalise<
1870 Policy,
1871 policies::promote_float<false>,
1872 policies::promote_double<false>,
1873 policies::discrete_quantile<>,
1874 policies::assert_undefined<> >::type forwarding_policy;
1875
1876 detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
1877
1878 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
1879 detail::gamma_incomplete_imp(static_cast<value_type>(a),
1880 static_cast<value_type>(z), false, false,
1881 forwarding_policy(), static_cast<value_type*>(0)), "tgamma_lower<%1%>(%1%, %1%)");
1882 }
1883 template <class T1, class T2>
1884 inline typename tools::promote_args<T1, T2>::type
tgamma_lower(T1 a,T2 z)1885 tgamma_lower(T1 a, T2 z)
1886 {
1887 return tgamma_lower(a, z, policies::policy<>());
1888 }
1889 //
1890 // Regularised upper incomplete gamma:
1891 //
1892 template <class T1, class T2, class Policy>
1893 inline typename tools::promote_args<T1, T2>::type
gamma_q(T1 a,T2 z,const Policy &)1894 gamma_q(T1 a, T2 z, const Policy& /* pol */)
1895 {
1896 BOOST_FPU_EXCEPTION_GUARD
1897 typedef typename tools::promote_args<T1, T2>::type result_type;
1898 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1899 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1900 typedef typename policies::normalise<
1901 Policy,
1902 policies::promote_float<false>,
1903 policies::promote_double<false>,
1904 policies::discrete_quantile<>,
1905 policies::assert_undefined<> >::type forwarding_policy;
1906
1907 detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
1908
1909 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
1910 detail::gamma_incomplete_imp(static_cast<value_type>(a),
1911 static_cast<value_type>(z), true, true,
1912 forwarding_policy(), static_cast<value_type*>(0)), "gamma_q<%1%>(%1%, %1%)");
1913 }
1914 template <class T1, class T2>
1915 inline typename tools::promote_args<T1, T2>::type
gamma_q(T1 a,T2 z)1916 gamma_q(T1 a, T2 z)
1917 {
1918 return gamma_q(a, z, policies::policy<>());
1919 }
1920 //
1921 // Regularised lower incomplete gamma:
1922 //
1923 template <class T1, class T2, class Policy>
1924 inline typename tools::promote_args<T1, T2>::type
gamma_p(T1 a,T2 z,const Policy &)1925 gamma_p(T1 a, T2 z, const Policy&)
1926 {
1927 BOOST_FPU_EXCEPTION_GUARD
1928 typedef typename tools::promote_args<T1, T2>::type result_type;
1929 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1930 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1931 typedef typename policies::normalise<
1932 Policy,
1933 policies::promote_float<false>,
1934 policies::promote_double<false>,
1935 policies::discrete_quantile<>,
1936 policies::assert_undefined<> >::type forwarding_policy;
1937
1938 detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
1939
1940 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
1941 detail::gamma_incomplete_imp(static_cast<value_type>(a),
1942 static_cast<value_type>(z), true, false,
1943 forwarding_policy(), static_cast<value_type*>(0)), "gamma_p<%1%>(%1%, %1%)");
1944 }
1945 template <class T1, class T2>
1946 inline typename tools::promote_args<T1, T2>::type
gamma_p(T1 a,T2 z)1947 gamma_p(T1 a, T2 z)
1948 {
1949 return gamma_p(a, z, policies::policy<>());
1950 }
1951
1952 // ratios of gamma functions:
1953 template <class T1, class T2, class Policy>
1954 inline typename tools::promote_args<T1, T2>::type
tgamma_delta_ratio(T1 z,T2 delta,const Policy &)1955 tgamma_delta_ratio(T1 z, T2 delta, const Policy& /* pol */)
1956 {
1957 BOOST_FPU_EXCEPTION_GUARD
1958 typedef typename tools::promote_args<T1, T2>::type result_type;
1959 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1960 typedef typename policies::normalise<
1961 Policy,
1962 policies::promote_float<false>,
1963 policies::promote_double<false>,
1964 policies::discrete_quantile<>,
1965 policies::assert_undefined<> >::type forwarding_policy;
1966
1967 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(z), static_cast<value_type>(delta), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
1968 }
1969 template <class T1, class T2>
1970 inline typename tools::promote_args<T1, T2>::type
tgamma_delta_ratio(T1 z,T2 delta)1971 tgamma_delta_ratio(T1 z, T2 delta)
1972 {
1973 return tgamma_delta_ratio(z, delta, policies::policy<>());
1974 }
1975 template <class T1, class T2, class Policy>
1976 inline typename tools::promote_args<T1, T2>::type
tgamma_ratio(T1 a,T2 b,const Policy &)1977 tgamma_ratio(T1 a, T2 b, const Policy&)
1978 {
1979 typedef typename tools::promote_args<T1, T2>::type result_type;
1980 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1981 typedef typename policies::normalise<
1982 Policy,
1983 policies::promote_float<false>,
1984 policies::promote_double<false>,
1985 policies::discrete_quantile<>,
1986 policies::assert_undefined<> >::type forwarding_policy;
1987
1988 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(b), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
1989 }
1990 template <class T1, class T2>
1991 inline typename tools::promote_args<T1, T2>::type
tgamma_ratio(T1 a,T2 b)1992 tgamma_ratio(T1 a, T2 b)
1993 {
1994 return tgamma_ratio(a, b, policies::policy<>());
1995 }
1996
1997 template <class T1, class T2, class Policy>
1998 inline typename tools::promote_args<T1, T2>::type
gamma_p_derivative(T1 a,T2 x,const Policy &)1999 gamma_p_derivative(T1 a, T2 x, const Policy&)
2000 {
2001 BOOST_FPU_EXCEPTION_GUARD
2002 typedef typename tools::promote_args<T1, T2>::type result_type;
2003 typedef typename policies::evaluation<result_type, Policy>::type value_type;
2004 typedef typename policies::normalise<
2005 Policy,
2006 policies::promote_float<false>,
2007 policies::promote_double<false>,
2008 policies::discrete_quantile<>,
2009 policies::assert_undefined<> >::type forwarding_policy;
2010
2011 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::gamma_p_derivative_imp(static_cast<value_type>(a), static_cast<value_type>(x), forwarding_policy()), "boost::math::gamma_p_derivative<%1%>(%1%, %1%)");
2012 }
2013 template <class T1, class T2>
2014 inline typename tools::promote_args<T1, T2>::type
gamma_p_derivative(T1 a,T2 x)2015 gamma_p_derivative(T1 a, T2 x)
2016 {
2017 return gamma_p_derivative(a, x, policies::policy<>());
2018 }
2019
2020 } // namespace math
2021 } // namespace boost
2022
2023 #ifdef BOOST_MSVC
2024 # pragma warning(pop)
2025 #endif
2026
2027 #include <boost/math/special_functions/detail/igamma_inverse.hpp>
2028 #include <boost/math/special_functions/detail/gamma_inva.hpp>
2029 #include <boost/math/special_functions/erf.hpp>
2030
2031 #endif // BOOST_MATH_SF_GAMMA_HPP
2032
2033
2034
2035
2036