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_BETA_HPP
7#define BOOST_MATH_SPECIAL_BETA_HPP
8
9#ifdef _MSC_VER
10#pragma once
11#endif
12
13#include <boost/math/special_functions/math_fwd.hpp>
14#include <boost/math/tools/config.hpp>
15#include <boost/math/special_functions/gamma.hpp>
16#include <boost/math/special_functions/binomial.hpp>
17#include <boost/math/special_functions/factorials.hpp>
18#include <boost/math/special_functions/erf.hpp>
19#include <boost/math/special_functions/log1p.hpp>
20#include <boost/math/special_functions/expm1.hpp>
21#include <boost/math/special_functions/trunc.hpp>
22#include <boost/math/tools/roots.hpp>
23#include <boost/math/tools/assert.hpp>
24#include <cmath>
25
26namespace boost{ namespace math{
27
28namespace detail{
29
30//
31// Implementation of Beta(a,b) using the Lanczos approximation:
32//
33template <class T, class Lanczos, class Policy>
34T beta_imp(T a, T b, const Lanczos&, const Policy& pol)
35{
36 BOOST_MATH_STD_USING // for ADL of std names
37
38 if(a <= 0)
39 return policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got a=%1%).", a, pol);
40 if(b <= 0)
41 return policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got b=%1%).", b, pol);
42
43 T result; // LCOV_EXCL_LINE
44
45 T prefix = 1;
46 T c = a + b;
47
48 // Special cases:
49 if((c == a) && (b < tools::epsilon<T>()))
50 return 1 / b;
51 else if((c == b) && (a < tools::epsilon<T>()))
52 return 1 / a;
53 if(b == 1)
54 return 1/a;
55 else if(a == 1)
56 return 1/b;
57 else if(c < tools::epsilon<T>())
58 {
59 result = c / a;
60 result /= b;
61 return result;
62 }
63
64 /*
65 //
66 // This code appears to be no longer necessary: it was
67 // used to offset errors introduced from the Lanczos
68 // approximation, but the current Lanczos approximations
69 // are sufficiently accurate for all z that we can ditch
70 // this. It remains in the file for future reference...
71 //
72 // If a or b are less than 1, shift to greater than 1:
73 if(a < 1)
74 {
75 prefix *= c / a;
76 c += 1;
77 a += 1;
78 }
79 if(b < 1)
80 {
81 prefix *= c / b;
82 c += 1;
83 b += 1;
84 }
85 */
86
87 if(a < b)
88 std::swap(a, b);
89
90 // Lanczos calculation:
91 T agh = static_cast<T>(a + Lanczos::g() - 0.5f);
92 T bgh = static_cast<T>(b + Lanczos::g() - 0.5f);
93 T cgh = static_cast<T>(c + Lanczos::g() - 0.5f);
94 result = Lanczos::lanczos_sum_expG_scaled(a) * (Lanczos::lanczos_sum_expG_scaled(b) / Lanczos::lanczos_sum_expG_scaled(c));
95 T ambh = a - 0.5f - b;
96 if((fabs(b * ambh) < (cgh * 100)) && (a > 100))
97 {
98 // Special case where the base of the power term is close to 1
99 // compute (1+x)^y instead:
100 result *= exp(ambh * boost::math::log1p(-b / cgh, pol));
101 }
102 else
103 {
104 result *= pow(agh / cgh, a - T(0.5) - b);
105 }
106 if(cgh > 1e10f)
107 // this avoids possible overflow, but appears to be marginally less accurate:
108 result *= pow((agh / cgh) * (bgh / cgh), b);
109 else
110 result *= pow((agh * bgh) / (cgh * cgh), b);
111 result *= sqrt(boost::math::constants::e<T>() / bgh);
112
113 // If a and b were originally less than 1 we need to scale the result:
114 result *= prefix;
115
116 return result;
117} // template <class T, class Lanczos> beta_imp(T a, T b, const Lanczos&)
118
119//
120// Generic implementation of Beta(a,b) without Lanczos approximation support
121// (Caution this is slow!!!):
122//
123template <class T, class Policy>
124T beta_imp(T a, T b, const lanczos::undefined_lanczos& l, const Policy& pol)
125{
126 BOOST_MATH_STD_USING
127
128 if(a <= 0)
129 return policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got a=%1%).", a, pol);
130 if(b <= 0)
131 return policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got b=%1%).", b, pol);
132
133 const T c = a + b;
134
135 // Special cases:
136 if ((c == a) && (b < tools::epsilon<T>()))
137 return 1 / b;
138 else if ((c == b) && (a < tools::epsilon<T>()))
139 return 1 / a;
140 if (b == 1)
141 return 1 / a;
142 else if (a == 1)
143 return 1 / b;
144 else if (c < tools::epsilon<T>())
145 {
146 T result = c / a;
147 result /= b;
148 return result;
149 }
150
151 // Regular cases start here:
152 const T min_sterling = minimum_argument_for_bernoulli_recursion<T>();
153
154 long shift_a = 0;
155 long shift_b = 0;
156
157 if(a < min_sterling)
158 shift_a = 1 + ltrunc(min_sterling - a);
159 if(b < min_sterling)
160 shift_b = 1 + ltrunc(min_sterling - b);
161 long shift_c = shift_a + shift_b;
162
163 if ((shift_a == 0) && (shift_b == 0))
164 {
165 return pow(a / c, a) * pow(b / c, b) * scaled_tgamma_no_lanczos(a, pol) * scaled_tgamma_no_lanczos(b, pol) / scaled_tgamma_no_lanczos(c, pol);
166 }
167 else if ((a < 1) && (b < 1))
168 {
169 return boost::math::tgamma(a, pol) * (boost::math::tgamma(b, pol) / boost::math::tgamma(c));
170 }
171 else if(a < 1)
172 return boost::math::tgamma(a, pol) * boost::math::tgamma_delta_ratio(b, a, pol);
173 else if(b < 1)
174 return boost::math::tgamma(b, pol) * boost::math::tgamma_delta_ratio(a, b, pol);
175 else
176 {
177 T result = beta_imp(T(a + shift_a), T(b + shift_b), l, pol);
178 //
179 // Recursion:
180 //
181 for (long i = 0; i < shift_c; ++i)
182 {
183 result *= c + i;
184 if (i < shift_a)
185 result /= a + i;
186 if (i < shift_b)
187 result /= b + i;
188 }
189 return result;
190 }
191
192} // template <class T>T beta_imp(T a, T b, const lanczos::undefined_lanczos& l)
193
194
195//
196// Compute the leading power terms in the incomplete Beta:
197//
198// (x^a)(y^b)/Beta(a,b) when normalised, and
199// (x^a)(y^b) otherwise.
200//
201// Almost all of the error in the incomplete beta comes from this
202// function: particularly when a and b are large. Computing large
203// powers are *hard* though, and using logarithms just leads to
204// horrendous cancellation errors.
205//
206template <class T, class Lanczos, class Policy>
207T ibeta_power_terms(T a,
208 T b,
209 T x,
210 T y,
211 const Lanczos&,
212 bool normalised,
213 const Policy& pol,
214 T prefix = 1,
215 const char* function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)")
216{
217 BOOST_MATH_STD_USING
218
219 if(!normalised)
220 {
221 // can we do better here?
222 return pow(x, a) * pow(y, b);
223 }
224
225 T result; // LCOV_EXCL_LINE
226
227 T c = a + b;
228
229 // combine power terms with Lanczos approximation:
230 T agh = static_cast<T>(a + Lanczos::g() - 0.5f);
231 T bgh = static_cast<T>(b + Lanczos::g() - 0.5f);
232 T cgh = static_cast<T>(c + Lanczos::g() - 0.5f);
233 if ((a < tools::min_value<T>()) || (b < tools::min_value<T>()))
234 result = 0; // denominator overflows in this case
235 else
236 result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));
237 result *= prefix;
238 // combine with the leftover terms from the Lanczos approximation:
239 result *= sqrt(bgh / boost::math::constants::e<T>());
240 result *= sqrt(agh / cgh);
241
242 // l1 and l2 are the base of the exponents minus one:
243 T l1 = (x * b - y * agh) / agh;
244 T l2 = (y * a - x * bgh) / bgh;
245 if(((std::min)(fabs(l1), fabs(l2)) < 0.2))
246 {
247 // when the base of the exponent is very near 1 we get really
248 // gross errors unless extra care is taken:
249 if((l1 * l2 > 0) || ((std::min)(a, b) < 1))
250 {
251 //
252 // This first branch handles the simple cases where either:
253 //
254 // * The two power terms both go in the same direction
255 // (towards zero or towards infinity). In this case if either
256 // term overflows or underflows, then the product of the two must
257 // do so also.
258 // *Alternatively if one exponent is less than one, then we
259 // can't productively use it to eliminate overflow or underflow
260 // from the other term. Problems with spurious overflow/underflow
261 // can't be ruled out in this case, but it is *very* unlikely
262 // since one of the power terms will evaluate to a number close to 1.
263 //
264 if(fabs(l1) < 0.1)
265 {
266 result *= exp(a * boost::math::log1p(l1, pol));
267 BOOST_MATH_INSTRUMENT_VARIABLE(result);
268 }
269 else
270 {
271 result *= pow((x * cgh) / agh, a);
272 BOOST_MATH_INSTRUMENT_VARIABLE(result);
273 }
274 if(fabs(l2) < 0.1)
275 {
276 result *= exp(b * boost::math::log1p(l2, pol));
277 BOOST_MATH_INSTRUMENT_VARIABLE(result);
278 }
279 else
280 {
281 result *= pow((y * cgh) / bgh, b);
282 BOOST_MATH_INSTRUMENT_VARIABLE(result);
283 }
284 }
285 else if((std::max)(fabs(l1), fabs(l2)) < 0.5)
286 {
287 //
288 // Both exponents are near one and both the exponents are
289 // greater than one and further these two
290 // power terms tend in opposite directions (one towards zero,
291 // the other towards infinity), so we have to combine the terms
292 // to avoid any risk of overflow or underflow.
293 //
294 // We do this by moving one power term inside the other, we have:
295 //
296 // (1 + l1)^a * (1 + l2)^b
297 // = ((1 + l1)*(1 + l2)^(b/a))^a
298 // = (1 + l1 + l3 + l1*l3)^a ; l3 = (1 + l2)^(b/a) - 1
299 // = exp((b/a) * log(1 + l2)) - 1
300 //
301 // The tricky bit is deciding which term to move inside :-)
302 // By preference we move the larger term inside, so that the
303 // size of the largest exponent is reduced. However, that can
304 // only be done as long as l3 (see above) is also small.
305 //
306 bool small_a = a < b;
307 T ratio = b / a;
308 if((small_a && (ratio * l2 < 0.1)) || (!small_a && (l1 / ratio > 0.1)))
309 {
310 T l3 = boost::math::expm1(ratio * boost::math::log1p(l2, pol), pol);
311 l3 = l1 + l3 + l3 * l1;
312 l3 = a * boost::math::log1p(l3, pol);
313 result *= exp(l3);
314 BOOST_MATH_INSTRUMENT_VARIABLE(result);
315 }
316 else
317 {
318 T l3 = boost::math::expm1(boost::math::log1p(l1, pol) / ratio, pol);
319 l3 = l2 + l3 + l3 * l2;
320 l3 = b * boost::math::log1p(l3, pol);
321 result *= exp(l3);
322 BOOST_MATH_INSTRUMENT_VARIABLE(result);
323 }
324 }
325 else if(fabs(l1) < fabs(l2))
326 {
327 // First base near 1 only:
328 T l = a * boost::math::log1p(l1, pol)
329 + b * log((y * cgh) / bgh);
330 if((l <= tools::log_min_value<T>()) || (l >= tools::log_max_value<T>()))
331 {
332 l += log(result);
333 if(l >= tools::log_max_value<T>())
334 return policies::raise_overflow_error<T>(function, nullptr, pol); // LCOV_EXCL_LINE we can probably never get here, probably!
335 result = exp(l);
336 }
337 else
338 result *= exp(l);
339 BOOST_MATH_INSTRUMENT_VARIABLE(result);
340 }
341 else
342 {
343 // Second base near 1 only:
344 T l = b * boost::math::log1p(l2, pol)
345 + a * log((x * cgh) / agh);
346 if((l <= tools::log_min_value<T>()) || (l >= tools::log_max_value<T>()))
347 {
348 l += log(result);
349 if(l >= tools::log_max_value<T>())
350 return policies::raise_overflow_error<T>(function, nullptr, pol); // LCOV_EXCL_LINE we can probably never get here, probably!
351 result = exp(l);
352 }
353 else
354 result *= exp(l);
355 BOOST_MATH_INSTRUMENT_VARIABLE(result);
356 }
357 }
358 else
359 {
360 // general case:
361 T b1 = (x * cgh) / agh;
362 T b2 = (y * cgh) / bgh;
363 l1 = a * log(b1);
364 l2 = b * log(b2);
365 BOOST_MATH_INSTRUMENT_VARIABLE(b1);
366 BOOST_MATH_INSTRUMENT_VARIABLE(b2);
367 BOOST_MATH_INSTRUMENT_VARIABLE(l1);
368 BOOST_MATH_INSTRUMENT_VARIABLE(l2);
369 if((l1 >= tools::log_max_value<T>())
370 || (l1 <= tools::log_min_value<T>())
371 || (l2 >= tools::log_max_value<T>())
372 || (l2 <= tools::log_min_value<T>())
373 )
374 {
375 // Oops, under/overflow, sidestep if we can:
376 if(a < b)
377 {
378 T p1 = pow(b2, b / a);
379 T l3 = (b1 != 0) && (p1 != 0) ? (a * (log(b1) + log(p1))) : tools::max_value<T>(); // arbitrary large value if the logs would fail!
380 if((l3 < tools::log_max_value<T>())
381 && (l3 > tools::log_min_value<T>()))
382 {
383 result *= pow(p1 * b1, a);
384 }
385 else
386 {
387 l2 += l1 + log(result);
388 if(l2 >= tools::log_max_value<T>())
389 return policies::raise_overflow_error<T>(function, nullptr, pol); // LCOV_EXCL_LINE we can probably never get here, probably!
390 result = exp(l2);
391 }
392 }
393 else
394 {
395 // This protects against spurious overflow in a/b:
396 T p1 = (b1 < 1) && (b < 1) && (tools::max_value<T>() * b < a) ? static_cast<T>(0) : static_cast<T>(pow(b1, a / b));
397 T l3 = (p1 != 0) && (b2 != 0) ? (log(p1) + log(b2)) * b : tools::max_value<T>(); // arbitrary large value if the logs would fail!
398 if((l3 < tools::log_max_value<T>())
399 && (l3 > tools::log_min_value<T>()))
400 {
401 result *= pow(p1 * b2, b);
402 }
403 else if(result != 0) // we can elude the calculation below if we're already going to be zero
404 {
405 l2 += l1 + log(result);
406 if(l2 >= tools::log_max_value<T>())
407 return policies::raise_overflow_error<T>(function, nullptr, pol); // LCOV_EXCL_LINE we can probably never get here, probably!
408 result = exp(l2);
409 }
410 }
411 BOOST_MATH_INSTRUMENT_VARIABLE(result);
412 }
413 else
414 {
415 // finally the normal case:
416 result *= pow(b1, a) * pow(b2, b);
417 BOOST_MATH_INSTRUMENT_VARIABLE(result);
418 }
419 }
420
421 BOOST_MATH_INSTRUMENT_VARIABLE(result);
422
423 if (0 == result)
424 {
425 if ((a > 1) && (x == 0))
426 return result; // true zero LCOV_EXCL_LINE we can probably never get here
427 if ((b > 1) && (y == 0))
428 return result; // true zero LCOV_EXCL_LINE we can probably never get here
429 return boost::math::policies::raise_underflow_error<T>(function, nullptr, pol);
430 }
431
432 return result;
433}
434//
435// Compute the leading power terms in the incomplete Beta:
436//
437// (x^a)(y^b)/Beta(a,b) when normalised, and
438// (x^a)(y^b) otherwise.
439//
440// Almost all of the error in the incomplete beta comes from this
441// function: particularly when a and b are large. Computing large
442// powers are *hard* though, and using logarithms just leads to
443// horrendous cancellation errors.
444//
445// This version is generic, slow, and does not use the Lanczos approximation.
446//
447template <class T, class Policy>
448T ibeta_power_terms(T a,
449 T b,
450 T x,
451 T y,
452 const boost::math::lanczos::undefined_lanczos& l,
453 bool normalised,
454 const Policy& pol,
455 T prefix = 1,
456 const char* = "boost::math::ibeta<%1%>(%1%, %1%, %1%)")
457{
458 BOOST_MATH_STD_USING
459
460 if(!normalised)
461 {
462 return prefix * pow(x, a) * pow(y, b);
463 }
464
465 T c = a + b;
466
467 const T min_sterling = minimum_argument_for_bernoulli_recursion<T>();
468
469 long shift_a = 0;
470 long shift_b = 0;
471
472 if (a < min_sterling)
473 shift_a = 1 + ltrunc(min_sterling - a);
474 if (b < min_sterling)
475 shift_b = 1 + ltrunc(min_sterling - b);
476
477 if ((shift_a == 0) && (shift_b == 0))
478 {
479 T power1, power2;
480 bool need_logs = false;
481 if (a < b)
482 {
483 BOOST_MATH_IF_CONSTEXPR(std::numeric_limits<T>::has_infinity)
484 {
485 power1 = pow((x * y * c * c) / (a * b), a);
486 power2 = pow((y * c) / b, b - a);
487 }
488 else
489 {
490 // We calculate these logs purely so we can check for overflow in the power functions
491 T l1 = log((x * y * c * c) / (a * b));
492 T l2 = log((y * c) / b);
493 if ((l1 * a > tools::log_min_value<T>()) && (l1 * a < tools::log_max_value<T>()) && (l2 * (b - a) < tools::log_max_value<T>()) && (l2 * (b - a) > tools::log_min_value<T>()))
494 {
495 power1 = pow((x * y * c * c) / (a * b), a);
496 power2 = pow((y * c) / b, b - a);
497 }
498 else
499 {
500 need_logs = true;
501 }
502 }
503 }
504 else
505 {
506 BOOST_MATH_IF_CONSTEXPR(std::numeric_limits<T>::has_infinity)
507 {
508 power1 = pow((x * y * c * c) / (a * b), b);
509 power2 = pow((x * c) / a, a - b);
510 }
511 else
512 {
513 // We calculate these logs purely so we can check for overflow in the power functions
514 T l1 = log((x * y * c * c) / (a * b)) * b;
515 T l2 = log((x * c) / a) * (a - b);
516 if ((l1 * a > tools::log_min_value<T>()) && (l1 * a < tools::log_max_value<T>()) && (l2 * (b - a) < tools::log_max_value<T>()) && (l2 * (b - a) > tools::log_min_value<T>()))
517 {
518 power1 = pow((x * y * c * c) / (a * b), b);
519 power2 = pow((x * c) / a, a - b);
520 }
521 else
522 need_logs = true;
523 }
524 }
525 BOOST_MATH_IF_CONSTEXPR(std::numeric_limits<T>::has_infinity)
526 {
527 if (!(boost::math::isnormal)(power1) || !(boost::math::isnormal)(power2))
528 {
529 need_logs = true;
530 }
531 }
532 if (need_logs)
533 {
534 //
535 // We want:
536 //
537 // (xc / a)^a (yc / b)^b
538 //
539 // But we know that one or other term will over / underflow and combining the logs will be next to useless as that will cause significant cancellation.
540 // If we assume b > a and express z ^ b as(z ^ b / a) ^ a with z = (yc / b) then we can move one power term inside the other :
541 //
542 // ((xc / a) * (yc / b)^(b / a))^a
543 //
544 // However, we're not quite there yet, as the term being exponentiated is quite likely to be close to unity, so let:
545 //
546 // xc / a = 1 + (xb - ya) / a
547 //
548 // analogously let :
549 //
550 // 1 + p = (yc / b) ^ (b / a) = 1 + expm1((b / a) * log1p((ya - xb) / b))
551 //
552 // so putting the two together we have :
553 //
554 // exp(a * log1p((xb - ya) / a + p + p(xb - ya) / a))
555 //
556 // Analogously, when a > b we can just swap all the terms around.
557 //
558 // Finally, there are a few cases (x or y is unity) when the above logic can't be used
559 // or where there is no logarithmic cancellation and accuracy is better just using
560 // the regular formula:
561 //
562 T xc_a = x * c / a;
563 T yc_b = y * c / b;
564 if ((x == 1) || (y == 1) || (fabs(xc_a - 1) > 0.25) || (fabs(yc_b - 1) > 0.25))
565 {
566 // The above logic fails, the result is almost certainly zero:
567 power1 = exp(log(xc_a) * a + log(yc_b) * b);
568 power2 = 1;
569 }
570 else if (b > a)
571 {
572 T p = boost::math::expm1((b / a) * boost::math::log1p((y * a - x * b) / b));
573 power1 = exp(a * boost::math::log1p((x * b - y * a) / a + p * (x * c / a)));
574 power2 = 1;
575 }
576 else
577 {
578 T p = boost::math::expm1((a / b) * boost::math::log1p((x * b - y * a) / a));
579 power1 = exp(b * boost::math::log1p((y * a - x * b) / b + p * (y * c / b)));
580 power2 = 1;
581 }
582 }
583 return prefix * power1 * power2 * scaled_tgamma_no_lanczos(c, pol) / (scaled_tgamma_no_lanczos(a, pol) * scaled_tgamma_no_lanczos(b, pol));
584 }
585
586 T power1 = pow(x, a);
587 T power2 = pow(y, b);
588 T bet = beta_imp(a, b, l, pol);
589
590 if(!(boost::math::isnormal)(power1) || !(boost::math::isnormal)(power2) || !(boost::math::isnormal)(bet))
591 {
592 int shift_c = shift_a + shift_b;
593 T result = ibeta_power_terms(T(a + shift_a), T(b + shift_b), x, y, l, normalised, pol, prefix);
594 if ((boost::math::isnormal)(result))
595 {
596 for (int i = 0; i < shift_c; ++i)
597 {
598 result /= c + i;
599 if (i < shift_a)
600 {
601 result *= a + i;
602 result /= x;
603 }
604 if (i < shift_b)
605 {
606 result *= b + i;
607 result /= y;
608 }
609 }
610 return prefix * result;
611 }
612 else
613 {
614 T log_result = log(x) * a + log(y) * b + log(prefix);
615 if ((boost::math::isnormal)(bet))
616 log_result -= log(bet);
617 else
618 log_result += boost::math::lgamma(c, pol) - boost::math::lgamma(a, pol) - boost::math::lgamma(b, pol);
619 return exp(log_result);
620 }
621 }
622 return prefix * power1 * (power2 / bet);
623}
624//
625// Series approximation to the incomplete beta:
626//
627template <class T>
628struct ibeta_series_t
629{
630 typedef T result_type;
631 ibeta_series_t(T a_, T b_, T x_, T mult) : result(mult), x(x_), apn(a_), poch(1-b_), n(1) {}
632 T operator()()
633 {
634 T r = result / apn;
635 apn += 1;
636 result *= poch * x / n;
637 ++n;
638 poch += 1;
639 return r;
640 }
641private:
642 T result, x, apn, poch;
643 int n;
644};
645
646template <class T, class Lanczos, class Policy>
647T ibeta_series(T a, T b, T x, T s0, const Lanczos&, bool normalised, T* p_derivative, T y, const Policy& pol)
648{
649 BOOST_MATH_STD_USING
650
651 T result;
652
653 BOOST_MATH_ASSERT((p_derivative == 0) || normalised);
654
655 if(normalised)
656 {
657 T c = a + b;
658
659 // incomplete beta power term, combined with the Lanczos approximation:
660 T agh = static_cast<T>(a + Lanczos::g() - 0.5f);
661 T bgh = static_cast<T>(b + Lanczos::g() - 0.5f);
662 T cgh = static_cast<T>(c + Lanczos::g() - 0.5f);
663 if ((a < tools::min_value<T>()) || (b < tools::min_value<T>()))
664 result = 0; // denorms cause overflow in the Lanzos series, result will be zero anyway
665 else
666 result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));
667
668 if (!(boost::math::isfinite)(result))
669 result = 0; // LCOV_EXCL_LINE we can probably never get here, covered already above?
670
671 T l1 = log(cgh / bgh) * (b - 0.5f);
672 T l2 = log(x * cgh / agh) * a;
673 //
674 // Check for over/underflow in the power terms:
675 //
676 if((l1 > tools::log_min_value<T>())
677 && (l1 < tools::log_max_value<T>())
678 && (l2 > tools::log_min_value<T>())
679 && (l2 < tools::log_max_value<T>()))
680 {
681 if(a * b < bgh * 10)
682 result *= exp((b - 0.5f) * boost::math::log1p(a / bgh, pol));
683 else
684 result *= pow(cgh / bgh, T(b - T(0.5)));
685 result *= pow(x * cgh / agh, a);
686 result *= sqrt(agh / boost::math::constants::e<T>());
687
688 if(p_derivative)
689 {
690 *p_derivative = result * pow(y, b);
691 BOOST_MATH_ASSERT(*p_derivative >= 0);
692 }
693 }
694 else
695 {
696 //
697 // Oh dear, we need logs, and this *will* cancel:
698 //
699 if (result != 0) // elude calculation when result will be zero.
700 {
701 result = log(result) + l1 + l2 + (log(agh) - 1) / 2;
702 if (p_derivative)
703 *p_derivative = exp(result + b * log(y));
704 result = exp(result);
705 }
706 }
707 }
708 else
709 {
710 // Non-normalised, just compute the power:
711 result = pow(x, a);
712 }
713 if(result < tools::min_value<T>())
714 return s0; // Safeguard: series can't cope with denorms.
715 ibeta_series_t<T> s(a, b, x, result);
716 std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
717 result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0);
718 policies::check_series_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (with lanczos)", max_iter, pol);
719 return result;
720}
721//
722// Incomplete Beta series again, this time without Lanczos support:
723//
724template <class T, class Policy>
725T ibeta_series(T a, T b, T x, T s0, const boost::math::lanczos::undefined_lanczos& l, bool normalised, T* p_derivative, T y, const Policy& pol)
726{
727 BOOST_MATH_STD_USING
728
729 T result;
730 BOOST_MATH_ASSERT((p_derivative == 0) || normalised);
731
732 if(normalised)
733 {
734 const T min_sterling = minimum_argument_for_bernoulli_recursion<T>();
735
736 long shift_a = 0;
737 long shift_b = 0;
738
739 if (a < min_sterling)
740 shift_a = 1 + ltrunc(min_sterling - a);
741 if (b < min_sterling)
742 shift_b = 1 + ltrunc(min_sterling - b);
743
744 T c = a + b;
745
746 if ((shift_a == 0) && (shift_b == 0))
747 {
748 result = pow(x * c / a, a) * pow(c / b, b) * scaled_tgamma_no_lanczos(c, pol) / (scaled_tgamma_no_lanczos(a, pol) * scaled_tgamma_no_lanczos(b, pol));
749 }
750 else if ((a < 1) && (b > 1))
751 result = pow(x, a) / (boost::math::tgamma(a, pol) * boost::math::tgamma_delta_ratio(b, a, pol));
752 else
753 {
754 T power = pow(x, a);
755 T bet = beta_imp(a, b, l, pol);
756 if (!(boost::math::isnormal)(power) || !(boost::math::isnormal)(bet))
757 {
758 result = exp(a * log(x) + boost::math::lgamma(c, pol) - boost::math::lgamma(a, pol) - boost::math::lgamma(b, pol));
759 }
760 else
761 result = power / bet;
762 }
763 if(p_derivative)
764 {
765 *p_derivative = result * pow(y, b);
766 BOOST_MATH_ASSERT(*p_derivative >= 0);
767 }
768 }
769 else
770 {
771 // Non-normalised, just compute the power:
772 result = pow(x, a);
773 }
774 if(result < tools::min_value<T>())
775 return s0; // Safeguard: series can't cope with denorms.
776 ibeta_series_t<T> s(a, b, x, result);
777 std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
778 result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0);
779 policies::check_series_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (without lanczos)", max_iter, pol);
780 return result;
781}
782
783//
784// Continued fraction for the incomplete beta:
785//
786template <class T>
787struct ibeta_fraction2_t
788{
789 typedef std::pair<T, T> result_type;
790
791 ibeta_fraction2_t(T a_, T b_, T x_, T y_) : a(a_), b(b_), x(x_), y(y_), m(0) {}
792
793 result_type operator()()
794 {
795 T aN = (a + m - 1) * (a + b + m - 1) * m * (b - m) * x * x;
796 T denom = (a + 2 * m - 1);
797 aN /= denom * denom;
798
799 T bN = static_cast<T>(m);
800 bN += (m * (b - m) * x) / (a + 2*m - 1);
801 bN += ((a + m) * (a * y - b * x + 1 + m *(2 - x))) / (a + 2*m + 1);
802
803 ++m;
804
805 return std::make_pair(aN, bN);
806 }
807
808private:
809 T a, b, x, y;
810 int m;
811};
812//
813// Evaluate the incomplete beta via the continued fraction representation:
814//
815template <class T, class Policy>
816inline T ibeta_fraction2(T a, T b, T x, T y, const Policy& pol, bool normalised, T* p_derivative)
817{
818 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
819 BOOST_MATH_STD_USING
820 T result = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol);
821 if(p_derivative)
822 {
823 *p_derivative = result;
824 BOOST_MATH_ASSERT(*p_derivative >= 0);
825 }
826 if(result == 0)
827 return result;
828
829 ibeta_fraction2_t<T> f(a, b, x, y);
830 T fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon<T, Policy>());
831 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
832 BOOST_MATH_INSTRUMENT_VARIABLE(result);
833 return result / fract;
834}
835//
836// Computes the difference between ibeta(a,b,x) and ibeta(a+k,b,x):
837//
838template <class T, class Policy>
839T ibeta_a_step(T a, T b, T x, T y, int k, const Policy& pol, bool normalised, T* p_derivative)
840{
841 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
842
843 BOOST_MATH_INSTRUMENT_VARIABLE(k);
844
845 T prefix = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol);
846 if(p_derivative)
847 {
848 *p_derivative = prefix;
849 BOOST_MATH_ASSERT(*p_derivative >= 0);
850 }
851 prefix /= a;
852 if(prefix == 0)
853 return prefix;
854 T sum = 1;
855 T term = 1;
856 // series summation from 0 to k-1:
857 for(int i = 0; i < k-1; ++i)
858 {
859 term *= (a+b+i) * x / (a+i+1);
860 sum += term;
861 }
862 prefix *= sum;
863
864 return prefix;
865}
866//
867// This function is only needed for the non-regular incomplete beta,
868// it computes the delta in:
869// beta(a,b,x) = prefix + delta * beta(a+k,b,x)
870// it is currently only called for small k.
871//
872template <class T>
873inline T rising_factorial_ratio(T a, T b, int k)
874{
875 // calculate:
876 // (a)(a+1)(a+2)...(a+k-1)
877 // _______________________
878 // (b)(b+1)(b+2)...(b+k-1)
879
880 // This is only called with small k, for large k
881 // it is grossly inefficient, do not use outside it's
882 // intended purpose!!!
883 BOOST_MATH_INSTRUMENT_VARIABLE(k);
884 BOOST_MATH_ASSERT(k > 0);
885
886 T result = 1;
887 for(int i = 0; i < k; ++i)
888 result *= (a+i) / (b+i);
889 return result;
890}
891//
892// Routine for a > 15, b < 1
893//
894// Begin by figuring out how large our table of Pn's should be,
895// quoted accuracies are "guesstimates" based on empirical observation.
896// Note that the table size should never exceed the size of our
897// tables of factorials.
898//
899template <class T>
900struct Pn_size
901{
902 // This is likely to be enough for ~35-50 digit accuracy
903 // but it's hard to quantify exactly:
904 static constexpr unsigned value =
905 ::boost::math::max_factorial<T>::value >= 100 ? 50
906 : ::boost::math::max_factorial<T>::value >= ::boost::math::max_factorial<double>::value ? 30
907 : ::boost::math::max_factorial<T>::value >= ::boost::math::max_factorial<float>::value ? 15 : 1;
908 static_assert(::boost::math::max_factorial<T>::value >= ::boost::math::max_factorial<float>::value, "Type does not provide for 35-50 digits of accuracy.");
909};
910template <>
911struct Pn_size<float>
912{
913 static constexpr unsigned value = 15; // ~8-15 digit accuracy
914 static_assert(::boost::math::max_factorial<float>::value >= 30, "Type does not provide for 8-15 digits of accuracy.");
915};
916template <>
917struct Pn_size<double>
918{
919 static constexpr unsigned value = 30; // 16-20 digit accuracy
920 static_assert(::boost::math::max_factorial<double>::value >= 60, "Type does not provide for 16-20 digits of accuracy.");
921};
922template <>
923struct Pn_size<long double>
924{
925 static constexpr unsigned value = 50; // ~35-50 digit accuracy
926 static_assert(::boost::math::max_factorial<long double>::value >= 100, "Type does not provide for ~35-50 digits of accuracy");
927};
928
929template <class T, class Policy>
930T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const Policy& pol, bool normalised)
931{
932 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
933 BOOST_MATH_STD_USING
934 //
935 // This is DiDonato and Morris's BGRAT routine, see Eq's 9 through 9.6.
936 //
937 // Some values we'll need later, these are Eq 9.1:
938 //
939 T bm1 = b - 1;
940 T t = a + bm1 / 2;
941 T lx, u; // LCOV_EXCL_LINE
942 if(y < 0.35)
943 lx = boost::math::log1p(-y, pol);
944 else
945 lx = log(x);
946 u = -t * lx;
947 // and from from 9.2:
948 T prefix; // LCOV_EXCL_LINE
949 T h = regularised_gamma_prefix(b, u, pol, lanczos_type());
950 if(h <= tools::min_value<T>())
951 return s0;
952 if(normalised)
953 {
954 prefix = h / boost::math::tgamma_delta_ratio(a, b, pol);
955 prefix /= pow(t, b);
956 }
957 else
958 {
959 prefix = full_igamma_prefix(b, u, pol) / pow(t, b);
960 }
961 prefix *= mult;
962 //
963 // now we need the quantity Pn, unfortunately this is computed
964 // recursively, and requires a full history of all the previous values
965 // so no choice but to declare a big table and hope it's big enough...
966 //
967 T p[ ::boost::math::detail::Pn_size<T>::value ] = { 1 }; // see 9.3.
968 //
969 // Now an initial value for J, see 9.6:
970 //
971 T j = boost::math::gamma_q(b, u, pol) / h;
972 //
973 // Now we can start to pull things together and evaluate the sum in Eq 9:
974 //
975 T sum = s0 + prefix * j; // Value at N = 0
976 // some variables we'll need:
977 unsigned tnp1 = 1; // 2*N+1
978 T lx2 = lx / 2;
979 lx2 *= lx2;
980 T lxp = 1;
981 T t4 = 4 * t * t;
982 T b2n = b;
983
984 for(unsigned n = 1; n < sizeof(p)/sizeof(p[0]); ++n)
985 {
986 /*
987 // debugging code, enable this if you want to determine whether
988 // the table of Pn's is large enough...
989 //
990 static int max_count = 2;
991 if(n > max_count)
992 {
993 max_count = n;
994 std::cerr << "Max iterations in BGRAT was " << n << std::endl;
995 }
996 */
997 //
998 // begin by evaluating the next Pn from Eq 9.4:
999 //
1000 tnp1 += 2;
1001 p[n] = 0;
1002 T mbn = b - n;
1003 unsigned tmp1 = 3;
1004 for(unsigned m = 1; m < n; ++m)
1005 {
1006 mbn = m * b - n;
1007 p[n] += mbn * p[n-m] / boost::math::unchecked_factorial<T>(tmp1);
1008 tmp1 += 2;
1009 }
1010 p[n] /= n;
1011 p[n] += bm1 / boost::math::unchecked_factorial<T>(tnp1);
1012 //
1013 // Now we want Jn from Jn-1 using Eq 9.6:
1014 //
1015 j = (b2n * (b2n + 1) * j + (u + b2n + 1) * lxp) / t4;
1016 lxp *= lx2;
1017 b2n += 2;
1018 //
1019 // pull it together with Eq 9:
1020 //
1021 T r = prefix * p[n] * j;
1022 sum += r;
1023 // r is always small:
1024 BOOST_MATH_ASSERT(tools::max_value<T>() * tools::epsilon<T>() > fabs(r));
1025 if(fabs(r / tools::epsilon<T>()) < fabs(sum))
1026 break;
1027 }
1028 return sum;
1029} // template <class T, class Lanczos>T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const Lanczos& l, bool normalised)
1030
1031//
1032// For integer arguments we can relate the incomplete beta to the
1033// complement of the binomial distribution cdf and use this finite sum.
1034//
1035template <class T, class Policy>
1036T binomial_ccdf(T n, T k, T x, T y, const Policy& pol)
1037{
1038 BOOST_MATH_STD_USING // ADL of std names
1039
1040 T result = pow(x, n);
1041
1042 if(result > tools::min_value<T>())
1043 {
1044 T term = result;
1045 for(unsigned i = itrunc(T(n - 1)); i > k; --i)
1046 {
1047 term *= ((i + 1) * y) / ((n - i) * x);
1048 result += term;
1049 }
1050 }
1051 else
1052 {
1053 // First term underflows so we need to start at the mode of the
1054 // distribution and work outwards:
1055 int start = itrunc(n * x);
1056 if(start <= k + 1)
1057 start = itrunc(k + 2);
1058 result = static_cast<T>(pow(x, T(start)) * pow(y, n - T(start)) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(start), pol));
1059 if(result == 0)
1060 {
1061 // OK, starting slightly above the mode didn't work,
1062 // we'll have to sum the terms the old fashioned way.
1063 // Very hard to get here, possibly only when exponent
1064 // range is very limited (as with type float):
1065 // LCOV_EXCL_START
1066 for(unsigned i = start - 1; i > k; --i)
1067 {
1068 result += static_cast<T>(pow(x, static_cast<T>(i)) * pow(y, n - i) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(i), pol));
1069 }
1070 // LCOV_EXCL_STOP
1071 }
1072 else
1073 {
1074 T term = result;
1075 T start_term = result;
1076 for(unsigned i = start - 1; i > k; --i)
1077 {
1078 term *= ((i + 1) * y) / ((n - i) * x);
1079 result += term;
1080 }
1081 term = start_term;
1082 for(unsigned i = start + 1; i <= n; ++i)
1083 {
1084 term *= (n - i + 1) * x / (i * y);
1085 result += term;
1086 }
1087 }
1088 }
1089
1090 return result;
1091}
1092
1093
1094//
1095// The incomplete beta function implementation:
1096// This is just a big bunch of spaghetti code to divide up the
1097// input range and select the right implementation method for
1098// each domain:
1099//
1100template <class T, class Policy>
1101T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_derivative)
1102{
1103 static const char* function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)";
1104 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1105 BOOST_MATH_STD_USING // for ADL of std math functions.
1106
1107 BOOST_MATH_INSTRUMENT_VARIABLE(a);
1108 BOOST_MATH_INSTRUMENT_VARIABLE(b);
1109 BOOST_MATH_INSTRUMENT_VARIABLE(x);
1110 BOOST_MATH_INSTRUMENT_VARIABLE(inv);
1111 BOOST_MATH_INSTRUMENT_VARIABLE(normalised);
1112
1113 bool invert = inv;
1114 T fract;
1115 T y = 1 - x;
1116
1117 BOOST_MATH_ASSERT((p_derivative == 0) || normalised);
1118
1119 if(!(boost::math::isfinite)(a))
1120 return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be finite (got a=%1%).", a, pol);
1121 if(!(boost::math::isfinite)(b))
1122 return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be finite (got b=%1%).", b, pol);
1123 if (!(0 <= x && x <= 1))
1124 return policies::raise_domain_error<T>(function, "The argument x to the incomplete beta function must be in [0,1] (got x=%1%).", x, pol);
1125
1126 if(p_derivative)
1127 *p_derivative = -1; // value not set.
1128
1129 if(normalised)
1130 {
1131 if(a < 0)
1132 return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol);
1133 if(b < 0)
1134 return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol);
1135 // extend to a few very special cases:
1136 if(a == 0)
1137 {
1138 if(b == 0)
1139 return policies::raise_domain_error<T>(function, "The arguments a and b to the incomplete beta function cannot both be zero, with x=%1%.", x, pol);
1140 if(b > 0)
1141 return static_cast<T>(inv ? 0 : 1);
1142 }
1143 else if(b == 0)
1144 {
1145 if(a > 0)
1146 return static_cast<T>(inv ? 1 : 0);
1147 }
1148 }
1149 else
1150 {
1151 if(a <= 0)
1152 return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a, pol);
1153 if(b <= 0)
1154 return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be greater than zero (got b=%1%).", b, pol);
1155 }
1156
1157 if(x == 0)
1158 {
1159 if(p_derivative)
1160 {
1161 *p_derivative = (a == 1) ? (T)1 : (a < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2);
1162 }
1163 return (invert ? (normalised ? T(1) : boost::math::beta(a, b, pol)) : T(0));
1164 }
1165 if(x == 1)
1166 {
1167 if(p_derivative)
1168 {
1169 *p_derivative = (b == 1) ? T(1) : (b < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2);
1170 }
1171 return (invert == 0 ? (normalised ? 1 : boost::math::beta(a, b, pol)) : 0);
1172 }
1173 if((a == 0.5f) && (b == 0.5f))
1174 {
1175 // We have an arcsine distribution:
1176 if(p_derivative)
1177 {
1178 *p_derivative = 1 / (constants::pi<T>() * sqrt(y * x));
1179 }
1180 T p = invert ? asin(sqrt(y)) / constants::half_pi<T>() : asin(sqrt(x)) / constants::half_pi<T>();
1181 if(!normalised)
1182 p *= constants::pi<T>();
1183 return p;
1184 }
1185 if(a == 1)
1186 {
1187 std::swap(a, b);
1188 std::swap(x, y);
1189 invert = !invert;
1190 }
1191 if(b == 1)
1192 {
1193 //
1194 // Special case see: http://functions.wolfram.com/GammaBetaErf/BetaRegularized/03/01/01/
1195 //
1196 if(a == 1)
1197 {
1198 if(p_derivative)
1199 *p_derivative = 1;
1200 return invert ? y : x;
1201 }
1202
1203 if(p_derivative)
1204 {
1205 *p_derivative = a * pow(x, a - 1);
1206 }
1207 T p; // LCOV_EXCL_LINE
1208 if(y < 0.5)
1209 p = invert ? T(-boost::math::expm1(a * boost::math::log1p(-y, pol), pol)) : T(exp(a * boost::math::log1p(-y, pol)));
1210 else
1211 p = invert ? T(-boost::math::powm1(x, a, pol)) : T(pow(x, a));
1212 if(!normalised)
1213 p /= a;
1214 return p;
1215 }
1216
1217 if((std::min)(a, b) <= 1)
1218 {
1219 if(x > 0.5)
1220 {
1221 std::swap(a, b);
1222 std::swap(x, y);
1223 invert = !invert;
1224 BOOST_MATH_INSTRUMENT_VARIABLE(invert);
1225 }
1226 if((std::max)(a, b) <= 1)
1227 {
1228 // Both a,b < 1:
1229 if((a >= (std::min)(T(0.2), b)) || (pow(x, a) <= 0.9))
1230 {
1231 if(!invert)
1232 {
1233 fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
1234 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1235 }
1236 else
1237 {
1238 fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1239 invert = false;
1240 fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1241 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1242 }
1243 }
1244 else
1245 {
1246 std::swap(a, b);
1247 std::swap(x, y);
1248 invert = !invert;
1249 if(y >= 0.3)
1250 {
1251 if(!invert)
1252 {
1253 fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
1254 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1255 }
1256 else
1257 {
1258 fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1259 invert = false;
1260 fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1261 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1262 }
1263 }
1264 else
1265 {
1266 // Sidestep on a, and then use the series representation:
1267 T prefix; // LCOV_EXCL_LINE
1268 if(!normalised)
1269 {
1270 prefix = rising_factorial_ratio(T(a+b), a, 20);
1271 }
1272 else
1273 {
1274 prefix = 1;
1275 }
1276 fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative);
1277 if(!invert)
1278 {
1279 fract = beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
1280 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1281 }
1282 else
1283 {
1284 fract -= (normalised ? 1 : boost::math::beta(a, b, pol));
1285 invert = false;
1286 fract = -beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
1287 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1288 }
1289 }
1290 }
1291 }
1292 else
1293 {
1294 // One of a, b < 1 only:
1295 if((b <= 1) || ((x < 0.1) && (pow(b * x, a) <= 0.7)))
1296 {
1297 if(!invert)
1298 {
1299 fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
1300 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1301 }
1302 else
1303 {
1304 fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1305 invert = false;
1306 fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1307 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1308 }
1309 }
1310 else
1311 {
1312 std::swap(a, b);
1313 std::swap(x, y);
1314 invert = !invert;
1315
1316 if(y >= 0.3)
1317 {
1318 if(!invert)
1319 {
1320 fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
1321 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1322 }
1323 else
1324 {
1325 fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1326 invert = false;
1327 fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1328 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1329 }
1330 }
1331 else if(a >= 15)
1332 {
1333 if(!invert)
1334 {
1335 fract = beta_small_b_large_a_series(a, b, x, y, T(0), T(1), pol, normalised);
1336 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1337 }
1338 else
1339 {
1340 fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1341 invert = false;
1342 fract = -beta_small_b_large_a_series(a, b, x, y, fract, T(1), pol, normalised);
1343 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1344 }
1345 }
1346 else
1347 {
1348 // Sidestep to improve errors:
1349 T prefix; // LCOV_EXCL_LINE
1350 if(!normalised)
1351 {
1352 prefix = rising_factorial_ratio(T(a+b), a, 20);
1353 }
1354 else
1355 {
1356 prefix = 1;
1357 }
1358 fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative);
1359 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1360 if(!invert)
1361 {
1362 fract = beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
1363 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1364 }
1365 else
1366 {
1367 fract -= (normalised ? 1 : boost::math::beta(a, b, pol));
1368 invert = false;
1369 fract = -beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
1370 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1371 }
1372 }
1373 }
1374 }
1375 }
1376 else
1377 {
1378 // Both a,b >= 1:
1379 T lambda; // LCOV_EXCL_LINE
1380 if(a < b)
1381 {
1382 lambda = a - (a + b) * x;
1383 }
1384 else
1385 {
1386 lambda = (a + b) * y - b;
1387 }
1388 if(lambda < 0)
1389 {
1390 std::swap(a, b);
1391 std::swap(x, y);
1392 invert = !invert;
1393 BOOST_MATH_INSTRUMENT_VARIABLE(invert);
1394 }
1395
1396 if(b < 40)
1397 {
1398 if((floor(a) == a) && (floor(b) == b) && (a < static_cast<T>((std::numeric_limits<int>::max)() - 100)) && (y != 1))
1399 {
1400 // relate to the binomial distribution and use a finite sum:
1401 T k = a - 1;
1402 T n = b + k;
1403 fract = binomial_ccdf(n, k, x, y, pol);
1404 if(!normalised)
1405 fract *= boost::math::beta(a, b, pol);
1406 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1407 }
1408 else if(b * x <= 0.7)
1409 {
1410 if(!invert)
1411 {
1412 fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
1413 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1414 }
1415 else
1416 {
1417 fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1418 invert = false;
1419 fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1420 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1421 }
1422 }
1423 else if(a > 15)
1424 {
1425 // sidestep so we can use the series representation:
1426 int n = itrunc(T(floor(b)), pol);
1427 if(n == b)
1428 --n;
1429 T bbar = b - n;
1430 T prefix; // LCOV_EXCL_LINE
1431 if(!normalised)
1432 {
1433 prefix = rising_factorial_ratio(T(a+bbar), bbar, n);
1434 }
1435 else
1436 {
1437 prefix = 1;
1438 }
1439 fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(nullptr));
1440 fract = beta_small_b_large_a_series(a, bbar, x, y, fract, T(1), pol, normalised);
1441 fract /= prefix;
1442 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1443 }
1444 else if(normalised)
1445 {
1446 // The formula here for the non-normalised case is tricky to figure
1447 // out (for me!!), and requires two pochhammer calculations rather
1448 // than one, so leave it for now and only use this in the normalized case....
1449 int n = itrunc(T(floor(b)), pol);
1450 T bbar = b - n;
1451 if(bbar <= 0)
1452 {
1453 --n;
1454 bbar += 1;
1455 }
1456 fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(nullptr));
1457 fract += ibeta_a_step(a, bbar, x, y, 20, pol, normalised, static_cast<T*>(nullptr));
1458 if(invert)
1459 fract -= 1; // Note this line would need changing if we ever enable this branch in non-normalized case
1460 fract = beta_small_b_large_a_series(T(a+20), bbar, x, y, fract, T(1), pol, normalised);
1461 if(invert)
1462 {
1463 fract = -fract;
1464 invert = false;
1465 }
1466 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1467 }
1468 else
1469 {
1470 fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
1471 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1472 }
1473 }
1474 else
1475 {
1476 fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
1477 BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1478 }
1479 }
1480 if(p_derivative)
1481 {
1482 if(*p_derivative < 0)
1483 {
1484 *p_derivative = ibeta_power_terms(a, b, x, y, lanczos_type(), true, pol);
1485 }
1486 T div = y * x;
1487
1488 if(*p_derivative != 0)
1489 {
1490 if((tools::max_value<T>() * div < *p_derivative))
1491 {
1492 // overflow, return an arbitrarily large value:
1493 *p_derivative = tools::max_value<T>() / 2; // LCOV_EXCL_LINE Probably can only get here with denormalized x.
1494 }
1495 else
1496 {
1497 *p_derivative /= div;
1498 }
1499 }
1500 }
1501 return invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) - fract : fract;
1502} // template <class T, class Lanczos>T ibeta_imp(T a, T b, T x, const Lanczos& l, bool inv, bool normalised)
1503
1504template <class T, class Policy>
1505inline T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised)
1506{
1507 return ibeta_imp(a, b, x, pol, inv, normalised, static_cast<T*>(nullptr));
1508}
1509
1510template <class T, class Policy>
1511T ibeta_derivative_imp(T a, T b, T x, const Policy& pol)
1512{
1513 static const char* function = "ibeta_derivative<%1%>(%1%,%1%,%1%)";
1514 //
1515 // start with the usual error checks:
1516 //
1517 if (!(boost::math::isfinite)(a))
1518 return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be finite (got a=%1%).", a, pol);
1519 if (!(boost::math::isfinite)(b))
1520 return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be finite (got b=%1%).", b, pol);
1521 if (!(0 <= x && x <= 1))
1522 return policies::raise_domain_error<T>(function, "The argument x to the incomplete beta function must be in [0,1] (got x=%1%).", x, pol);
1523
1524 if(a <= 0)
1525 return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a, pol);
1526 if(b <= 0)
1527 return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be greater than zero (got b=%1%).", b, pol);
1528 //
1529 // Now the corner cases:
1530 //
1531 if(x == 0)
1532 {
1533 return (a > 1) ? 0 :
1534 (a == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error<T>(function, nullptr, pol);
1535 }
1536 else if(x == 1)
1537 {
1538 return (b > 1) ? 0 :
1539 (b == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error<T>(function, nullptr, pol);
1540 }
1541 //
1542 // Now the regular cases:
1543 //
1544 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1545 T y = (1 - x) * x;
1546 T f1;
1547 if (!(boost::math::isinf)(1 / y))
1548 {
1549 f1 = ibeta_power_terms<T>(a, b, x, 1 - x, lanczos_type(), true, pol, 1 / y, function);
1550 }
1551 else
1552 {
1553 return (a > 1) ? 0 : (a == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error<T>(function, nullptr, pol);
1554 }
1555
1556 return f1;
1557}
1558//
1559// Some forwarding functions that disambiguate the third argument type:
1560//
1561template <class RT1, class RT2, class Policy>
1562inline typename tools::promote_args<RT1, RT2>::type
1563 beta(RT1 a, RT2 b, const Policy&, const std::true_type*)
1564{
1565 BOOST_FPU_EXCEPTION_GUARD
1566 typedef typename tools::promote_args<RT1, RT2>::type result_type;
1567 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1568 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1569 typedef typename policies::normalise<
1570 Policy,
1571 policies::promote_float<false>,
1572 policies::promote_double<false>,
1573 policies::discrete_quantile<>,
1574 policies::assert_undefined<> >::type forwarding_policy;
1575
1576 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::beta_imp(static_cast<value_type>(a), static_cast<value_type>(b), evaluation_type(), forwarding_policy()), "boost::math::beta<%1%>(%1%,%1%)");
1577}
1578template <class RT1, class RT2, class RT3>
1579inline typename tools::promote_args<RT1, RT2, RT3>::type
1580 beta(RT1 a, RT2 b, RT3 x, const std::false_type*)
1581{
1582 return boost::math::beta(a, b, x, policies::policy<>());
1583}
1584} // namespace detail
1585
1586//
1587// The actual function entry-points now follow, these just figure out
1588// which Lanczos approximation to use
1589// and forward to the implementation functions:
1590//
1591template <class RT1, class RT2, class A>
1592inline typename tools::promote_args<RT1, RT2, A>::type
1593 beta(RT1 a, RT2 b, A arg)
1594{
1595 using tag = typename policies::is_policy<A>::type;
1596 using ReturnType = tools::promote_args_t<RT1, RT2, A>;
1597 return static_cast<ReturnType>(boost::math::detail::beta(a, b, arg, static_cast<tag*>(nullptr)));
1598}
1599
1600template <class RT1, class RT2>
1601inline typename tools::promote_args<RT1, RT2>::type
1602 beta(RT1 a, RT2 b)
1603{
1604 return boost::math::beta(a, b, policies::policy<>());
1605}
1606
1607template <class RT1, class RT2, class RT3, class Policy>
1608inline typename tools::promote_args<RT1, RT2, RT3>::type
1609 beta(RT1 a, RT2 b, RT3 x, const Policy&)
1610{
1611 BOOST_FPU_EXCEPTION_GUARD
1612 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1613 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1614 typedef typename policies::normalise<
1615 Policy,
1616 policies::promote_float<false>,
1617 policies::promote_double<false>,
1618 policies::discrete_quantile<>,
1619 policies::assert_undefined<> >::type forwarding_policy;
1620
1621 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), false, false), "boost::math::beta<%1%>(%1%,%1%,%1%)");
1622}
1623
1624template <class RT1, class RT2, class RT3, class Policy>
1625inline typename tools::promote_args<RT1, RT2, RT3>::type
1626 betac(RT1 a, RT2 b, RT3 x, const Policy&)
1627{
1628 BOOST_FPU_EXCEPTION_GUARD
1629 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1630 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1631 typedef typename policies::normalise<
1632 Policy,
1633 policies::promote_float<false>,
1634 policies::promote_double<false>,
1635 policies::discrete_quantile<>,
1636 policies::assert_undefined<> >::type forwarding_policy;
1637
1638 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), true, false), "boost::math::betac<%1%>(%1%,%1%,%1%)");
1639}
1640template <class RT1, class RT2, class RT3>
1641inline typename tools::promote_args<RT1, RT2, RT3>::type
1642 betac(RT1 a, RT2 b, RT3 x)
1643{
1644 return boost::math::betac(a, b, x, policies::policy<>());
1645}
1646
1647template <class RT1, class RT2, class RT3, class Policy>
1648inline typename tools::promote_args<RT1, RT2, RT3>::type
1649 ibeta(RT1 a, RT2 b, RT3 x, const Policy&)
1650{
1651 BOOST_FPU_EXCEPTION_GUARD
1652 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1653 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1654 typedef typename policies::normalise<
1655 Policy,
1656 policies::promote_float<false>,
1657 policies::promote_double<false>,
1658 policies::discrete_quantile<>,
1659 policies::assert_undefined<> >::type forwarding_policy;
1660
1661 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), false, true), "boost::math::ibeta<%1%>(%1%,%1%,%1%)");
1662}
1663template <class RT1, class RT2, class RT3>
1664inline typename tools::promote_args<RT1, RT2, RT3>::type
1665 ibeta(RT1 a, RT2 b, RT3 x)
1666{
1667 return boost::math::ibeta(a, b, x, policies::policy<>());
1668}
1669
1670template <class RT1, class RT2, class RT3, class Policy>
1671inline typename tools::promote_args<RT1, RT2, RT3>::type
1672 ibetac(RT1 a, RT2 b, RT3 x, const Policy&)
1673{
1674 BOOST_FPU_EXCEPTION_GUARD
1675 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1676 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1677 typedef typename policies::normalise<
1678 Policy,
1679 policies::promote_float<false>,
1680 policies::promote_double<false>,
1681 policies::discrete_quantile<>,
1682 policies::assert_undefined<> >::type forwarding_policy;
1683
1684 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), true, true), "boost::math::ibetac<%1%>(%1%,%1%,%1%)");
1685}
1686template <class RT1, class RT2, class RT3>
1687inline typename tools::promote_args<RT1, RT2, RT3>::type
1688 ibetac(RT1 a, RT2 b, RT3 x)
1689{
1690 return boost::math::ibetac(a, b, x, policies::policy<>());
1691}
1692
1693template <class RT1, class RT2, class RT3, class Policy>
1694inline typename tools::promote_args<RT1, RT2, RT3>::type
1695 ibeta_derivative(RT1 a, RT2 b, RT3 x, const Policy&)
1696{
1697 BOOST_FPU_EXCEPTION_GUARD
1698 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1699 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1700 typedef typename policies::normalise<
1701 Policy,
1702 policies::promote_float<false>,
1703 policies::promote_double<false>,
1704 policies::discrete_quantile<>,
1705 policies::assert_undefined<> >::type forwarding_policy;
1706
1707 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_derivative_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy()), "boost::math::ibeta_derivative<%1%>(%1%,%1%,%1%)");
1708}
1709template <class RT1, class RT2, class RT3>
1710inline typename tools::promote_args<RT1, RT2, RT3>::type
1711 ibeta_derivative(RT1 a, RT2 b, RT3 x)
1712{
1713 return boost::math::ibeta_derivative(a, b, x, policies::policy<>());
1714}
1715
1716} // namespace math
1717} // namespace boost
1718
1719#include <boost/math/special_functions/detail/ibeta_inverse.hpp>
1720#include <boost/math/special_functions/detail/ibeta_inv_ab.hpp>
1721
1722#endif // BOOST_MATH_SPECIAL_BETA_HPP
1723

source code of boost/libs/math/include/boost/math/special_functions/beta.hpp