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

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