1// (C) Copyright John Maddock 2006.
2// Use, modification and distribution are subject to the
3// Boost Software License, Version 1.0. (See accompanying file
4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5
6#ifndef BOOST_MATH_SPECIAL_FUNCTIONS_IGAMMA_INVERSE_HPP
7#define BOOST_MATH_SPECIAL_FUNCTIONS_IGAMMA_INVERSE_HPP
8
9#ifdef _MSC_VER
10#pragma once
11#endif
12
13#include <boost/math/tools/tuple.hpp>
14#include <boost/math/special_functions/gamma.hpp>
15#include <boost/math/special_functions/sign.hpp>
16#include <boost/math/tools/roots.hpp>
17#include <boost/math/policies/error_handling.hpp>
18
19namespace boost{ namespace math{
20
21namespace detail{
22
23template <class T>
24T find_inverse_s(T p, T q)
25{
26 //
27 // Computation of the Incomplete Gamma Function Ratios and their Inverse
28 // ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR.
29 // ACM Transactions on Mathematical Software, Vol. 12, No. 4,
30 // December 1986, Pages 377-393.
31 //
32 // See equation 32.
33 //
34 BOOST_MATH_STD_USING
35 T t;
36 if(p < T(0.5))
37 {
38 t = sqrt(-2 * log(p));
39 }
40 else
41 {
42 t = sqrt(-2 * log(q));
43 }
44 static const double a[4] = { 3.31125922108741, 11.6616720288968, 4.28342155967104, 0.213623493715853 };
45 static const double b[5] = { 1, 6.61053765625462, 6.40691597760039, 1.27364489782223, 0.3611708101884203e-1 };
46 T s = t - tools::evaluate_polynomial(a, t) / tools::evaluate_polynomial(b, t);
47 if(p < T(0.5))
48 s = -s;
49 return s;
50}
51
52template <class T>
53T didonato_SN(T a, T x, unsigned N, T tolerance = 0)
54{
55 //
56 // Computation of the Incomplete Gamma Function Ratios and their Inverse
57 // ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR.
58 // ACM Transactions on Mathematical Software, Vol. 12, No. 4,
59 // December 1986, Pages 377-393.
60 //
61 // See equation 34.
62 //
63 T sum = 1;
64 if(N >= 1)
65 {
66 T partial = x / (a + 1);
67 sum += partial;
68 for(unsigned i = 2; i <= N; ++i)
69 {
70 partial *= x / (a + i);
71 sum += partial;
72 if(partial < tolerance)
73 break;
74 }
75 }
76 return sum;
77}
78
79template <class T, class Policy>
80inline T didonato_FN(T p, T a, T x, unsigned N, T tolerance, const Policy& pol)
81{
82 //
83 // Computation of the Incomplete Gamma Function Ratios and their Inverse
84 // ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR.
85 // ACM Transactions on Mathematical Software, Vol. 12, No. 4,
86 // December 1986, Pages 377-393.
87 //
88 // See equation 34.
89 //
90 BOOST_MATH_STD_USING
91 T u = log(p) + boost::math::lgamma(a + 1, pol);
92 return exp((u + x - log(didonato_SN(a, x, N, tolerance))) / a);
93}
94
95template <class T, class Policy>
96T find_inverse_gamma(T a, T p, T q, const Policy& pol, bool* p_has_10_digits)
97{
98 //
99 // In order to understand what's going on here, you will
100 // need to refer to:
101 //
102 // Computation of the Incomplete Gamma Function Ratios and their Inverse
103 // ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR.
104 // ACM Transactions on Mathematical Software, Vol. 12, No. 4,
105 // December 1986, Pages 377-393.
106 //
107 BOOST_MATH_STD_USING
108
109 T result;
110 *p_has_10_digits = false;
111
112 if(a == 1)
113 {
114 result = -log(q);
115 BOOST_MATH_INSTRUMENT_VARIABLE(result);
116 }
117 else if(a < 1)
118 {
119 T g = boost::math::tgamma(a, pol);
120 T b = q * g;
121 BOOST_MATH_INSTRUMENT_VARIABLE(g);
122 BOOST_MATH_INSTRUMENT_VARIABLE(b);
123 if((b >T(0.6)) || ((b >= T(0.45)) && (a >= T(0.3))))
124 {
125 // DiDonato & Morris Eq 21:
126 //
127 // There is a slight variation from DiDonato and Morris here:
128 // the first form given here is unstable when p is close to 1,
129 // making it impossible to compute the inverse of Q(a,x) for small
130 // q. Fortunately the second form works perfectly well in this case.
131 //
132 T u;
133 if((b * q > T(1e-8)) && (q > T(1e-5)))
134 {
135 u = pow(p * g * a, 1 / a);
136 BOOST_MATH_INSTRUMENT_VARIABLE(u);
137 }
138 else
139 {
140 u = exp((-q / a) - constants::euler<T>());
141 BOOST_MATH_INSTRUMENT_VARIABLE(u);
142 }
143 result = u / (1 - (u / (a + 1)));
144 BOOST_MATH_INSTRUMENT_VARIABLE(result);
145 }
146 else if((a < 0.3) && (b >= 0.35))
147 {
148 // DiDonato & Morris Eq 22:
149 T t = exp(-constants::euler<T>() - b);
150 T u = t * exp(t);
151 result = t * exp(u);
152 BOOST_MATH_INSTRUMENT_VARIABLE(result);
153 }
154 else if((b > 0.15) || (a >= 0.3))
155 {
156 // DiDonato & Morris Eq 23:
157 T y = -log(b);
158 T u = y - (1 - a) * log(y);
159 result = y - (1 - a) * log(u) - log(1 + (1 - a) / (1 + u));
160 BOOST_MATH_INSTRUMENT_VARIABLE(result);
161 }
162 else if (b > 0.1)
163 {
164 // DiDonato & Morris Eq 24:
165 T y = -log(b);
166 T u = y - (1 - a) * log(y);
167 result = y - (1 - a) * log(u) - log((u * u + 2 * (3 - a) * u + (2 - a) * (3 - a)) / (u * u + (5 - a) * u + 2));
168 BOOST_MATH_INSTRUMENT_VARIABLE(result);
169 }
170 else
171 {
172 // DiDonato & Morris Eq 25:
173 T y = -log(b);
174 T c1 = (a - 1) * log(y);
175 T c1_2 = c1 * c1;
176 T c1_3 = c1_2 * c1;
177 T c1_4 = c1_2 * c1_2;
178 T a_2 = a * a;
179 T a_3 = a_2 * a;
180
181 T c2 = (a - 1) * (1 + c1);
182 T c3 = (a - 1) * (-(c1_2 / 2) + (a - 2) * c1 + (3 * a - 5) / 2);
183 T c4 = (a - 1) * ((c1_3 / 3) - (3 * a - 5) * c1_2 / 2 + (a_2 - 6 * a + 7) * c1 + (11 * a_2 - 46 * a + 47) / 6);
184 T c5 = (a - 1) * (-(c1_4 / 4)
185 + (11 * a - 17) * c1_3 / 6
186 + (-3 * a_2 + 13 * a -13) * c1_2
187 + (2 * a_3 - 25 * a_2 + 72 * a - 61) * c1 / 2
188 + (25 * a_3 - 195 * a_2 + 477 * a - 379) / 12);
189
190 T y_2 = y * y;
191 T y_3 = y_2 * y;
192 T y_4 = y_2 * y_2;
193 result = y + c1 + (c2 / y) + (c3 / y_2) + (c4 / y_3) + (c5 / y_4);
194 BOOST_MATH_INSTRUMENT_VARIABLE(result);
195 if(b < 1e-28f)
196 *p_has_10_digits = true;
197 }
198 }
199 else
200 {
201 // DiDonato and Morris Eq 31:
202 T s = find_inverse_s(p, q);
203
204 BOOST_MATH_INSTRUMENT_VARIABLE(s);
205
206 T s_2 = s * s;
207 T s_3 = s_2 * s;
208 T s_4 = s_2 * s_2;
209 T s_5 = s_4 * s;
210 T ra = sqrt(a);
211
212 BOOST_MATH_INSTRUMENT_VARIABLE(ra);
213
214 T w = a + s * ra + (s * s -1) / 3;
215 w += (s_3 - 7 * s) / (36 * ra);
216 w -= (3 * s_4 + 7 * s_2 - 16) / (810 * a);
217 w += (9 * s_5 + 256 * s_3 - 433 * s) / (38880 * a * ra);
218
219 BOOST_MATH_INSTRUMENT_VARIABLE(w);
220
221 if((a >= 500) && (fabs(1 - w / a) < 1e-6))
222 {
223 result = w;
224 *p_has_10_digits = true;
225 BOOST_MATH_INSTRUMENT_VARIABLE(result);
226 }
227 else if (p > 0.5)
228 {
229 if(w < 3 * a)
230 {
231 result = w;
232 BOOST_MATH_INSTRUMENT_VARIABLE(result);
233 }
234 else
235 {
236 T D = (std::max)(T(2), T(a * (a - 1)));
237 T lg = boost::math::lgamma(a, pol);
238 T lb = log(q) + lg;
239 if(lb < -D * T(2.3))
240 {
241 // DiDonato and Morris Eq 25:
242 T y = -lb;
243 T c1 = (a - 1) * log(y);
244 T c1_2 = c1 * c1;
245 T c1_3 = c1_2 * c1;
246 T c1_4 = c1_2 * c1_2;
247 T a_2 = a * a;
248 T a_3 = a_2 * a;
249
250 T c2 = (a - 1) * (1 + c1);
251 T c3 = (a - 1) * (-(c1_2 / 2) + (a - 2) * c1 + (3 * a - 5) / 2);
252 T c4 = (a - 1) * ((c1_3 / 3) - (3 * a - 5) * c1_2 / 2 + (a_2 - 6 * a + 7) * c1 + (11 * a_2 - 46 * a + 47) / 6);
253 T c5 = (a - 1) * (-(c1_4 / 4)
254 + (11 * a - 17) * c1_3 / 6
255 + (-3 * a_2 + 13 * a -13) * c1_2
256 + (2 * a_3 - 25 * a_2 + 72 * a - 61) * c1 / 2
257 + (25 * a_3 - 195 * a_2 + 477 * a - 379) / 12);
258
259 T y_2 = y * y;
260 T y_3 = y_2 * y;
261 T y_4 = y_2 * y_2;
262 result = y + c1 + (c2 / y) + (c3 / y_2) + (c4 / y_3) + (c5 / y_4);
263 BOOST_MATH_INSTRUMENT_VARIABLE(result);
264 }
265 else
266 {
267 // DiDonato and Morris Eq 33:
268 T u = -lb + (a - 1) * log(w) - log(1 + (1 - a) / (1 + w));
269 result = -lb + (a - 1) * log(u) - log(1 + (1 - a) / (1 + u));
270 BOOST_MATH_INSTRUMENT_VARIABLE(result);
271 }
272 }
273 }
274 else
275 {
276 T z = w;
277 T ap1 = a + 1;
278 T ap2 = a + 2;
279 if(w < 0.15f * ap1)
280 {
281 // DiDonato and Morris Eq 35:
282 T v = log(p) + boost::math::lgamma(ap1, pol);
283 z = exp((v + w) / a);
284 s = boost::math::log1p(z / ap1 * (1 + z / ap2), pol);
285 z = exp((v + z - s) / a);
286 s = boost::math::log1p(z / ap1 * (1 + z / ap2), pol);
287 z = exp((v + z - s) / a);
288 s = boost::math::log1p(z / ap1 * (1 + z / ap2 * (1 + z / (a + 3))), pol);
289 z = exp((v + z - s) / a);
290 BOOST_MATH_INSTRUMENT_VARIABLE(z);
291 }
292
293 if((z <= 0.01 * ap1) || (z > 0.7 * ap1))
294 {
295 result = z;
296 if(z <= T(0.002) * ap1)
297 *p_has_10_digits = true;
298 BOOST_MATH_INSTRUMENT_VARIABLE(result);
299 }
300 else
301 {
302 // DiDonato and Morris Eq 36:
303 T ls = log(didonato_SN(a, z, 100, T(1e-4)));
304 T v = log(p) + boost::math::lgamma(ap1, pol);
305 z = exp((v + z - ls) / a);
306 result = z * (1 - (a * log(z) - z - v + ls) / (a - z));
307
308 BOOST_MATH_INSTRUMENT_VARIABLE(result);
309 }
310 }
311 }
312 return result;
313}
314
315template <class T, class Policy>
316struct gamma_p_inverse_func
317{
318 gamma_p_inverse_func(T a_, T p_, bool inv) : a(a_), p(p_), invert(inv)
319 {
320 //
321 // If p is too near 1 then P(x) - p suffers from cancellation
322 // errors causing our root-finding algorithms to "thrash", better
323 // to invert in this case and calculate Q(x) - (1-p) instead.
324 //
325 // Of course if p is *very* close to 1, then the answer we get will
326 // be inaccurate anyway (because there's not enough information in p)
327 // but at least we will converge on the (inaccurate) answer quickly.
328 //
329 if(p > T(0.9))
330 {
331 p = 1 - p;
332 invert = !invert;
333 }
334 }
335
336 boost::math::tuple<T, T, T> operator()(const T& x)const
337 {
338 BOOST_FPU_EXCEPTION_GUARD
339 //
340 // Calculate P(x) - p and the first two derivates, or if the invert
341 // flag is set, then Q(x) - q and it's derivatives.
342 //
343 typedef typename policies::evaluation<T, Policy>::type value_type;
344 // typedef typename lanczos::lanczos<T, Policy>::type evaluation_type;
345 typedef typename policies::normalise<
346 Policy,
347 policies::promote_float<false>,
348 policies::promote_double<false>,
349 policies::discrete_quantile<>,
350 policies::assert_undefined<> >::type forwarding_policy;
351
352 BOOST_MATH_STD_USING // For ADL of std functions.
353
354 T f, f1;
355 value_type ft;
356 f = static_cast<T>(boost::math::detail::gamma_incomplete_imp(
357 static_cast<value_type>(a),
358 static_cast<value_type>(x),
359 true, invert,
360 forwarding_policy(), &ft));
361 f1 = static_cast<T>(ft);
362 T f2;
363 T div = (a - x - 1) / x;
364 f2 = f1;
365 if(fabs(div) > 1)
366 {
367 // split if statement to address M1 mac clang bug;
368 // see issue 826
369 if (tools::max_value<T>() / fabs(div) < f2)
370 {
371 // overflow:
372 f2 = -tools::max_value<T>() / 2;
373 }
374 else
375 {
376 f2 *= div;
377 }
378 }
379 else
380 {
381 f2 *= div;
382 }
383
384 if(invert)
385 {
386 f1 = -f1;
387 f2 = -f2;
388 }
389
390 return boost::math::make_tuple(static_cast<T>(f - p), f1, f2);
391 }
392private:
393 T a, p;
394 bool invert;
395};
396
397template <class T, class Policy>
398T gamma_p_inv_imp(T a, T p, const Policy& pol)
399{
400 BOOST_MATH_STD_USING // ADL of std functions.
401
402 static const char* function = "boost::math::gamma_p_inv<%1%>(%1%, %1%)";
403
404 BOOST_MATH_INSTRUMENT_VARIABLE(a);
405 BOOST_MATH_INSTRUMENT_VARIABLE(p);
406
407 if(a <= 0)
408 return policies::raise_domain_error<T>(function, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a, pol);
409 if((p < 0) || (p > 1))
410 return policies::raise_domain_error<T>(function, "Probability must be in the range [0,1] in the incomplete gamma function inverse (got p=%1%).", p, pol);
411 if(p == 1)
412 return policies::raise_overflow_error<T>(function, nullptr, Policy());
413 if(p == 0)
414 return 0;
415 bool has_10_digits;
416 T guess = detail::find_inverse_gamma<T>(a, p, 1 - p, pol, &has_10_digits);
417 if((policies::digits<T, Policy>() <= 36) && has_10_digits)
418 return guess;
419 T lower = tools::min_value<T>();
420 if(guess <= lower)
421 guess = tools::min_value<T>();
422 BOOST_MATH_INSTRUMENT_VARIABLE(guess);
423 //
424 // Work out how many digits to converge to, normally this is
425 // 2/3 of the digits in T, but if the first derivative is very
426 // large convergence is slow, so we'll bump it up to full
427 // precision to prevent premature termination of the root-finding routine.
428 //
429 unsigned digits = policies::digits<T, Policy>();
430 if(digits < 30)
431 {
432 digits *= 2;
433 digits /= 3;
434 }
435 else
436 {
437 digits /= 2;
438 digits -= 1;
439 }
440 if((a < T(0.125)) && (fabs(gamma_p_derivative(a, guess, pol)) > 1 / sqrt(tools::epsilon<T>())))
441 digits = policies::digits<T, Policy>() - 2;
442 //
443 // Go ahead and iterate:
444 //
445 std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
446 guess = tools::halley_iterate(
447 detail::gamma_p_inverse_func<T, Policy>(a, p, false),
448 guess,
449 lower,
450 tools::max_value<T>(),
451 digits,
452 max_iter);
453 policies::check_root_iterations<T>(function, max_iter, pol);
454 BOOST_MATH_INSTRUMENT_VARIABLE(guess);
455 if(guess == lower)
456 guess = policies::raise_underflow_error<T>(function, "Expected result known to be non-zero, but is smaller than the smallest available number.", pol);
457 return guess;
458}
459
460template <class T, class Policy>
461T gamma_q_inv_imp(T a, T q, const Policy& pol)
462{
463 BOOST_MATH_STD_USING // ADL of std functions.
464
465 static const char* function = "boost::math::gamma_q_inv<%1%>(%1%, %1%)";
466
467 if(a <= 0)
468 return policies::raise_domain_error<T>(function, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a, pol);
469 if((q < 0) || (q > 1))
470 return policies::raise_domain_error<T>(function, "Probability must be in the range [0,1] in the incomplete gamma function inverse (got q=%1%).", q, pol);
471 if(q == 0)
472 return policies::raise_overflow_error<T>(function, nullptr, Policy());
473 if(q == 1)
474 return 0;
475 bool has_10_digits;
476 T guess = detail::find_inverse_gamma<T>(a, 1 - q, q, pol, &has_10_digits);
477 if((policies::digits<T, Policy>() <= 36) && has_10_digits)
478 return guess;
479 T lower = tools::min_value<T>();
480 if(guess <= lower)
481 guess = tools::min_value<T>();
482 //
483 // Work out how many digits to converge to, normally this is
484 // 2/3 of the digits in T, but if the first derivative is very
485 // large convergence is slow, so we'll bump it up to full
486 // precision to prevent premature termination of the root-finding routine.
487 //
488 unsigned digits = policies::digits<T, Policy>();
489 if(digits < 30)
490 {
491 digits *= 2;
492 digits /= 3;
493 }
494 else
495 {
496 digits /= 2;
497 digits -= 1;
498 }
499 if((a < 0.125) && (fabs(gamma_p_derivative(a, guess, pol)) > 1 / sqrt(tools::epsilon<T>())))
500 digits = policies::digits<T, Policy>();
501 //
502 // Go ahead and iterate:
503 //
504 std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
505 guess = tools::halley_iterate(
506 detail::gamma_p_inverse_func<T, Policy>(a, q, true),
507 guess,
508 lower,
509 tools::max_value<T>(),
510 digits,
511 max_iter);
512 policies::check_root_iterations<T>(function, max_iter, pol);
513 if(guess == lower)
514 guess = policies::raise_underflow_error<T>(function, "Expected result known to be non-zero, but is smaller than the smallest available number.", pol);
515 return guess;
516}
517
518} // namespace detail
519
520template <class T1, class T2, class Policy>
521inline typename tools::promote_args<T1, T2>::type
522 gamma_p_inv(T1 a, T2 p, const Policy& pol)
523{
524 typedef typename tools::promote_args<T1, T2>::type result_type;
525 return detail::gamma_p_inv_imp(
526 static_cast<result_type>(a),
527 static_cast<result_type>(p), pol);
528}
529
530template <class T1, class T2, class Policy>
531inline typename tools::promote_args<T1, T2>::type
532 gamma_q_inv(T1 a, T2 p, const Policy& pol)
533{
534 typedef typename tools::promote_args<T1, T2>::type result_type;
535 return detail::gamma_q_inv_imp(
536 static_cast<result_type>(a),
537 static_cast<result_type>(p), pol);
538}
539
540template <class T1, class T2>
541inline typename tools::promote_args<T1, T2>::type
542 gamma_p_inv(T1 a, T2 p)
543{
544 return gamma_p_inv(a, p, policies::policy<>());
545}
546
547template <class T1, class T2>
548inline typename tools::promote_args<T1, T2>::type
549 gamma_q_inv(T1 a, T2 p)
550{
551 return gamma_q_inv(a, p, policies::policy<>());
552}
553
554} // namespace math
555} // namespace boost
556
557#endif // BOOST_MATH_SPECIAL_FUNCTIONS_IGAMMA_INVERSE_HPP
558
559
560
561

source code of include/boost/math/special_functions/detail/igamma_inverse.hpp