Chris@16: //  (C) Copyright John Maddock 2006.
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_MATH_SPECIAL_BETA_HPP
Chris@16: #define BOOST_MATH_SPECIAL_BETA_HPP
Chris@16: 
Chris@16: #ifdef _MSC_VER
Chris@16: #pragma once
Chris@16: #endif
Chris@16: 
Chris@16: #include <boost/math/special_functions/math_fwd.hpp>
Chris@16: #include <boost/math/tools/config.hpp>
Chris@16: #include <boost/math/special_functions/gamma.hpp>
Chris@101: #include <boost/math/special_functions/binomial.hpp>
Chris@16: #include <boost/math/special_functions/factorials.hpp>
Chris@16: #include <boost/math/special_functions/erf.hpp>
Chris@16: #include <boost/math/special_functions/log1p.hpp>
Chris@16: #include <boost/math/special_functions/expm1.hpp>
Chris@16: #include <boost/math/special_functions/trunc.hpp>
Chris@16: #include <boost/math/tools/roots.hpp>
Chris@16: #include <boost/static_assert.hpp>
Chris@16: #include <boost/config/no_tr1/cmath.hpp>
Chris@16: 
Chris@16: namespace boost{ namespace math{
Chris@16: 
Chris@16: namespace detail{
Chris@16: 
Chris@16: //
Chris@16: // Implementation of Beta(a,b) using the Lanczos approximation:
Chris@16: //
Chris@16: template <class T, class Lanczos, class Policy>
Chris@16: T beta_imp(T a, T b, const Lanczos&, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING  // for ADL of std names
Chris@16: 
Chris@16:    if(a <= 0)
Chris@101:       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);
Chris@16:    if(b <= 0)
Chris@101:       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);
Chris@16: 
Chris@16:    T result;
Chris@16: 
Chris@16:    T prefix = 1;
Chris@16:    T c = a + b;
Chris@16: 
Chris@16:    // Special cases:
Chris@16:    if((c == a) && (b < tools::epsilon<T>()))
Chris@16:       return boost::math::tgamma(b, pol);
Chris@16:    else if((c == b) && (a < tools::epsilon<T>()))
Chris@16:       return boost::math::tgamma(a, pol);
Chris@16:    if(b == 1)
Chris@16:       return 1/a;
Chris@16:    else if(a == 1)
Chris@16:       return 1/b;
Chris@16: 
Chris@16:    /*
Chris@16:    //
Chris@16:    // This code appears to be no longer necessary: it was
Chris@16:    // used to offset errors introduced from the Lanczos
Chris@16:    // approximation, but the current Lanczos approximations
Chris@16:    // are sufficiently accurate for all z that we can ditch
Chris@16:    // this.  It remains in the file for future reference...
Chris@16:    //
Chris@16:    // If a or b are less than 1, shift to greater than 1:
Chris@16:    if(a < 1)
Chris@16:    {
Chris@16:       prefix *= c / a;
Chris@16:       c += 1;
Chris@16:       a += 1;
Chris@16:    }
Chris@16:    if(b < 1)
Chris@16:    {
Chris@16:       prefix *= c / b;
Chris@16:       c += 1;
Chris@16:       b += 1;
Chris@16:    }
Chris@16:    */
Chris@16: 
Chris@16:    if(a < b)
Chris@16:       std::swap(a, b);
Chris@16: 
Chris@16:    // Lanczos calculation:
Chris@16:    T agh = a + Lanczos::g() - T(0.5);
Chris@16:    T bgh = b + Lanczos::g() - T(0.5);
Chris@16:    T cgh = c + Lanczos::g() - T(0.5);
Chris@16:    result = Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b) / Lanczos::lanczos_sum_expG_scaled(c);
Chris@16:    T ambh = a - T(0.5) - b;
Chris@16:    if((fabs(b * ambh) < (cgh * 100)) && (a > 100))
Chris@16:    {
Chris@16:       // Special case where the base of the power term is close to 1
Chris@16:       // compute (1+x)^y instead:
Chris@16:       result *= exp(ambh * boost::math::log1p(-b / cgh, pol));
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       result *= pow(agh / cgh, a - T(0.5) - b);
Chris@16:    }
Chris@16:    if(cgh > 1e10f)
Chris@16:       // this avoids possible overflow, but appears to be marginally less accurate:
Chris@16:       result *= pow((agh / cgh) * (bgh / cgh), b);
Chris@16:    else
Chris@16:       result *= pow((agh * bgh) / (cgh * cgh), b);
Chris@16:    result *= sqrt(boost::math::constants::e<T>() / bgh);
Chris@16: 
Chris@16:    // If a and b were originally less than 1 we need to scale the result:
Chris@16:    result *= prefix;
Chris@16: 
Chris@16:    return result;
Chris@16: } // template <class T, class Lanczos> beta_imp(T a, T b, const Lanczos&)
Chris@16: 
Chris@16: //
Chris@16: // Generic implementation of Beta(a,b) without Lanczos approximation support
Chris@16: // (Caution this is slow!!!):
Chris@16: //
Chris@16: template <class T, class Policy>
Chris@16: T beta_imp(T a, T b, const lanczos::undefined_lanczos& /* l */, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    if(a <= 0)
Chris@101:       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);
Chris@16:    if(b <= 0)
Chris@101:       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);
Chris@16: 
Chris@16:    T result;
Chris@16: 
Chris@16:    T prefix = 1;
Chris@16:    T c = a + b;
Chris@16: 
Chris@16:    // special cases:
Chris@16:    if((c == a) && (b < tools::epsilon<T>()))
Chris@16:       return boost::math::tgamma(b, pol);
Chris@16:    else if((c == b) && (a < tools::epsilon<T>()))
Chris@16:       return boost::math::tgamma(a, pol);
Chris@16:    if(b == 1)
Chris@16:       return 1/a;
Chris@16:    else if(a == 1)
Chris@16:       return 1/b;
Chris@16: 
Chris@16:    // shift to a and b > 1 if required:
Chris@16:    if(a < 1)
Chris@16:    {
Chris@16:       prefix *= c / a;
Chris@16:       c += 1;
Chris@16:       a += 1;
Chris@16:    }
Chris@16:    if(b < 1)
Chris@16:    {
Chris@16:       prefix *= c / b;
Chris@16:       c += 1;
Chris@16:       b += 1;
Chris@16:    }
Chris@16:    if(a < b)
Chris@16:       std::swap(a, b);
Chris@16: 
Chris@16:    // set integration limits:
Chris@16:    T la = (std::max)(T(10), a);
Chris@16:    T lb = (std::max)(T(10), b);
Chris@16:    T lc = (std::max)(T(10), T(a+b));
Chris@16: 
Chris@16:    // calculate the fraction parts:
Chris@16:    T sa = detail::lower_gamma_series(a, la, pol) / a;
Chris@16:    sa += detail::upper_gamma_fraction(a, la, ::boost::math::policies::get_epsilon<T, Policy>());
Chris@16:    T sb = detail::lower_gamma_series(b, lb, pol) / b;
Chris@16:    sb += detail::upper_gamma_fraction(b, lb, ::boost::math::policies::get_epsilon<T, Policy>());
Chris@16:    T sc = detail::lower_gamma_series(c, lc, pol) / c;
Chris@16:    sc += detail::upper_gamma_fraction(c, lc, ::boost::math::policies::get_epsilon<T, Policy>());
Chris@16: 
Chris@16:    // and the exponent part:
Chris@16:    result = exp(lc - la - lb) * pow(la/lc, a) * pow(lb/lc, b);
Chris@16: 
Chris@16:    // and combine:
Chris@16:    result *= sa * sb / sc;
Chris@16: 
Chris@16:    // if a and b were originally less than 1 we need to scale the result:
Chris@16:    result *= prefix;
Chris@16: 
Chris@16:    return result;
Chris@16: } // template <class T>T beta_imp(T a, T b, const lanczos::undefined_lanczos& l)
Chris@16: 
Chris@16: 
Chris@16: //
Chris@16: // Compute the leading power terms in the incomplete Beta:
Chris@16: //
Chris@16: // (x^a)(y^b)/Beta(a,b) when normalised, and
Chris@16: // (x^a)(y^b) otherwise.
Chris@16: //
Chris@16: // Almost all of the error in the incomplete beta comes from this
Chris@16: // function: particularly when a and b are large. Computing large
Chris@16: // powers are *hard* though, and using logarithms just leads to
Chris@16: // horrendous cancellation errors.
Chris@16: //
Chris@16: template <class T, class Lanczos, class Policy>
Chris@16: T ibeta_power_terms(T a,
Chris@16:                         T b,
Chris@16:                         T x,
Chris@16:                         T y,
Chris@16:                         const Lanczos&,
Chris@16:                         bool normalised,
Chris@16:                         const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    if(!normalised)
Chris@16:    {
Chris@16:       // can we do better here?
Chris@16:       return pow(x, a) * pow(y, b);
Chris@16:    }
Chris@16: 
Chris@16:    T result;
Chris@16: 
Chris@16:    T prefix = 1;
Chris@16:    T c = a + b;
Chris@16: 
Chris@16:    // combine power terms with Lanczos approximation:
Chris@16:    T agh = a + Lanczos::g() - T(0.5);
Chris@16:    T bgh = b + Lanczos::g() - T(0.5);
Chris@16:    T cgh = c + Lanczos::g() - T(0.5);
Chris@16:    result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));
Chris@16: 
Chris@16:    // l1 and l2 are the base of the exponents minus one:
Chris@16:    T l1 = (x * b - y * agh) / agh;
Chris@16:    T l2 = (y * a - x * bgh) / bgh;
Chris@16:    if(((std::min)(fabs(l1), fabs(l2)) < 0.2))
Chris@16:    {
Chris@16:       // when the base of the exponent is very near 1 we get really
Chris@16:       // gross errors unless extra care is taken:
Chris@16:       if((l1 * l2 > 0) || ((std::min)(a, b) < 1))
Chris@16:       {
Chris@16:          //
Chris@16:          // This first branch handles the simple cases where either: 
Chris@16:          //
Chris@16:          // * The two power terms both go in the same direction 
Chris@16:          // (towards zero or towards infinity).  In this case if either 
Chris@16:          // term overflows or underflows, then the product of the two must 
Chris@16:          // do so also.  
Chris@16:          // *Alternatively if one exponent is less than one, then we 
Chris@16:          // can't productively use it to eliminate overflow or underflow 
Chris@16:          // from the other term.  Problems with spurious overflow/underflow 
Chris@16:          // can't be ruled out in this case, but it is *very* unlikely 
Chris@16:          // since one of the power terms will evaluate to a number close to 1.
Chris@16:          //
Chris@16:          if(fabs(l1) < 0.1)
Chris@16:          {
Chris@16:             result *= exp(a * boost::math::log1p(l1, pol));
Chris@16:             BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:          }
Chris@16:          else
Chris@16:          {
Chris@16:             result *= pow((x * cgh) / agh, a);
Chris@16:             BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:          }
Chris@16:          if(fabs(l2) < 0.1)
Chris@16:          {
Chris@16:             result *= exp(b * boost::math::log1p(l2, pol));
Chris@16:             BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:          }
Chris@16:          else
Chris@16:          {
Chris@16:             result *= pow((y * cgh) / bgh, b);
Chris@16:             BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:          }
Chris@16:       }
Chris@16:       else if((std::max)(fabs(l1), fabs(l2)) < 0.5)
Chris@16:       {
Chris@16:          //
Chris@16:          // Both exponents are near one and both the exponents are 
Chris@16:          // greater than one and further these two 
Chris@16:          // power terms tend in opposite directions (one towards zero, 
Chris@16:          // the other towards infinity), so we have to combine the terms 
Chris@16:          // to avoid any risk of overflow or underflow.
Chris@16:          //
Chris@16:          // We do this by moving one power term inside the other, we have:
Chris@16:          //
Chris@16:          //    (1 + l1)^a * (1 + l2)^b
Chris@16:          //  = ((1 + l1)*(1 + l2)^(b/a))^a
Chris@16:          //  = (1 + l1 + l3 + l1*l3)^a   ;  l3 = (1 + l2)^(b/a) - 1
Chris@16:          //                                    = exp((b/a) * log(1 + l2)) - 1
Chris@16:          //
Chris@16:          // The tricky bit is deciding which term to move inside :-)
Chris@16:          // By preference we move the larger term inside, so that the
Chris@16:          // size of the largest exponent is reduced.  However, that can
Chris@16:          // only be done as long as l3 (see above) is also small.
Chris@16:          //
Chris@16:          bool small_a = a < b;
Chris@16:          T ratio = b / a;
Chris@16:          if((small_a && (ratio * l2 < 0.1)) || (!small_a && (l1 / ratio > 0.1)))
Chris@16:          {
Chris@16:             T l3 = boost::math::expm1(ratio * boost::math::log1p(l2, pol), pol);
Chris@16:             l3 = l1 + l3 + l3 * l1;
Chris@16:             l3 = a * boost::math::log1p(l3, pol);
Chris@16:             result *= exp(l3);
Chris@16:             BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:          }
Chris@16:          else
Chris@16:          {
Chris@16:             T l3 = boost::math::expm1(boost::math::log1p(l1, pol) / ratio, pol);
Chris@16:             l3 = l2 + l3 + l3 * l2;
Chris@16:             l3 = b * boost::math::log1p(l3, pol);
Chris@16:             result *= exp(l3);
Chris@16:             BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:          }
Chris@16:       }
Chris@16:       else if(fabs(l1) < fabs(l2))
Chris@16:       {
Chris@16:          // First base near 1 only:
Chris@16:          T l = a * boost::math::log1p(l1, pol)
Chris@16:             + b * log((y * cgh) / bgh);
Chris@16:          result *= exp(l);
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          // Second base near 1 only:
Chris@16:          T l = b * boost::math::log1p(l2, pol)
Chris@16:             + a * log((x * cgh) / agh);
Chris@16:          result *= exp(l);
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:       }
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       // general case:
Chris@16:       T b1 = (x * cgh) / agh;
Chris@16:       T b2 = (y * cgh) / bgh;
Chris@16:       l1 = a * log(b1);
Chris@16:       l2 = b * log(b2);
Chris@16:       BOOST_MATH_INSTRUMENT_VARIABLE(b1);
Chris@16:       BOOST_MATH_INSTRUMENT_VARIABLE(b2);
Chris@16:       BOOST_MATH_INSTRUMENT_VARIABLE(l1);
Chris@16:       BOOST_MATH_INSTRUMENT_VARIABLE(l2);
Chris@16:       if((l1 >= tools::log_max_value<T>())
Chris@16:          || (l1 <= tools::log_min_value<T>())
Chris@16:          || (l2 >= tools::log_max_value<T>())
Chris@16:          || (l2 <= tools::log_min_value<T>())
Chris@16:          )
Chris@16:       {
Chris@16:          // Oops, overflow, sidestep:
Chris@16:          if(a < b)
Chris@16:             result *= pow(pow(b2, b/a) * b1, a);
Chris@16:          else
Chris@16:             result *= pow(pow(b1, a/b) * b2, b);
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          // finally the normal case:
Chris@16:          result *= pow(b1, a) * pow(b2, b);
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:       }
Chris@16:    }
Chris@16:    // combine with the leftover terms from the Lanczos approximation:
Chris@16:    result *= sqrt(bgh / boost::math::constants::e<T>());
Chris@16:    result *= sqrt(agh / cgh);
Chris@16:    result *= prefix;
Chris@16: 
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16: 
Chris@16:    return result;
Chris@16: }
Chris@16: //
Chris@16: // Compute the leading power terms in the incomplete Beta:
Chris@16: //
Chris@16: // (x^a)(y^b)/Beta(a,b) when normalised, and
Chris@16: // (x^a)(y^b) otherwise.
Chris@16: //
Chris@16: // Almost all of the error in the incomplete beta comes from this
Chris@16: // function: particularly when a and b are large. Computing large
Chris@16: // powers are *hard* though, and using logarithms just leads to
Chris@16: // horrendous cancellation errors.
Chris@16: //
Chris@16: // This version is generic, slow, and does not use the Lanczos approximation.
Chris@16: //
Chris@16: template <class T, class Policy>
Chris@16: T ibeta_power_terms(T a,
Chris@16:                         T b,
Chris@16:                         T x,
Chris@16:                         T y,
Chris@16:                         const boost::math::lanczos::undefined_lanczos&,
Chris@16:                         bool normalised,
Chris@16:                         const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    if(!normalised)
Chris@16:    {
Chris@16:       return pow(x, a) * pow(y, b);
Chris@16:    }
Chris@16: 
Chris@16:    T result= 0; // assignment here silences warnings later
Chris@16: 
Chris@16:    T c = a + b;
Chris@16: 
Chris@16:    // integration limits for the gamma functions:
Chris@16:    //T la = (std::max)(T(10), a);
Chris@16:    //T lb = (std::max)(T(10), b);
Chris@16:    //T lc = (std::max)(T(10), a+b);
Chris@16:    T la = a + 5;
Chris@16:    T lb = b + 5;
Chris@16:    T lc = a + b + 5;
Chris@16:    // gamma function partials:
Chris@16:    T sa = detail::lower_gamma_series(a, la, pol) / a;
Chris@16:    sa += detail::upper_gamma_fraction(a, la, ::boost::math::policies::get_epsilon<T, Policy>());
Chris@16:    T sb = detail::lower_gamma_series(b, lb, pol) / b;
Chris@16:    sb += detail::upper_gamma_fraction(b, lb, ::boost::math::policies::get_epsilon<T, Policy>());
Chris@16:    T sc = detail::lower_gamma_series(c, lc, pol) / c;
Chris@16:    sc += detail::upper_gamma_fraction(c, lc, ::boost::math::policies::get_epsilon<T, Policy>());
Chris@16:    // gamma function powers combined with incomplete beta powers:
Chris@16: 
Chris@16:    T b1 = (x * lc) / la;
Chris@16:    T b2 = (y * lc) / lb;
Chris@16:    T e1 = lc - la - lb;
Chris@16:    T lb1 = a * log(b1);
Chris@16:    T lb2 = b * log(b2);
Chris@16: 
Chris@16:    if((lb1 >= tools::log_max_value<T>())
Chris@16:       || (lb1 <= tools::log_min_value<T>())
Chris@16:       || (lb2 >= tools::log_max_value<T>())
Chris@16:       || (lb2 <= tools::log_min_value<T>())
Chris@16:       || (e1 >= tools::log_max_value<T>())
Chris@16:       || (e1 <= tools::log_min_value<T>())
Chris@16:       )
Chris@16:    {
Chris@16:       result = exp(lb1 + lb2 - e1);
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       T p1, p2;
Chris@16:       if((fabs(b1 - 1) * a < 10) && (a > 1))
Chris@16:          p1 = exp(a * boost::math::log1p((x * b - y * la) / la, pol));
Chris@16:       else
Chris@16:          p1 = pow(b1, a);
Chris@16:       if((fabs(b2 - 1) * b < 10) && (b > 1))
Chris@16:          p2 = exp(b * boost::math::log1p((y * a - x * lb) / lb, pol));
Chris@16:       else
Chris@16:          p2 = pow(b2, b);
Chris@16:       T p3 = exp(e1);
Chris@16:       result = p1 * p2 / p3;
Chris@16:    }
Chris@16:    // and combine with the remaining gamma function components:
Chris@16:    result /= sa * sb / sc;
Chris@16: 
Chris@16:    return result;
Chris@16: }
Chris@16: //
Chris@16: // Series approximation to the incomplete beta:
Chris@16: //
Chris@16: template <class T>
Chris@16: struct ibeta_series_t
Chris@16: {
Chris@16:    typedef T result_type;
Chris@16:    ibeta_series_t(T a_, T b_, T x_, T mult) : result(mult), x(x_), apn(a_), poch(1-b_), n(1) {}
Chris@16:    T operator()()
Chris@16:    {
Chris@16:       T r = result / apn;
Chris@16:       apn += 1;
Chris@16:       result *= poch * x / n;
Chris@16:       ++n;
Chris@16:       poch += 1;
Chris@16:       return r;
Chris@16:    }
Chris@16: private:
Chris@16:    T result, x, apn, poch;
Chris@16:    int n;
Chris@16: };
Chris@16: 
Chris@16: template <class T, class Lanczos, class Policy>
Chris@16: T ibeta_series(T a, T b, T x, T s0, const Lanczos&, bool normalised, T* p_derivative, T y, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    T result;
Chris@16: 
Chris@16:    BOOST_ASSERT((p_derivative == 0) || normalised);
Chris@16: 
Chris@16:    if(normalised)
Chris@16:    {
Chris@16:       T c = a + b;
Chris@16: 
Chris@16:       // incomplete beta power term, combined with the Lanczos approximation:
Chris@16:       T agh = a + Lanczos::g() - T(0.5);
Chris@16:       T bgh = b + Lanczos::g() - T(0.5);
Chris@16:       T cgh = c + Lanczos::g() - T(0.5);
Chris@16:       result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));
Chris@16:       if(a * b < bgh * 10)
Chris@16:          result *= exp((b - 0.5f) * boost::math::log1p(a / bgh, pol));
Chris@16:       else
Chris@16:          result *= pow(cgh / bgh, b - 0.5f);
Chris@16:       result *= pow(x * cgh / agh, a);
Chris@16:       result *= sqrt(agh / boost::math::constants::e<T>());
Chris@16: 
Chris@16:       if(p_derivative)
Chris@16:       {
Chris@16:          *p_derivative = result * pow(y, b);
Chris@16:          BOOST_ASSERT(*p_derivative >= 0);
Chris@16:       }
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       // Non-normalised, just compute the power:
Chris@16:       result = pow(x, a);
Chris@16:    }
Chris@16:    if(result < tools::min_value<T>())
Chris@16:       return s0; // Safeguard: series can't cope with denorms.
Chris@16:    ibeta_series_t<T> s(a, b, x, result);
Chris@16:    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
Chris@16:    result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0);
Chris@16:    policies::check_series_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (with lanczos)", max_iter, pol);
Chris@16:    return result;
Chris@16: }
Chris@16: //
Chris@16: // Incomplete Beta series again, this time without Lanczos support:
Chris@16: //
Chris@16: template <class T, class Policy>
Chris@16: 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)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    T result;
Chris@16:    BOOST_ASSERT((p_derivative == 0) || normalised);
Chris@16: 
Chris@16:    if(normalised)
Chris@16:    {
Chris@16:       T c = a + b;
Chris@16: 
Chris@16:       // figure out integration limits for the gamma function:
Chris@16:       //T la = (std::max)(T(10), a);
Chris@16:       //T lb = (std::max)(T(10), b);
Chris@16:       //T lc = (std::max)(T(10), a+b);
Chris@16:       T la = a + 5;
Chris@16:       T lb = b + 5;
Chris@16:       T lc = a + b + 5;
Chris@16: 
Chris@16:       // calculate the gamma parts:
Chris@16:       T sa = detail::lower_gamma_series(a, la, pol) / a;
Chris@16:       sa += detail::upper_gamma_fraction(a, la, ::boost::math::policies::get_epsilon<T, Policy>());
Chris@16:       T sb = detail::lower_gamma_series(b, lb, pol) / b;
Chris@16:       sb += detail::upper_gamma_fraction(b, lb, ::boost::math::policies::get_epsilon<T, Policy>());
Chris@16:       T sc = detail::lower_gamma_series(c, lc, pol) / c;
Chris@16:       sc += detail::upper_gamma_fraction(c, lc, ::boost::math::policies::get_epsilon<T, Policy>());
Chris@16: 
Chris@16:       // and their combined power-terms:
Chris@16:       T b1 = (x * lc) / la;
Chris@16:       T b2 = lc/lb;
Chris@16:       T e1 = lc - la - lb;
Chris@16:       T lb1 = a * log(b1);
Chris@16:       T lb2 = b * log(b2);
Chris@16: 
Chris@16:       if((lb1 >= tools::log_max_value<T>())
Chris@16:          || (lb1 <= tools::log_min_value<T>())
Chris@16:          || (lb2 >= tools::log_max_value<T>())
Chris@16:          || (lb2 <= tools::log_min_value<T>())
Chris@16:          || (e1 >= tools::log_max_value<T>())
Chris@16:          || (e1 <= tools::log_min_value<T>()) )
Chris@16:       {
Chris@16:          T p = lb1 + lb2 - e1;
Chris@16:          result = exp(p);
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          result = pow(b1, a);
Chris@16:          if(a * b < lb * 10)
Chris@16:             result *= exp(b * boost::math::log1p(a / lb, pol));
Chris@16:          else
Chris@16:             result *= pow(b2, b);
Chris@16:          result /= exp(e1);
Chris@16:       }
Chris@16:       // and combine the results:
Chris@16:       result /= sa * sb / sc;
Chris@16: 
Chris@16:       if(p_derivative)
Chris@16:       {
Chris@16:          *p_derivative = result * pow(y, b);
Chris@16:          BOOST_ASSERT(*p_derivative >= 0);
Chris@16:       }
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       // Non-normalised, just compute the power:
Chris@16:       result = pow(x, a);
Chris@16:    }
Chris@16:    if(result < tools::min_value<T>())
Chris@16:       return s0; // Safeguard: series can't cope with denorms.
Chris@16:    ibeta_series_t<T> s(a, b, x, result);
Chris@16:    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
Chris@16:    result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0);
Chris@16:    policies::check_series_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (without lanczos)", max_iter, pol);
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: //
Chris@16: // Continued fraction for the incomplete beta:
Chris@16: //
Chris@16: template <class T>
Chris@16: struct ibeta_fraction2_t
Chris@16: {
Chris@16:    typedef std::pair<T, T> result_type;
Chris@16: 
Chris@16:    ibeta_fraction2_t(T a_, T b_, T x_, T y_) : a(a_), b(b_), x(x_), y(y_), m(0) {}
Chris@16: 
Chris@16:    result_type operator()()
Chris@16:    {
Chris@16:       T aN = (a + m - 1) * (a + b + m - 1) * m * (b - m) * x * x;
Chris@16:       T denom = (a + 2 * m - 1);
Chris@16:       aN /= denom * denom;
Chris@16: 
Chris@16:       T bN = m;
Chris@16:       bN += (m * (b - m) * x) / (a + 2*m - 1);
Chris@16:       bN += ((a + m) * (a * y - b * x + 1 + m *(2 - x))) / (a + 2*m + 1);
Chris@16: 
Chris@16:       ++m;
Chris@16: 
Chris@16:       return std::make_pair(aN, bN);
Chris@16:    }
Chris@16: 
Chris@16: private:
Chris@16:    T a, b, x, y;
Chris@16:    int m;
Chris@16: };
Chris@16: //
Chris@16: // Evaluate the incomplete beta via the continued fraction representation:
Chris@16: //
Chris@16: template <class T, class Policy>
Chris@16: inline T ibeta_fraction2(T a, T b, T x, T y, const Policy& pol, bool normalised, T* p_derivative)
Chris@16: {
Chris@16:    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    T result = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol);
Chris@16:    if(p_derivative)
Chris@16:    {
Chris@16:       *p_derivative = result;
Chris@16:       BOOST_ASSERT(*p_derivative >= 0);
Chris@16:    }
Chris@16:    if(result == 0)
Chris@16:       return result;
Chris@16: 
Chris@16:    ibeta_fraction2_t<T> f(a, b, x, y);
Chris@16:    T fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon<T, Policy>());
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:    return result / fract;
Chris@16: }
Chris@16: //
Chris@16: // Computes the difference between ibeta(a,b,x) and ibeta(a+k,b,x):
Chris@16: //
Chris@16: template <class T, class Policy>
Chris@16: T ibeta_a_step(T a, T b, T x, T y, int k, const Policy& pol, bool normalised, T* p_derivative)
Chris@16: {
Chris@16:    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
Chris@16: 
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(k);
Chris@16: 
Chris@16:    T prefix = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol);
Chris@16:    if(p_derivative)
Chris@16:    {
Chris@16:       *p_derivative = prefix;
Chris@16:       BOOST_ASSERT(*p_derivative >= 0);
Chris@16:    }
Chris@16:    prefix /= a;
Chris@16:    if(prefix == 0)
Chris@16:       return prefix;
Chris@16:    T sum = 1;
Chris@16:    T term = 1;
Chris@16:    // series summation from 0 to k-1:
Chris@16:    for(int i = 0; i < k-1; ++i)
Chris@16:    {
Chris@16:       term *= (a+b+i) * x / (a+i+1);
Chris@16:       sum += term;
Chris@16:    }
Chris@16:    prefix *= sum;
Chris@16: 
Chris@16:    return prefix;
Chris@16: }
Chris@16: //
Chris@16: // This function is only needed for the non-regular incomplete beta,
Chris@16: // it computes the delta in:
Chris@16: // beta(a,b,x) = prefix + delta * beta(a+k,b,x)
Chris@16: // it is currently only called for small k.
Chris@16: //
Chris@16: template <class T>
Chris@16: inline T rising_factorial_ratio(T a, T b, int k)
Chris@16: {
Chris@16:    // calculate:
Chris@16:    // (a)(a+1)(a+2)...(a+k-1)
Chris@16:    // _______________________
Chris@16:    // (b)(b+1)(b+2)...(b+k-1)
Chris@16: 
Chris@16:    // This is only called with small k, for large k
Chris@16:    // it is grossly inefficient, do not use outside it's
Chris@16:    // intended purpose!!!
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(k);
Chris@16:    if(k == 0)
Chris@16:       return 1;
Chris@16:    T result = 1;
Chris@16:    for(int i = 0; i < k; ++i)
Chris@16:       result *= (a+i) / (b+i);
Chris@16:    return result;
Chris@16: }
Chris@16: //
Chris@16: // Routine for a > 15, b < 1
Chris@16: //
Chris@16: // Begin by figuring out how large our table of Pn's should be,
Chris@16: // quoted accuracies are "guestimates" based on empiracal observation.
Chris@16: // Note that the table size should never exceed the size of our
Chris@16: // tables of factorials.
Chris@16: //
Chris@16: template <class T>
Chris@16: struct Pn_size
Chris@16: {
Chris@16:    // This is likely to be enough for ~35-50 digit accuracy
Chris@16:    // but it's hard to quantify exactly:
Chris@16:    BOOST_STATIC_CONSTANT(unsigned, value = 50);
Chris@16:    BOOST_STATIC_ASSERT(::boost::math::max_factorial<T>::value >= 100);
Chris@16: };
Chris@16: template <>
Chris@16: struct Pn_size<float>
Chris@16: {
Chris@16:    BOOST_STATIC_CONSTANT(unsigned, value = 15); // ~8-15 digit accuracy
Chris@16:    BOOST_STATIC_ASSERT(::boost::math::max_factorial<float>::value >= 30);
Chris@16: };
Chris@16: template <>
Chris@16: struct Pn_size<double>
Chris@16: {
Chris@16:    BOOST_STATIC_CONSTANT(unsigned, value = 30); // 16-20 digit accuracy
Chris@16:    BOOST_STATIC_ASSERT(::boost::math::max_factorial<double>::value >= 60);
Chris@16: };
Chris@16: template <>
Chris@16: struct Pn_size<long double>
Chris@16: {
Chris@16:    BOOST_STATIC_CONSTANT(unsigned, value = 50); // ~35-50 digit accuracy
Chris@16:    BOOST_STATIC_ASSERT(::boost::math::max_factorial<long double>::value >= 100);
Chris@16: };
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const Policy& pol, bool normalised)
Chris@16: {
Chris@16:    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    //
Chris@16:    // This is DiDonato and Morris's BGRAT routine, see Eq's 9 through 9.6.
Chris@16:    //
Chris@16:    // Some values we'll need later, these are Eq 9.1:
Chris@16:    //
Chris@16:    T bm1 = b - 1;
Chris@16:    T t = a + bm1 / 2;
Chris@16:    T lx, u;
Chris@16:    if(y < 0.35)
Chris@16:       lx = boost::math::log1p(-y, pol);
Chris@16:    else
Chris@16:       lx = log(x);
Chris@16:    u = -t * lx;
Chris@16:    // and from from 9.2:
Chris@16:    T prefix;
Chris@16:    T h = regularised_gamma_prefix(b, u, pol, lanczos_type());
Chris@16:    if(h <= tools::min_value<T>())
Chris@16:       return s0;
Chris@16:    if(normalised)
Chris@16:    {
Chris@16:       prefix = h / boost::math::tgamma_delta_ratio(a, b, pol);
Chris@16:       prefix /= pow(t, b);
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       prefix = full_igamma_prefix(b, u, pol) / pow(t, b);
Chris@16:    }
Chris@16:    prefix *= mult;
Chris@16:    //
Chris@16:    // now we need the quantity Pn, unfortunatately this is computed
Chris@16:    // recursively, and requires a full history of all the previous values
Chris@16:    // so no choice but to declare a big table and hope it's big enough...
Chris@16:    //
Chris@16:    T p[ ::boost::math::detail::Pn_size<T>::value ] = { 1 };  // see 9.3.
Chris@16:    //
Chris@16:    // Now an initial value for J, see 9.6:
Chris@16:    //
Chris@16:    T j = boost::math::gamma_q(b, u, pol) / h;
Chris@16:    //
Chris@16:    // Now we can start to pull things together and evaluate the sum in Eq 9:
Chris@16:    //
Chris@16:    T sum = s0 + prefix * j;  // Value at N = 0
Chris@16:    // some variables we'll need:
Chris@16:    unsigned tnp1 = 1; // 2*N+1
Chris@16:    T lx2 = lx / 2;
Chris@16:    lx2 *= lx2;
Chris@16:    T lxp = 1;
Chris@16:    T t4 = 4 * t * t;
Chris@16:    T b2n = b;
Chris@16: 
Chris@16:    for(unsigned n = 1; n < sizeof(p)/sizeof(p[0]); ++n)
Chris@16:    {
Chris@16:       /*
Chris@16:       // debugging code, enable this if you want to determine whether
Chris@16:       // the table of Pn's is large enough...
Chris@16:       //
Chris@16:       static int max_count = 2;
Chris@16:       if(n > max_count)
Chris@16:       {
Chris@16:          max_count = n;
Chris@16:          std::cerr << "Max iterations in BGRAT was " << n << std::endl;
Chris@16:       }
Chris@16:       */
Chris@16:       //
Chris@16:       // begin by evaluating the next Pn from Eq 9.4:
Chris@16:       //
Chris@16:       tnp1 += 2;
Chris@16:       p[n] = 0;
Chris@16:       T mbn = b - n;
Chris@16:       unsigned tmp1 = 3;
Chris@16:       for(unsigned m = 1; m < n; ++m)
Chris@16:       {
Chris@16:          mbn = m * b - n;
Chris@16:          p[n] += mbn * p[n-m] / boost::math::unchecked_factorial<T>(tmp1);
Chris@16:          tmp1 += 2;
Chris@16:       }
Chris@16:       p[n] /= n;
Chris@16:       p[n] += bm1 / boost::math::unchecked_factorial<T>(tnp1);
Chris@16:       //
Chris@16:       // Now we want Jn from Jn-1 using Eq 9.6:
Chris@16:       //
Chris@16:       j = (b2n * (b2n + 1) * j + (u + b2n + 1) * lxp) / t4;
Chris@16:       lxp *= lx2;
Chris@16:       b2n += 2;
Chris@16:       //
Chris@16:       // pull it together with Eq 9:
Chris@16:       //
Chris@16:       T r = prefix * p[n] * j;
Chris@16:       sum += r;
Chris@16:       if(r > 1)
Chris@16:       {
Chris@16:          if(fabs(r) < fabs(tools::epsilon<T>() * sum))
Chris@16:             break;
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          if(fabs(r / tools::epsilon<T>()) < fabs(sum))
Chris@16:             break;
Chris@16:       }
Chris@16:    }
Chris@16:    return sum;
Chris@16: } // 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)
Chris@16: 
Chris@16: //
Chris@16: // For integer arguments we can relate the incomplete beta to the
Chris@16: // complement of the binomial distribution cdf and use this finite sum.
Chris@16: //
Chris@16: template <class T>
Chris@101: T binomial_ccdf(T n, T k, T x, T y)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // ADL of std names
Chris@101: 
Chris@16:    T result = pow(x, n);
Chris@101: 
Chris@101:    if(result > tools::min_value<T>())
Chris@16:    {
Chris@101:       T term = result;
Chris@101:       for(unsigned i = itrunc(T(n - 1)); i > k; --i)
Chris@101:       {
Chris@101:          term *= ((i + 1) * y) / ((n - i) * x);
Chris@101:          result += term;
Chris@101:       }
Chris@101:    }
Chris@101:    else
Chris@101:    {
Chris@101:       // First term underflows so we need to start at the mode of the
Chris@101:       // distribution and work outwards:
Chris@101:       int start = itrunc(n * x);
Chris@101:       if(start <= k + 1)
Chris@101:          start = itrunc(k + 2);
Chris@101:       result = pow(x, start) * pow(y, n - start) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(start));
Chris@101:       if(result == 0)
Chris@101:       {
Chris@101:          // OK, starting slightly above the mode didn't work, 
Chris@101:          // we'll have to sum the terms the old fashioned way:
Chris@101:          for(unsigned i = start - 1; i > k; --i)
Chris@101:          {
Chris@101:             result += pow(x, (int)i) * pow(y, n - i) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(i));
Chris@101:          }
Chris@101:       }
Chris@101:       else
Chris@101:       {
Chris@101:          T term = result;
Chris@101:          T start_term = result;
Chris@101:          for(unsigned i = start - 1; i > k; --i)
Chris@101:          {
Chris@101:             term *= ((i + 1) * y) / ((n - i) * x);
Chris@101:             result += term;
Chris@101:          }
Chris@101:          term = start_term;
Chris@101:          for(unsigned i = start + 1; i <= n; ++i)
Chris@101:          {
Chris@101:             term *= (n - i + 1) * x / (i * y);
Chris@101:             result += term;
Chris@101:          }
Chris@101:       }
Chris@16:    }
Chris@16: 
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: 
Chris@16: //
Chris@16: // The incomplete beta function implementation:
Chris@16: // This is just a big bunch of spagetti code to divide up the
Chris@16: // input range and select the right implementation method for
Chris@16: // each domain:
Chris@16: //
Chris@16: template <class T, class Policy>
Chris@16: T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_derivative)
Chris@16: {
Chris@16:    static const char* function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)";
Chris@16:    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
Chris@16:    BOOST_MATH_STD_USING // for ADL of std math functions.
Chris@16: 
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(a);
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(b);
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(x);
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(inv);
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(normalised);
Chris@16: 
Chris@16:    bool invert = inv;
Chris@16:    T fract;
Chris@16:    T y = 1 - x;
Chris@16: 
Chris@16:    BOOST_ASSERT((p_derivative == 0) || normalised);
Chris@16: 
Chris@16:    if(p_derivative)
Chris@16:       *p_derivative = -1; // value not set.
Chris@16: 
Chris@16:    if((x < 0) || (x > 1))
Chris@101:       return policies::raise_domain_error<T>(function, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x, pol);
Chris@16: 
Chris@16:    if(normalised)
Chris@16:    {
Chris@16:       if(a < 0)
Chris@101:          return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol);
Chris@16:       if(b < 0)
Chris@101:          return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol);
Chris@16:       // extend to a few very special cases:
Chris@16:       if(a == 0)
Chris@16:       {
Chris@16:          if(b == 0)
Chris@101:             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);
Chris@16:          if(b > 0)
Chris@16:             return inv ? 0 : 1;
Chris@16:       }
Chris@16:       else if(b == 0)
Chris@16:       {
Chris@16:          if(a > 0)
Chris@16:             return inv ? 1 : 0;
Chris@16:       }
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       if(a <= 0)
Chris@101:          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);
Chris@16:       if(b <= 0)
Chris@101:          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);
Chris@16:    }
Chris@16: 
Chris@16:    if(x == 0)
Chris@16:    {
Chris@16:       if(p_derivative)
Chris@16:       {
Chris@16:          *p_derivative = (a == 1) ? (T)1 : (a < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2);
Chris@16:       }
Chris@16:       return (invert ? (normalised ? T(1) : boost::math::beta(a, b, pol)) : T(0));
Chris@16:    }
Chris@16:    if(x == 1)
Chris@16:    {
Chris@16:       if(p_derivative)
Chris@16:       {
Chris@16:          *p_derivative = (b == 1) ? T(1) : (b < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2);
Chris@16:       }
Chris@16:       return (invert == 0 ? (normalised ? 1 : boost::math::beta(a, b, pol)) : 0);
Chris@16:    }
Chris@16:    if((a == 0.5f) && (b == 0.5f))
Chris@16:    {
Chris@16:       // We have an arcsine distribution:
Chris@16:       if(p_derivative)
Chris@16:       {
Chris@101:          *p_derivative = 1 / constants::pi<T>() * sqrt(y * x);
Chris@16:       }
Chris@16:       T p = invert ? asin(sqrt(y)) / constants::half_pi<T>() : asin(sqrt(x)) / constants::half_pi<T>();
Chris@16:       if(!normalised)
Chris@16:          p *= constants::pi<T>();
Chris@16:       return p;
Chris@16:    }
Chris@16:    if(a == 1)
Chris@16:    {
Chris@16:       std::swap(a, b);
Chris@16:       std::swap(x, y);
Chris@16:       invert = !invert;
Chris@16:    }
Chris@16:    if(b == 1)
Chris@16:    {
Chris@16:       //
Chris@16:       // Special case see: http://functions.wolfram.com/GammaBetaErf/BetaRegularized/03/01/01/
Chris@16:       //
Chris@16:       if(a == 1)
Chris@16:       {
Chris@16:          if(p_derivative)
Chris@101:             *p_derivative = 1;
Chris@16:          return invert ? y : x;
Chris@16:       }
Chris@16:       
Chris@16:       if(p_derivative)
Chris@16:       {
Chris@101:          *p_derivative = a * pow(x, a - 1);
Chris@16:       }
Chris@16:       T p;
Chris@16:       if(y < 0.5)
Chris@101:          p = invert ? T(-boost::math::expm1(a * boost::math::log1p(-y, pol), pol)) : T(exp(a * boost::math::log1p(-y, pol)));
Chris@16:       else
Chris@101:          p = invert ? T(-boost::math::powm1(x, a, pol)) : T(pow(x, a));
Chris@16:       if(!normalised)
Chris@16:          p /= a;
Chris@16:       return p;
Chris@16:    }
Chris@16: 
Chris@16:    if((std::min)(a, b) <= 1)
Chris@16:    {
Chris@16:       if(x > 0.5)
Chris@16:       {
Chris@16:          std::swap(a, b);
Chris@16:          std::swap(x, y);
Chris@16:          invert = !invert;
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(invert);
Chris@16:       }
Chris@16:       if((std::max)(a, b) <= 1)
Chris@16:       {
Chris@16:          // Both a,b < 1:
Chris@16:          if((a >= (std::min)(T(0.2), b)) || (pow(x, a) <= 0.9))
Chris@16:          {
Chris@16:             if(!invert)
Chris@16:             {
Chris@16:                fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
Chris@16:                BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:             }
Chris@16:             else
Chris@16:             {
Chris@16:                fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
Chris@16:                invert = false;
Chris@16:                fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
Chris@16:                BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:             }
Chris@16:          }
Chris@16:          else
Chris@16:          {
Chris@16:             std::swap(a, b);
Chris@16:             std::swap(x, y);
Chris@16:             invert = !invert;
Chris@16:             if(y >= 0.3)
Chris@16:             {
Chris@16:                if(!invert)
Chris@16:                {
Chris@16:                   fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
Chris@16:                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:                }
Chris@16:                else
Chris@16:                {
Chris@16:                   fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
Chris@16:                   invert = false;
Chris@16:                   fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
Chris@16:                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:                }
Chris@16:             }
Chris@16:             else
Chris@16:             {
Chris@16:                // Sidestep on a, and then use the series representation:
Chris@16:                T prefix;
Chris@16:                if(!normalised)
Chris@16:                {
Chris@16:                   prefix = rising_factorial_ratio(T(a+b), a, 20);
Chris@16:                }
Chris@16:                else
Chris@16:                {
Chris@16:                   prefix = 1;
Chris@16:                }
Chris@16:                fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative);
Chris@16:                if(!invert)
Chris@16:                {
Chris@16:                   fract = beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
Chris@16:                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:                }
Chris@16:                else
Chris@16:                {
Chris@16:                   fract -= (normalised ? 1 : boost::math::beta(a, b, pol));
Chris@16:                   invert = false;
Chris@16:                   fract = -beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
Chris@16:                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:                }
Chris@16:             }
Chris@16:          }
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          // One of a, b < 1 only:
Chris@16:          if((b <= 1) || ((x < 0.1) && (pow(b * x, a) <= 0.7)))
Chris@16:          {
Chris@16:             if(!invert)
Chris@16:             {
Chris@16:                fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
Chris@16:                BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:             }
Chris@16:             else
Chris@16:             {
Chris@16:                fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
Chris@16:                invert = false;
Chris@16:                fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
Chris@16:                BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:             }
Chris@16:          }
Chris@16:          else
Chris@16:          {
Chris@16:             std::swap(a, b);
Chris@16:             std::swap(x, y);
Chris@16:             invert = !invert;
Chris@16: 
Chris@16:             if(y >= 0.3)
Chris@16:             {
Chris@16:                if(!invert)
Chris@16:                {
Chris@16:                   fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
Chris@16:                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:                }
Chris@16:                else
Chris@16:                {
Chris@16:                   fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
Chris@16:                   invert = false;
Chris@16:                   fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
Chris@16:                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:                }
Chris@16:             }
Chris@16:             else if(a >= 15)
Chris@16:             {
Chris@16:                if(!invert)
Chris@16:                {
Chris@16:                   fract = beta_small_b_large_a_series(a, b, x, y, T(0), T(1), pol, normalised);
Chris@16:                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:                }
Chris@16:                else
Chris@16:                {
Chris@16:                   fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
Chris@16:                   invert = false;
Chris@16:                   fract = -beta_small_b_large_a_series(a, b, x, y, fract, T(1), pol, normalised);
Chris@16:                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:                }
Chris@16:             }
Chris@16:             else
Chris@16:             {
Chris@16:                // Sidestep to improve errors:
Chris@16:                T prefix;
Chris@16:                if(!normalised)
Chris@16:                {
Chris@16:                   prefix = rising_factorial_ratio(T(a+b), a, 20);
Chris@16:                }
Chris@16:                else
Chris@16:                {
Chris@16:                   prefix = 1;
Chris@16:                }
Chris@16:                fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative);
Chris@16:                BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:                if(!invert)
Chris@16:                {
Chris@16:                   fract = beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
Chris@16:                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:                }
Chris@16:                else
Chris@16:                {
Chris@16:                   fract -= (normalised ? 1 : boost::math::beta(a, b, pol));
Chris@16:                   invert = false;
Chris@16:                   fract = -beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
Chris@16:                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:                }
Chris@16:             }
Chris@16:          }
Chris@16:       }
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       // Both a,b >= 1:
Chris@16:       T lambda;
Chris@16:       if(a < b)
Chris@16:       {
Chris@16:          lambda = a - (a + b) * x;
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          lambda = (a + b) * y - b;
Chris@16:       }
Chris@16:       if(lambda < 0)
Chris@16:       {
Chris@16:          std::swap(a, b);
Chris@16:          std::swap(x, y);
Chris@16:          invert = !invert;
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(invert);
Chris@16:       }
Chris@16:       
Chris@16:       if(b < 40)
Chris@16:       {
Chris@16:          if((floor(a) == a) && (floor(b) == b) && (a < (std::numeric_limits<int>::max)() - 100))
Chris@16:          {
Chris@16:             // relate to the binomial distribution and use a finite sum:
Chris@16:             T k = a - 1;
Chris@16:             T n = b + k;
Chris@16:             fract = binomial_ccdf(n, k, x, y);
Chris@16:             if(!normalised)
Chris@16:                fract *= boost::math::beta(a, b, pol);
Chris@16:             BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:          }
Chris@16:          else if(b * x <= 0.7)
Chris@16:          {
Chris@16:             if(!invert)
Chris@16:             {
Chris@16:                fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
Chris@16:                BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:             }
Chris@16:             else
Chris@16:             {
Chris@16:                fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
Chris@16:                invert = false;
Chris@16:                fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
Chris@16:                BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:             }
Chris@16:          }
Chris@16:          else if(a > 15)
Chris@16:          {
Chris@16:             // sidestep so we can use the series representation:
Chris@16:             int n = itrunc(T(floor(b)), pol);
Chris@16:             if(n == b)
Chris@16:                --n;
Chris@16:             T bbar = b - n;
Chris@16:             T prefix;
Chris@16:             if(!normalised)
Chris@16:             {
Chris@16:                prefix = rising_factorial_ratio(T(a+bbar), bbar, n);
Chris@16:             }
Chris@16:             else
Chris@16:             {
Chris@16:                prefix = 1;
Chris@16:             }
Chris@16:             fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(0));
Chris@16:             fract = beta_small_b_large_a_series(a,  bbar, x, y, fract, T(1), pol, normalised);
Chris@16:             fract /= prefix;
Chris@16:             BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:          }
Chris@16:          else if(normalised)
Chris@16:          {
Chris@101:             // The formula here for the non-normalised case is tricky to figure
Chris@16:             // out (for me!!), and requires two pochhammer calculations rather
Chris@101:             // than one, so leave it for now and only use this in the normalized case....
Chris@16:             int n = itrunc(T(floor(b)), pol);
Chris@16:             T bbar = b - n;
Chris@16:             if(bbar <= 0)
Chris@16:             {
Chris@16:                --n;
Chris@16:                bbar += 1;
Chris@16:             }
Chris@16:             fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(0));
Chris@16:             fract += ibeta_a_step(a, bbar, x, y, 20, pol, normalised, static_cast<T*>(0));
Chris@16:             if(invert)
Chris@101:                fract -= 1;  // Note this line would need changing if we ever enable this branch in non-normalized case
Chris@16:             fract = beta_small_b_large_a_series(T(a+20),  bbar, x, y, fract, T(1), pol, normalised);
Chris@16:             if(invert)
Chris@16:             {
Chris@16:                fract = -fract;
Chris@16:                invert = false;
Chris@16:             }
Chris@16:             BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:          }
Chris@16:          else
Chris@16:          {
Chris@16:             fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
Chris@16:             BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:          }
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(fract);
Chris@16:       }
Chris@16:    }
Chris@16:    if(p_derivative)
Chris@16:    {
Chris@16:       if(*p_derivative < 0)
Chris@16:       {
Chris@16:          *p_derivative = ibeta_power_terms(a, b, x, y, lanczos_type(), true, pol);
Chris@16:       }
Chris@16:       T div = y * x;
Chris@16: 
Chris@16:       if(*p_derivative != 0)
Chris@16:       {
Chris@16:          if((tools::max_value<T>() * div < *p_derivative))
Chris@16:          {
Chris@16:             // overflow, return an arbitarily large value:
Chris@16:             *p_derivative = tools::max_value<T>() / 2;
Chris@16:          }
Chris@16:          else
Chris@16:          {
Chris@16:             *p_derivative /= div;
Chris@16:          }
Chris@16:       }
Chris@16:    }
Chris@16:    return invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) - fract : fract;
Chris@16: } // template <class T, class Lanczos>T ibeta_imp(T a, T b, T x, const Lanczos& l, bool inv, bool normalised)
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised)
Chris@16: {
Chris@16:    return ibeta_imp(a, b, x, pol, inv, normalised, static_cast<T*>(0));
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: T ibeta_derivative_imp(T a, T b, T x, const Policy& pol)
Chris@16: {
Chris@16:    static const char* function = "ibeta_derivative<%1%>(%1%,%1%,%1%)";
Chris@16:    //
Chris@16:    // start with the usual error checks:
Chris@16:    //
Chris@16:    if(a <= 0)
Chris@101:       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);
Chris@16:    if(b <= 0)
Chris@101:       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);
Chris@16:    if((x < 0) || (x > 1))
Chris@101:       return policies::raise_domain_error<T>(function, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x, pol);
Chris@16:    //
Chris@16:    // Now the corner cases:
Chris@16:    //
Chris@16:    if(x == 0)
Chris@16:    {
Chris@16:       return (a > 1) ? 0 : 
Chris@16:          (a == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error<T>(function, 0, pol);
Chris@16:    }
Chris@16:    else if(x == 1)
Chris@16:    {
Chris@16:       return (b > 1) ? 0 :
Chris@16:          (b == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error<T>(function, 0, pol);
Chris@16:    }
Chris@16:    //
Chris@16:    // Now the regular cases:
Chris@16:    //
Chris@16:    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
Chris@16:    T f1 = ibeta_power_terms<T>(a, b, x, 1 - x, lanczos_type(), true, pol);
Chris@16:    T y = (1 - x) * x;
Chris@16: 
Chris@16:    if(f1 == 0)
Chris@16:       return 0;
Chris@16:    
Chris@16:    if((tools::max_value<T>() * y < f1))
Chris@16:    {
Chris@16:       // overflow:
Chris@16:       return policies::raise_overflow_error<T>(function, 0, pol);
Chris@16:    }
Chris@16: 
Chris@16:    f1 /= y;
Chris@16: 
Chris@16:    return f1;
Chris@16: }
Chris@16: //
Chris@16: // Some forwarding functions that dis-ambiguate the third argument type:
Chris@16: //
Chris@16: template <class RT1, class RT2, class Policy>
Chris@16: inline typename tools::promote_args<RT1, RT2>::type 
Chris@16:    beta(RT1 a, RT2 b, const Policy&, const mpl::true_*)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename tools::promote_args<RT1, RT2>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    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%)");
Chris@16: }
Chris@16: template <class RT1, class RT2, class RT3>
Chris@16: inline typename tools::promote_args<RT1, RT2, RT3>::type 
Chris@16:    beta(RT1 a, RT2 b, RT3 x, const mpl::false_*)
Chris@16: {
Chris@16:    return boost::math::beta(a, b, x, policies::policy<>());
Chris@16: }
Chris@16: } // namespace detail
Chris@16: 
Chris@16: //
Chris@16: // The actual function entry-points now follow, these just figure out
Chris@16: // which Lanczos approximation to use
Chris@16: // and forward to the implementation functions:
Chris@16: //
Chris@16: template <class RT1, class RT2, class A>
Chris@16: inline typename tools::promote_args<RT1, RT2, A>::type 
Chris@16:    beta(RT1 a, RT2 b, A arg)
Chris@16: {
Chris@16:    typedef typename policies::is_policy<A>::type tag;
Chris@16:    return boost::math::detail::beta(a, b, arg, static_cast<tag*>(0));
Chris@16: }
Chris@16: 
Chris@16: template <class RT1, class RT2>
Chris@16: inline typename tools::promote_args<RT1, RT2>::type 
Chris@16:    beta(RT1 a, RT2 b)
Chris@16: {
Chris@16:    return boost::math::beta(a, b, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class RT1, class RT2, class RT3, class Policy>
Chris@16: inline typename tools::promote_args<RT1, RT2, RT3>::type 
Chris@16:    beta(RT1 a, RT2 b, RT3 x, const Policy&)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    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%)");
Chris@16: }
Chris@16: 
Chris@16: template <class RT1, class RT2, class RT3, class Policy>
Chris@16: inline typename tools::promote_args<RT1, RT2, RT3>::type 
Chris@16:    betac(RT1 a, RT2 b, RT3 x, const Policy&)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    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%)");
Chris@16: }
Chris@16: template <class RT1, class RT2, class RT3>
Chris@16: inline typename tools::promote_args<RT1, RT2, RT3>::type 
Chris@16:    betac(RT1 a, RT2 b, RT3 x)
Chris@16: {
Chris@16:    return boost::math::betac(a, b, x, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class RT1, class RT2, class RT3, class Policy>
Chris@16: inline typename tools::promote_args<RT1, RT2, RT3>::type 
Chris@16:    ibeta(RT1 a, RT2 b, RT3 x, const Policy&)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    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%)");
Chris@16: }
Chris@16: template <class RT1, class RT2, class RT3>
Chris@16: inline typename tools::promote_args<RT1, RT2, RT3>::type 
Chris@16:    ibeta(RT1 a, RT2 b, RT3 x)
Chris@16: {
Chris@16:    return boost::math::ibeta(a, b, x, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class RT1, class RT2, class RT3, class Policy>
Chris@16: inline typename tools::promote_args<RT1, RT2, RT3>::type 
Chris@16:    ibetac(RT1 a, RT2 b, RT3 x, const Policy&)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    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%)");
Chris@16: }
Chris@16: template <class RT1, class RT2, class RT3>
Chris@16: inline typename tools::promote_args<RT1, RT2, RT3>::type 
Chris@16:    ibetac(RT1 a, RT2 b, RT3 x)
Chris@16: {
Chris@16:    return boost::math::ibetac(a, b, x, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class RT1, class RT2, class RT3, class Policy>
Chris@16: inline typename tools::promote_args<RT1, RT2, RT3>::type 
Chris@16:    ibeta_derivative(RT1 a, RT2 b, RT3 x, const Policy&)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    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%)");
Chris@16: }
Chris@16: template <class RT1, class RT2, class RT3>
Chris@16: inline typename tools::promote_args<RT1, RT2, RT3>::type 
Chris@16:    ibeta_derivative(RT1 a, RT2 b, RT3 x)
Chris@16: {
Chris@16:    return boost::math::ibeta_derivative(a, b, x, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: } // namespace math
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #include <boost/math/special_functions/detail/ibeta_inverse.hpp>
Chris@16: #include <boost/math/special_functions/detail/ibeta_inv_ab.hpp>
Chris@16: 
Chris@16: #endif // BOOST_MATH_SPECIAL_BETA_HPP
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: