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

Provided by KDAB

Privacy Policy
Update your C++ knowledge – Modern C++11/14/17 Training
Find out more

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