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/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
  26 namespace boost{ namespace math{
  27
  28 namespace detail{
  29
  30 //
  31 // Implementation of Beta(a,b) using the Lanczos approximation:
  32 //
  33 template <class T, class Lanczos, class Policy>
  34 T 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 //
 123 template <class T, class Policy>
 124 T 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 //
 206 template <class T, class Lanczos, class Policy>
 207 T 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 //
 434 template <class T, class Policy>
 435 T 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 //
 526 template <class T>
 527 struct 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    }
 540 private:
 541    T result, x, apn, poch;
 542    int n;
 543 };
 544
 545 template <class T, class Lanczos, class Policy>
 546 T 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 //
 614 template <class T, class Policy>
 615 T 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 //
 676 template <class T>
 677 struct 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
 698 private:
 699    T a, b, x, y;
 700    int m;
 701 };
 702 //
 703 // Evaluate the incomplete beta via the continued fraction representation:
 704 //
 705 template <class T, class Policy>
 706 inline 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 //
 728 template <class T, class Policy>
 729 T 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 //
 762 template <class T>
 763 inline 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 "guestimates" based on empiracal observation.
 786 // Note that the table size should never exceed the size of our
 787 // tables of factorials.
 788 //
 789 template <class T>
 790 struct 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 };
 800 template <>
 801 struct 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 };
 806 template <>
 807 struct 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 };
 812 template <>
 813 struct 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
 819 template <class T, class Policy>
 820 T 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, unfortunatately 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 //
 931 template <class T>
 932 T 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 spagetti code to divide up the
 989 // input range and select the right implementation method for
 990 // each domain:
 991 //
 992 template <class T, class Policy>
 993 T 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 arbitarily 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
1392 template <class T, class Policy>
1393 inline 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
1398 template <class T, class Policy>
1399 T 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 //
1435 template <class RT1, class RT2, class Policy>
1436 inline typename tools::promote_args<RT1, RT2>::type
1437    beta(RT1 a, RT2 b, const Policy&, const mpl::true_*)
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 }
1452 template <class RT1, class RT2, class RT3>
1453 inline typename tools::promote_args<RT1, RT2, RT3>::type
1454    beta(RT1 a, RT2 b, RT3 x, const mpl::false_*)
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 //
1465 template <class RT1, class RT2, class A>
1466 inline 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
1473 template <class RT1, class RT2>
1474 inline 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
1480 template <class RT1, class RT2, class RT3, class Policy>
1481 inline 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
1497 template <class RT1, class RT2, class RT3, class Policy>
1498 inline 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 }
1513 template <class RT1, class RT2, class RT3>
1514 inline 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
1520 template <class RT1, class RT2, class RT3, class Policy>
1521 inline 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 }
1536 template <class RT1, class RT2, class RT3>
1537 inline 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
1543 template <class RT1, class RT2, class RT3, class Policy>
1544 inline 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 }
1559 template <class RT1, class RT2, class RT3>
1560 inline 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
1566 template <class RT1, class RT2, class RT3, class Policy>
1567 inline 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 }
1582 template <class RT1, class RT2, class RT3>
1583 inline 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