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