boost/math/special_functions/beta.hpp

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